Controllability:

1. General result:

A linear continuous-time time-invariant system:

is controllable iff: for any initial value

of the state at the initial time

and for a given final time

, there exists a control signal trajectory

such that

Let n be the system order (

). A necessary and sufficient condition for the system to be controllable (condition

) is:

𝒞 is called the controllability matrix .

2. Proof ():

The general solution of the state equation is:

Applied to our problem (

), it comes:

or:

3. Proof of the necessary condition ():

Cayley-Hamilton theorem (see: CayleyHamilton) allows to express:

Considering, without loss of generality a single input system (

), equation (1) can be rewritten as :

or:

(

) then there exists a vector

such that

. Then pre-multiplying (3) by

, it comes:

Thus the system (3) is undetermined and it is not possible de compute

and thus it is not possible to compute

to bring back the state to

( then the system is uncontrollable (

)) except if

is in the controllolability sub-space

characterized by:

where

is the null-space of the matrix 𝒞.

4. Proof of the sufficient condition ():

For the sufficient condition, a new condition

is introduced:

Then

. Indeed, the particular control signal:

satisfies the problem (equation (1)). And thus the system is controllable (

(Indeed: (3) in (1) gives:

To prove

, it can be proven that

is not invertible (

), then considering a vector

, it comes:

Remark : in the single input case:

and

is a scalar quadratic function of τ and thus is always positive or null (

). For its intregral to be null, it requires that:

. Thus:

or using the Cayley-Hamilton expresion of

(equation (2)):

Thus:

5. Discrete-time case:

A linear discrete-time time-invariant system:

is controllable iff: for any initial value

of the state at the initial time

and for a given final time N, there exists a control signal sequence

such that

Let n be the system order (

). A necessary and sufficient condition for the system to be controllable (condition

) is:

is called the controllability matrix .

Proof: the proof in the discrete-time case is simpler than in the continuous-time case. Note that in the discrete-time case, there is a constraint on the final time:

. Indeed, the direct integration of the state equation with

reads:

⋮

Thus, one can solve this last equation in

if and only if

. Then the solution reads:

6. A particular solution in the continuous-time case:

Considering the continuous-time system

with an initial state

, the objective is to find a control signal trajectory

such that

for a given final time

The discrete-time case can provide a solution to this problem considering that

is the output of a zero order hold (ZOH) sampled at

(see the following figure):

The general solution of the continuous-time state equation:

with:

and the change of variable

gives the discrete-time state-space representation of the system between

and

Then the solution (5) can be computed and applied to the input of the zero order hold such that

. Solution (5) provides

, of course:

That is illustrated in the following MATLAB session for a 4-th order continuous-time system with an arbitrary initial condition and a given final time

7. Illustration:

Given:

a continous-time system: ,
an arbitrary initial condition ,
a given final time ,

an open-loop control profile

is computed such that

The proposed example is the serie-connection of:

the plant model which corresponds to an unstable aerospace vehicule: ,
the actuator model: .

and is described by a SIMILINK file:

Numerical application:

J=1;a2=1;f=0.01;xi=0.7;w=10; % Numerical application.

open_system('model_4th_order');

snapshotModel('model_4th_order');

sys=linmod('model_4th_order');
sys.StateName
ans = 4×1 cell array
'model_4th_order/theta'        
'model_4th_order/u_r_dot//w'   
'model_4th_order/u_r'          
'model_4th_order/theta_dot'    

From now, we know that the state vector provided by linmod is:

A=sys.a;B=sys.b;
n=size(A,1);
Gc=ss(A,B,eye(n),0);       % continuous-time system
tf=3;                      % final time
Gd=c2d(Gc,tf/n);           % discrete-time system
Ad=Gd.a;Bd=Gd.b;
Cd=Bd;                     % controllability matrix
for ii=1:n-1;
    Cd=[Ad*Cd(:,1) Cd];    % here : Cd=[Ad^(n-1)Bd,..., AdBd, Bd].
end
rank(Cd)                   % check the system is controllable
ans = 4

We consider an initial condition on the attitude θ:

x0=[1;0;0;0]; % initial condition

u_n=-inv(Cd)*Ad^n*x0 % [u_0, u_1, u_2, u_3]^T

u_n = 4×1

% Simulation:

N=1000;

t=[0:tf/n/N:tf]';

U1=[];

for ii=1:n,

U1=[U1;u_n(ii)*ones(N,1)];

end

U1=[U1;0];

figure

lsim(Gc,U1,t,x0,'zoh');

grid

the stair-shape signal

due to the ZOH sampled at

can be seen on the various plot (grey lines).

8. Another solution in continuous-time:

This solution is directly computed in continuous-time and is given by equation (4). Indeed

is the solution of the Lyapunov equation:

Proof:

Thus

can be easily computed using the MATLAB function lyap.

Then the control signal bringing the state to

(

) is:

That is illustrated in the following MATLAB sequence:

Qc_tf=lyap(A,-B*B'+expm(-A*tf)*B*B'*expm(-A'*tf))
Qc_tf = 4×4
1019 ×
    0.0001    0.0011    0.0075   -0.0012
    0.0011    0.4137   -0.4917    0.0270
    0.0075   -0.4917    1.5170   -0.1627
   -0.0012    0.0270   -0.1627    0.0201

Clearly

is ill-conditionned !! Such a solution is thus very sensitive to the condition number of

and works ony for very small values of

(this is due to fact that the system is unstable):

tf=0.5;

Qc_tf=lyap(A,-B*B'+expm(-A*tf)*B*B'*expm(-A'*tf))

Qc_tf = 4×4

10³ ×

0019    0.1112    0.0024   -0.0123
1112    8.1208   -1.8612   -0.5845
0024   -1.8612    2.3109   -0.1542
   -0.0123   -0.5845   -0.1542    0.0863

Qc_tf_m1=inv(Qc_tf);

t=[0:tf/1000:tf]; % requires a very small step!!

U2=[];

for ii=1:length(t),

U2=[U2;-B'*expm(-A'*t(ii))*Qc_tf_m1*x0];

end

Y=lsim(Gc,U2,t,x0);

figure

subplot(5,1,1)

plot(t,Y(:,1));ylabel('theta');

subplot(5,1,2);

plot(t,Y(:,2));ylabel('u_r_dot/w');

subplot(5,1,3)

plot(t,Y(:,3));ylabel('u_r');

subplot(5,1,4);

plot(t,Y(:,4));ylabel('theta_dot');

subplot(5,1,5);

plot(t,U2);ylabel('u');

One can see that the magnitude of the control signal (last plot) is very big (

9. A closed-loop solution:

In the two previous illustrations, the control signal

is applied in open-loop. Although it meets the constraint

, the system is still unstable for

. From a practical point of view a closed-loop solution is preferred.

Then, one can use a variation of the controllability definition (for linear system): a linear system

is controllable iff one can compute a state feedback assigning the n closed-loop eigenvalues to any prescribed values

. Of course one will choose these n eigenvalues stable and fast enough to bring back the state to

efficiently. But this is an exponential converge to

, contrary to the initial definition.

This new definition allows to introduce the Controllability Canonical Form (CCF) (

) of the linear system:

One can recognize:

the coefficient of the characteristic polynomial of the open-loop system () on the last row of the matrix (with a minus sign and by increasing order of power of s),
the matrix with 0 every where except 1 in the last row.

The objective is thus to find a state feedback control law:

such that the closed-loop dynamics (

) fits the prescribed dynamics:

Then the solution is obvious :

Indeed:

It is shown that an arbitrary linear system

can be tansformed into the CCN by a regular change of state vector

iff:

Proof: It is well known that:

and

. We are going to identify

(column per column, starting from the last column ) form the relations

and

and then check that

is invertible:

The identification of

provides the last column of matrix

The identification of the last column of

provides

The identification of the column

provides

⋮

The identification of the column 2 of

provides

Since all the colmuns

of matrix

are now defined, we have to check that the identification of the column 1 of

works. That is confirmed thanks to the Cayley-Hamilton theorem (see: CayleyHamilton). Indeed, we get:

Finally, we have to check that

is invertible. Indeed, since all its columns

are linear combinations of

, then:

Once

and

are known, one can compute the state feedback gain

to be applied on the initial physical state

. Indeed:

Illustration: for a given n-th order system

and a given vector of n specifed eigenvalues

, the computation of

and

are embedded in the MATLAB function acker (see also: https://en.wikipedia.org/wiki/Ackermann%27s_formula ):

help acker
 acker  Pole placement gain selection using Ackermann's formula.
 
    K = acker(A,B,P)  calculates the feedback gain matrix K such that
    the single input system
            .
            x = Ax + Bu 
 
    with a feedback law of  u = -Kx  has closed loop poles at the 
    values specified in vector P, i.e.,  P = eig(A-B*K).
 
    Note: This algorithm uses Ackermann's formula.  This method
    is NOT numerically reliable and starts to break down rapidly
    for problems of order greater than 10, or for weakly controllable
    systems.  A warning message is printed if the nonzero closed-loop
    poles are greater than 10% from the desired locations specified 
    in P.
 
    See also  place.

Considering the previous example

, one can choose to assign the 4 closed loop eigenvalues to:

: i.e. the low-frequency closed-loop dynamic corresponds to a stable 2-nd order with a damping ratio of and a frequency of ,
the 2 roots of : i.e. the actuator dynamics is unchanged.

P=[-1+sqrt(-1);-1-sqrt(-1);roots([1 14 100])];

K=acker(A,B,P);

damp(A-B*K)

Pole Damping Frequency Time Constant (rad/TimeUnit) (TimeUnit) -7.00e+00 + 7.14e+00i 7.00e-01 1.00e+01 1.43e-01 -7.00e+00 - 7.14e+00i 7.00e-01 1.00e+01 1.43e-01 -1.00e+00 + 1.00e+00i 7.07e-01 1.41e+00 1.00e+00 -1.00e+00 - 1.00e+00i 7.07e-01 1.41e+00 1.00e+00

% Closed-loop system with the 4 states and the control signal on the

% output:

G_CL=ss(A-B*K,B,[eye(4);-K],0);

% Simulation

figure

initial(G_CL,x0)