Control software analysis, Part I
Open-loop properties

We consider a standard control system, in this case a one-mass, one spring system. such as that shown in Fig. 1.

The system is modeled via the following differential equation

$$\begin{array}[]{rcl}\displaystyle\frac{d}{dt}\left[\begin{array}[]{c}x_{p}\\ \dot{x}_{p}\end{array}\right]&=&\left[\begin{array}[]{cc}0&1\\ -1&0\end{array}\right]\left[\begin{array}[]{c}x_{p}\\ \dot{x}_{p}\end{array}\right]+\left[\begin{array}[]{c}0\\ 1\end{array}\right]u,\;\;x(0)=x_{0},\dot{x}_{p}(0)=\dot{x}_{p,0}\\[10.0pt] y&=&\left[0\;\;1\right]\left[\begin{array}[]{c}x_{p}\\ \dot{x}_{p}\end{array}\right].\end{array}$$

(1)

In these equations, $x_{p}$ and $\dot{x}_{p}$ are measured in meters and meters per second, respectively, while $u$ is expressed in Newtons. We assume that the output $y$ is available for feedback purposes. The chosen controller is of the lead-lag type,

$$\begin{array}[]{l}\tilde{y}(t)={\bf SAT}(y(t)),\\ u(s)=\displaystyle-128\frac{s+1}{s+0.1}\frac{s/5+1}{s/50+1}\tilde{y}(s),\end{array}$$

(2)

which guarantees good closed-loop system bandwidth (in excess of 20 rad/sec) and good steady-state disturbance rejection properties; following common practices, the input to the controller $y$ is first passed through a saturation function to avoid variable overflow in the controller. The essential properties of the closed-loop system may be easily extracted from the Bode plot of the open-loop transfer function of the plant mounted in series with the controller, shown in Fig. 2.

For example, the phase margin can be easily seen to exceed 50 degrees. This controller is then implemented on a computer by considering a sample-and-hold architecture running at 100 Hz, by implementing the state-space model

$$\begin{array}[]{rcl}x_{k+1}&=&\left[\begin{array}[]{cc}0.499&-0.050\\ 0.010&1.000\end{array}\right]x_{k}+\left[\begin{array}[]{c}1\\ 0\end{array}\right]{\bf SAT}(y_{k})\\ &=&Ax_{k}+B{\bf SAT}(y_{k})\\[10.0pt] u_{k}&=&\left[564.48\;\;0\right]x_{k}-1280\>{\bf SAT}(y_{k})\\ &=&Cx_{k}+D{\bf SAT}(y_{k})\end{array}$$

(3)

which is no more than a first-order Euler discretization of the following state-space realization of the controller (2)

$$\begin{array}[]{rcl}\displaystyle\frac{d}{dt}x&=&\left[\begin{array}[]{cc}-50.1&-5.0\\ 1.0&0.0\end{array}\right]x+\left[\begin{array}[]{c}100\\ 0\end{array}\right]{\bf SAT}(y)\\[10.0pt] u&=&\left[564.48\;\;0\right]x-1280\>{\bf SAT}(y).\end{array}$$

In the expression for the compensator, $y_{k}$ is $y(0.01k)$, the output of the system sampled at time $0.01k$. Likewise, $u_{k}$ is the compensator output fed to the plant during the time interval $\left[0.01k\;\;0.01(k+1)\right)$, resulting in the closed-loop description of the system as shown in Fig. 3.

This figure includes the sampling process ${\cal S}$ and the zero-order hold ZOH. An implementation of the compensator may then be realized by means of the simplified, but functional C code shown in Fig. 4.

1: int main(int argc, char *argv[])
2: {
3:   double *x[2], *xb[2], y, u;
4:   x[0] = 0;
5:   x[1] = 0;
6:   while(1){
7:     fscanf(stdin,"%f",&y);
8:     if (y >1){
9:       y=1
10:Ψ   }
11:Ψ   if (y<-1){
12:      y=-1
13:Ψ   }
14:    u = 564.48*x[0]-1280*y;
15:    fprintf(stdout,"%f\n",u);
16:    xb[0]=x[0];
17:    xb[1]=x[1];
18:    x[0]:= 0.4990*xb[0]-0.0500*xb[1]+y;
19:    x[1]:= 0.01*xb[0]+xb[1];
20:  }
21:}

The final step involves the compilation of the C code to executable code, and the eventual operation of the system in closed-loop.

Considering the entire system, made of the physical plant and its controller, we are led to asking whether “the system works”. In this paper, we are interested in asking this question about the software implementation of the control law shown in Fig. 4 and the corresponding executable binary code.

Many control engineers and researchers may wonder whether this question makes any sense at all: After all, control system textbooks suggest that the task of verifying whether the closed-loop system “works” amounts to establishing, in some way or some other, that the closed-loop “system” that consists of the system (1) combined with either the continuous-time controller (2) or the discrete-time controller (3) is stable, which requires no more than computing the roots of a low-order polynomial. Further, control systems manufacturers may require the controller itself to be open-loop stable. Other questions about system performance, including disturbance rejection properties and robustness against modeling inaccuracies, may also be easily answered by using algebraic or graphical methods.

Note, however, that the process of generating C code from the continuous or discrete control laws is often done automatically, by means of commercial autocoders whose implementation details are not fully known. Likewise, proprietary or public-domain compilers are then used to produce binary executables from the C code. While autocoders and compilers can be assumed to operate in good faith, they cannot just be trusted blindly. In addition, the codes generated out of the specifications often contain detailed computations whose details are not known at the level of control law specifications. Thus the produced C code and the binary executable must be examined in detail for adequate performance by following a process that is independent from the original code generation process. The amount of work involved in such an endeavour may appear formidable. However, the alternative which is to certify the correctness of autocoders and compilers is even more daunting [Rin99, AR04, RM99].

The rest of this paper is devoted to bringing forward the elements that make the static analysis of code that implements control system specification not only possible, but also necessary given how easy it is to perform.

Analyzing the software that implements a control system consists of answering several questions:

• •

  Does the control function keep executing properly and within acceptable execution time limits?

• •

  Is the implemented controller stable and what is the range of values explored by each variable?

• •

  Is the closed-loop system stable?

• •

  Does the software implement the specifications?

These questions may be extended to include further concerns about system performance. However, the problem remains fundamentally the same: That is, to reach the same level of confidence about system stability and performance when considering the controller code in Fig. 4 together with the plant (1) as that reached when analyzing the controller specification (3), together with the plant (1). To reach that goal, two options exist:

The first option is to follow a “macroscopic pattern matching” approach, whereby the controller specification (3) is reconstructed from the code in Fig. 4. This solution, adopted in [CCF⁺05, Fer04], is particularly advantageous when the code structure is relatively repeatable from one instance to the other. However, pattern matching can become more challenging if codes that implement controllers do not have well-identified and rigid structures.

The second option, which we call code-level analysis, is to develop proofs that directly identify the code semantics. Code-level analysis is, in principle, supported by an elegant logical environment that closely matches available control oriented methodologies. This environment, developed by Floyd and Hoare in the late 1960s [Flo67, Hoa69], intends to describe the code semantics by concatenating elementary knowledge about each line of code [Pel01, Ch.7]: Each line of code is annotated with one precondition pre satisfied by the program state prior to the line execution, and one postcondition post satisfied by the program state after the execution of the line. Thus each line of code loc is transformed into a Hoare triple of the kind {pre}loc{post }, where the pre and post are compatible with loc. Global program properties can then be inferred by concatenating Hoare triples, as follows: Considering two consecutive Hoare triples {pre_1}loc_1{post_1 } and {pre_2}loc_2{post_2 } and the property {post_1} $\rightarrow$ {pre_2}, whe have {pre_1}loc_1, loc_2{post_2 }. Further developments of Hoare logic will be used thereafter.

We now examine how Hoare logic can be used to formalize the analysis of control system stability requirements for the control system specifications and their implementation. Informally speaking, we know that the system that implements the controller is stable and all states are bounded because the matrix

$$A=\left[\begin{array}[]{cc}0.499&-0.050\\ 0.010&1.000\end{array}\right]$$

is stable, that is, its two eigenvalues have modulus less than one, and the input is effectively bounded in absolute value by one. While exact, this proof does not carry over easily to the software implementation. Moreover, it is somewhat incomplete since it does not tell anything related to the range covered by the variables involved. These issues may be addressed by means of Lyapunov stability theory: Consider the Lyapunov function $V(x)=x^{T}Px$, where

$${P}=\left[\begin{array}[]{cc}0.0300&0.2000\\ 0.2000&10.0000\end{array}\right]$$

and $P$ is a positive definite matrix. Then it is easy to show that the value of $x^{T}Px$ decays along the trajectories of

$$x_{k+1}=Ax_{k},$$

that is $V(x_{k}+1)-V(x_{k})<0$ if $y_{k}=0$ and $x_{k}\neq 0$. That alone indicates stability of the controller, yet it does not bring any indication about the size of the states of the system. For that purpose, and following well-known practices [BEFB94], it is possible to show that $P$ defines not only a valid Lyapunov function, but also an invariant ellipsoid

$${\cal E}_{P}=\left\{x\in{\bf R}^{2}\;|\;x^{T}Px\leq 1\right\}.$$

for any input $y$, and this ellipsoid defines a bounded region for $x_{k}$. In particular, if $x_{0}$ is included in ${\cal E}_{{P}}$ (and it does since $x_{0}=0$), then the state will stay in the ellipsoid forever. A simple eigenvalue computation on $P$ shows that each entry of the state $x$ is bounded in absolute value by 7.

Let us now examine how this knowledge can be included in the controller’s executable code so as to establish the code’s correctness. We may begin our effort with analyzing a Matlab script that implements the controller, given in Fig. 5. The reason for doing so is that a Matlab code structure is sufficiently close to typical control system specifications to make the correspondence between traditional proofs and code proof transparent, while remaining executable by Matlab’s interpreter.

1:  A = [0.4990, -0.0500; 0.0100, 1.0000];
2:  C = [564.48, 0];
3:  B = [1;0];D = -1280;
4:  x = zeros(2,1);
5:  while 1
6:   y = fscanf(stdin,"%f");
7:   y = max(min(y,1),-1);
8:   u = C*x + D*y;
9:   fprintf(stdout,"%f\n",u);
10:   x = A*x + B*y;
11: end

Identifying $x$ as the controller’s state, we now proceed with establishing the stability of the controller’s Matlab implementation by using Hoare logic. For that, we may begin the program annotation at any line of choice, for example at line 8: u = C*x + Dy. This assignment does not influence the state $x$ but determines the value of $u$. It also introduces the controller input $y$. Thus we propose the following triple:

$$\begin{array}[]{l}\left\{x\in{\cal E}_{{P}},\;y^{2}\leq 1\right\}\\ \mbox{\tt 8: u = C*x+D*y;}\\ \left\{x\in{\cal E}_{{P}},\;y^{2}\leq 1,\;u^{2}\leq 2(C{P}^{-1}C^{T}+D^{2})\right\}\end{array}$$

(4)

Much about this statement may appear to be rather arbitrary at first sight, and indeed requires global knowledge of the program (or its specifications, described earlier) to be established. However, the statement $\left\{u^{2}\leq 2(C{P}^{-1}C^{T}+D^{2})\right\}$ can be easily figured out from the value of $u$ and the statements $\left\{x\in{\cal E}_{{P}},\;y^{2}\leq 1\right\}$. From there on, it is possible to obtain statements about the evolution of the state-space by means of backwards propagation [Pel01]. To begin with, line 7: y = max(min(y,1),-1) implicitly indicates that $y^{2}\leq 1$ after this assignment, and can be commented via the triple

{$x\in{\cal E}_{{P}}$}

7: y = max(min(y,1),-1);

$\left\{x\in{\cal E}_{{P}},\;y^{2}\leq 1\right\}$

At this point, the assignment in line 6: y = fscanf(stdin,"%f") can simply be surrounded by the same statement $\left\{x\in{\cal E}_{{P}}\right\}$ to yield the triple

$\{x\in{\cal E}_{{P}}\}$

6: y = fscanf(stdin,"%f")

$\{x\in{\cal E}_{{P}}\}$

Line 5: while 1 may not be separated from line 11: end. Following standard Hoare logic [Pel01, Mon03], these two lines may be documented as follows:

$\{x\in{\cal E}_{P}\}$

5: while 1

$\{x\in{\cal E}_{P}\}$

⋮

$\{x\in{\cal E}_{P}\}$

11: end

$\{false\}$

Moving backwards, we see that the two lines that must be worked on are lines 4 and 10. Considering first the assignment 4: x = zeros(2,1), we are led to the triple

$\{true\}$

4: x = zeros(2,1)

$\{x\in{\cal E}_{P}\}$.

Moving further backwards yields no changes concerning the state variable $x$, and triples may be constructed for lines 1 through 3 according to the template

$\{true\}$

loc

$\{true\}$.

Line 10: x = A*x + B*y introduces the controller input $y$ and performs an important transformation on the state $x$. Using insight from line 7: y = max(min(y,1),-1) we introduce the condition $y^{2}\leq 1$. We also propagate the statement $\left\{x\in{\cal E}_{{P}}\right\}$ backwards: Replacing $x$ by its value $Ax+By$ yields the statement $\left\{Ax+By\in{\cal E}_{P}\right\}$. Thus the triple

$\left\{Ax+By\in{\cal E}_{{P}},\;y^{2}\leq 1\right\}$

10: x = A*x + B*y;

$\left\{x\in{\cal E}_{{P}}\right\}$

Moving to line 9: fprintf(stdout,"%f$\backslash$n",u), we see that it really does not act in any way on the state of the program. However, it introduces the controller output value $u$, whose numerical range must be known. At this point, we postulate the property $u^{2}\leq 2(CP^{-1}C^{T}+D^{2})$ and we are led to the triple

$\{Ax+By\in{\cal E}_{P},\;u^{2}\leq 2(CP^{-1}C^{T}+D^{2}),\;y^{2}\leq 1\}$

9: fprintf(stdout,"%f$\backslash$n",u)

$\{Ax+By\in{\cal E}_{P},\;y^{2}\leq 1\}$.

For this triple to be “compatible” with the triple (4) associated with line 8, we must prove the statement

$$\{x\in{\cal E}_{P},\;y^{2}\leq 1\}\Rightarrow\left\{Ax+By\in{\cal E}_{{P}},\;y^{2}\leq 1\right\}.$$

(5)

For the specific instances of $A$, $B$ and $P$ under consideration, this statement is indeed true, yet not trivial. A simple proof relies on the ${\cal S}$-procedure [AG64], which amounts to introducing the auxiliary (and true) statement

$$\forall(x,y)\;(Ax+By)^{T}P(Ax+By)-0.01x^{T}Px-0.99y^{2}\leq 0.$$

(6)

leading to

$$\{x\in{\cal E}_{P},\;y^{2}\leq 1,\;(Ax+By)^{T}P(Ax+By)-0.01x^{T}Px-0.99y^{2}\leq 0\}\Rightarrow\left\{Ax+By\in{\cal E}_{{P}},\;y^{2}\leq 1\right\},$$

which may now be checked much more easily than (5) alone. Thus a preferrable triple for line 8 is

$\left\{x\in{\cal E}_{{P}},\;y^{2}\leq 1\right\}$
8: u = C*x+D*y;
$\left\{x\in{\cal E}_{{P}},\;y^{2}\leq 1,\right.$
$\left.(Ax+By)^{T}P(Ax+By)-0.01x^{T}Px-0.99y^{2}\leq 0,\;u^{2}\leq 2(C{P}^{-1}C^{T}+D^{2})\right\}$

Concatenating statements about all individual instructions together leads us to the following annotated Matlab code:

$\{true\}$

1: A = [0.4990, -0.0500; 0.0100, 1.0000];

$\{true\}$

2: C = [-564.48, 0];

$\{true\}$

3: B = [1;0];D=1280

$\{true\}$

4: x = zeros(2,1);

$\left\{x\in{\cal E}_{{P}}\right\}$

5: while 1

$\{x\in{\cal E}_{{P}}\}$

6: y = fscanf(stdin,"%f")

$\{x\in{\cal E}_{{P}}\}$

7: y = max(min(y,1),-1);

$\left\{x\in{\cal E}_{{P}},\;y^{2}\leq 1\right\}$

8: u = C*x+D*y;

$\left\{x\in{\cal E}_{{P}},\;y^{2}\leq 1,\right.$

$\left.(Ax+By)^{T}P(Ax+By)-0.01x^{T}Px-0.99y^{2}\leq 0,\;u^{2}\leq 2(C{P}^{-1}C^{T}+D^{2})\right\}$

skip

$\{Ax+By\in{\cal E}_{P},\;u^{2}\leq 2(CP^{-1}C^{T}+D^{2}),\;y^{2}\leq 1\}$

9: fprintf(stdout,"%f$\backslash$n",u)

$\left\{Ax+By\in{\cal E}_{{P}},\;y^{2}\leq 1\right\}$

10: x = A*x + B*y;

$\left\{x\in{\cal E}_{{P}}\right\}$

11: end

$\{false\}$

This annotated program contains the instruction $skip$ whose purpose is only to separate two consecutive logical statements.

We will perform only the forward analysis for the sake of brevity. Skipping variable declarations, the first relevant line of code is 4: x[0]=0. We may therefore write the triple

$\left\{true\right\}$

4: x[0] = 0

$\{x[0]=0\}$

Line ${\tt 5:x[1]=0}$ completes the initialization of $x=\left[\mbox{x[0] x[1]}\right]^{T}$, leading to the triple

$\{x[0]=0\}$

5: x[1] = 0

$\{x\in{\cal E}_{P}\}$.

Doing as precedently, line 6: while(1){, together with line 20: } leads to the set of commented lines

$\{x\in{\cal E}_{P}\}$

6: while(1){

$\{x\in{\cal E}_{P}\}$

$\vdots$

$\{x\in{\cal E}_{P}\}$

20: }

$\{false\}$

Line 7: fscanf(stdin,"$\backslash$%f",&y); likewise leads to the triple

$\{x\in{\cal E}_{P}\}$

7: fscanf(stdin,"$\backslash$%f",&y);

$\{x\in{\cal E}_{P}\}$.

The test 8: if (y >1){ $\ldots$ 10: } may be commented using insight from line 9 as

$\{x\in{\cal E}_{P}\}$

8: if (y >1){

$\{x\in{\cal E}_{P},\;y>1\}$

⋮

$\{x\in{\cal E}_{P},\;y=1\}$

10: }

$\{x\in{\cal E}_{P},\;y<=1\}$.

which is of course compatible with the assignment 9: y = 1

$\{x\in{\cal E}_{P},\;y>1\}$

9: y = 1

$\{x\in{\cal E}_{P},\;y=1\}$.

The second test 11: if (y<-1){ $\ldots$ 13: } is documented the same way

$\{x\in{\cal E}_{P},\;y<=1\}$

11: if (y<-1){

$\{x\in{\cal E}_{P},\;y<-1\}$

⋮

$\{x\in{\cal E}_{P},\;y=-1\}$

13: }

$\{x\in{\cal E}_{P},\;y^{2}\leq 1$

and line 12: y=-1 yields the triple

$\{x\in{\cal E}_{P},\;y<-1\}$

12: y=-1

$\{x\in{\cal E}_{P},\;y=-1\}$.

The latter postcondition is clearly the same as that of line 7 of the Matlab program shown in Fig. 5. Following the same process as for the Matlab program, the precondition for line 14: u = 564.48*x[0]-1280*y is chosen to be weaker than the post-condition of line 13; we choose the precondition

$$\left\{\left[\begin{array}[]{l}x\\ y\end{array}\right]\in{\cal E}_{Q}\right\},\;\;\mbox{ with }Q=\left[\begin{array}[]{cc}0.01P&0\\ 0&0.99\end{array}\right],$$

to generate the triple

$\left\{\left[\begin{array}[]{l}x\\ y\end{array}\right]\in{\cal E}_{Q}\right\}$

14: u = 564.48*x[0]-1280*y

$\left\{\left[\begin{array}[]{l}x\\ y\end{array}\right]\in{\cal E}_{Q},\;u^{2}\leq\left[C\;D\right]Q^{-1}\left[C\;D\right]^{T}\right\}$

Line 15: fprintf(stdout,"%f$\backslash$n",u); may in turn be instrumented as the triple

$\left\{\left[\begin{array}[]{l}x\\ y\end{array}\right]\in{\cal E}_{Q},\;u^{2}\leq\left[C\;D\right]Q^{-1}\left[C\;D\right]^{T}\right\}$

15: fprintf(stdout,”%f$\backslash$n”,u);

$\left\{\left[\begin{array}[]{l}x\\ y\end{array}\right]\in{\cal E}_{Q}\right\}$

Lines 16 and 17 utilise a buffer variable aimed at performing the controller state update. The effect of these lines is easier to figure out by means of matrix operations: For example, line 16: xb[0]=x[0] can be seen as the matrix operation

$$\left[\begin{array}[]{cc}x\\ y\\ \mbox{xb}[1]\end{array}\right]=\left[\begin{array}[]{ccc}1&0&0\\ 0&1&0\\ 0&0&1\\ 1&0&0\end{array}\right]\left[\begin{array}[]{c}x\\ y\end{array}\right]$$

Given this matrix operation, it is easy to figure out the image of the ellipsoid ${\cal E}_{Q}$ through this mapping. The parameterization of the resulting ellipsoid needs changing, however, because it is singular (it is flat along one of its semi-axes). Whenever the matrix $R$ is invertible, the following define the same ellipsoid:

$$\left\{x\;|\;x^{T}R^{-1}x\leq 1\right\}=\left\{x\;|\;\left[\begin{array}[]{cc}1&x^{T}\\ x&R\end{array}\right]\geq 0\right\}$$

When $R$ is not invertible, then only one of the statements above defines an ellipsoid. For the sake of conciseness, we will name ${\cal G}_{R}$ the ellipsoid defined as

$${\cal G}_{R}=\left\{x\;|\;\left[\begin{array}[]{cc}1&x^{T}\\ x&R\end{array}\right]\geq 0\right\}.$$

With such an alternative ellipsoid definition we may write a concise triple for line 16 as

$\left\{\left[\begin{array}[]{l}x\\ y\end{array}\right]\in{\cal E}_{Q}\right\}$

16: xb[0]=x[0]

$\left\{\left[\begin{array}[]{l}x\\ y\\ \mbox{xb}[1]\end{array}\right]\in{\cal G}_{R}\right\}$

with

$$R=\left[\begin{array}[]{ccc}1&0&0\\ 0&1&0\\ 0&0&1\\ 1&0&0\end{array}\right]Q^{-1}\left[\begin{array}[]{ccc}1&0&0\\ 0&1&0\\ 0&0&1\\ 1&0&0\end{array}\right]^{T}$$

Line 17: xb[1]=x[1]; may be treated likewise, and yield the triple

$\left\{\left[\begin{array}[]{l}x\\ y\\ \mbox{xb}[1]\end{array}\right]\in{\cal G}_{R}\right\}$

17: xb[1]=x[1];

$\left\{\left[\begin{array}[]{l}x\\ y\\ xb\end{array}\right]\in{\cal G}_{S}\right\}$

where

$$S=\left[\begin{array}[]{cccc}1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1\\ 0&1&0&0\end{array}\right]R\left[\begin{array}[]{cccc}1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1\\ 1&0&0&0\end{array}\right]^{T}\mbox{ and }xb=\left[\begin{array}[]{l}\mbox{xb}[1]\\ \mbox{xb}[2]\end{array}\right].$$

Continuing along the same lines the triple for 18: x[0]:= 0.4990*xb[0]-0.0500*xb[1]+y; is

$\left\{\left[\begin{array}[]{l}x\\ y\\ xb\end{array}\right]\in{\cal G}_{S}\right\}$

18: x[0]:= 0.4990*xb[0]-0.0500*xb[1]+y;

$\left\{\left[\begin{array}[]{l}x\\ y\\ xb\end{array}\right]\in{\cal G}_{T}\right\}$

with

$$T=\left[\begin{array}[]{ccccc}0&0&1&0.4990&-0.0500\\ 0&1&0&0&0\\ 0&0&1&0&0\\ 0&0&0&1&0\\ 0&0&0&0&1\end{array}\right]S\left[\begin{array}[]{ccccc}0&0&1&0.4990&-0.0500\\ 0&1&0&0&0\\ 0&0&1&0&0\\ 0&0&0&1&0\\ 0&0&0&0&1\end{array}\right]^{T},$$

and the triple for line 19: x[1]:= 0.01*xb[0]+xb[1]; is

$\left\{\left[\begin{array}[]{l}x\\ y\\ xb\end{array}\right]\in{\cal G}_{T}\right\}$

19: x[1]:= 0.01*xb[0]+xb[1];

$\left\{x\in{\cal G}_{U}\right\}$

with

$$U=\left[\begin{array}[]{ccccc}1&0&0&0&0\\ 0&0&0&0.01&1\end{array}\right]T\left[\begin{array}[]{ccccc}1&0&0&0&0\\ 0&0&0&0.01&1\end{array}\right]^{T}.$$

Explicit calculations then show $U^{-1}\geq P$, which ensures consistency between the postcondition of line 19 and the precondition of line 20. Moreover, it is easy to show that $U=\tilde{P}^{-1}$, where $\tilde{P}$ is defined in Eq. (7).

Bringing the forward analyses of the Matlab code, on the one hand, and of the C code, on the other hand, shows how the higher level-of-detail present in the C code instructions yields considerably more detailed information about the behavior of each variable. This hgher level of detail is precisely what makes the analysis of the C code necessary beyond that of the Matlab code. The C code analysis is, however, not complete and still contains ambiguities. Consider for example Line 18: x[0]:= 0.4990*xb[0]-0.0500*xb[1]+y;. Depending on the compiler used, this line may be executed as x[0]:= (0.4990*xb[0]+(-0.0500*xb[1]+y)); or x[0]:=((0.4990*xb[0]-0.0500*xb[1])+y);, for example. The resulting intermediate values may differ widely, yet must be computed in advance to offer a complete analysis of the real time computations that take place during comntroller operations. This justifies a further, more detailed analysis of the code that is produced by compiling that given in Fig. 4.

In addition, it is clear that the analyses shown on a simple second-order controller example easily extend to linear controllers of arbitrary order (possibly subject to input saturation).