% Exercise Session 1

clear
close all

load ex1_data

%
% Question 1.A
% Poles

eig(A)

%
% Question 1.B
% Stability?
% There are

abs(eig(A)) > 1

% Yes, the system is unstable
%
% Qusetion 1.C
% Transmission Zeros

tzero(A,B,C,D)

%
% Question 1.D
% Choose a trasmission zero, call it \zetaand solve the corrisiponding homogenous system (corse p.84)
% \left\lbrack \begin{array}{cc}
% \zeta I-A & -B\\
% C & D
% \end{array}\right\rbrack\left\lbrack \begin{array}{c}
% x_0 \\
% u_0
% \end{array}\right\rbrack = 0 \to \mathrm{Mz}=0
% To solve this equation and find the M matrix, there are two approaches. The first one is solving the generalized eigen value problem and the second one is using the null command as following.
% The theory (course p. 82) teaches us that if \zeta is present a trasmission zero, then the matrix M is rank-deficient. Therefore, it is possible to find a non-trivial solution z (z \not= 0) of the aforementioned homogenous system. To become this solution, look for the vectors in the right zero space of M which can be done with the following commands:

zz=tzero(A,B,C,D)
zeta=zz(1); % choose one of the transmission zeros
M=[zeta*eye(length(A))-A -B; C D];
z=null(M);
x0=z(1:length(A),1);
u0=z(length(A)+1:length(z),1);

% When the zero space is higher dimensional, then z contains multiple columns. Each linear combination of the columns of z is a solution for the system.
% To prove that the inital state and input transmission zero drive the output to zero, we need, firstly, define the input as u[k] = u_0 \zeta^k, where k \in \mathbb{N}.

T = 10
k = 0:1:T
sys = ss(A,B,C,D, 1);
u = u0*zeta.^k
lsim(sys, u, k, x0)

% Since the selected transmission zero is not inside the unit circle, we can verify that due to the effect of the input power to k (u_0\zeta^k), there is an explosion after a while in the output. This phenomenon is not visible, if we choose a transmission zero inside the circle as the following example.
%

zeta=zz(6);
M=[zeta*eye(length(A))-A -B; C D];
z=null(M);
x0=z(1:length(A),1);
u0=z(length(A)+1:length(z),1);
u = u0*zeta.^k
lsim(sys, u, k, x0)

% As the plot shows, the output stays close to the zero axes. The randomness of the value is due to a round error.
%
% Question 1.E
% Transfer function from input 1 to output 2

[num,den]=ss2tf(A,B,C,D,1);

num(1,:) % numerator of the unknown transfer function
den(1,:) % denumenator
roots(num(1,:)) % zeros
roots(den(1,:)) % poles = to question 1A

% The relation between the zeros of the SISO-subsystem and the transmission zeros of the MIMO-system is not that clear. Not all transmission zeros have to be zeros of num(1; :) and also reverse, not all SISO-zeros have to be transmission zeros.
%
% Question 2.A
% Controllability Matrix Rank Test

Cont = ctrb(A, B)
length(A) - rank(Cont) % uncontrollable state

% PBH Test

[V, d] = eig(A')
P = B'*V

% Check the columns in P equal to zero. (If any) The eienvalues corresponds to these columns (look to the corrisponding elements in the diaonal matrix d) refer to the uncontrollable modes.
%
% Question 2.B
% Observability Matrix Rank Test

OO=obsv(A,C);
rank(OO) % < length(A) ?
length(A)

% PBH Test

[V,d]=eig(A);
d
C*V

% Also here, loo for the columns that contain only zeros and then look at the corresponding eigenvalues
%
% Question 2.C
% The stability of the uncontrollable/unobservable modes is determined by the eigenvalues that correspond to the zero columns in the PBH test. Alternatively, the Kalman decomposition also allows you to find the controllable and/or observable modes (course pages 72 to 74).
% If there exist unstable uncontrollable modes, the system is not stabilizable. If there are unstable unobservable modes, then the system is not detectable.
%
% Question 2. D
% Minimal Realization

sys=ss(A,B,C,D,-1);
sysm=minreal(sys);
[Am,Bm,Cm,Dm]=ssdata(sysm);

% Now, only the controllable and observable subsystem remains. You can get the same solution up to a state transformation (see course p. 75) from the Kalman decomposition. Calculating the Kalman decomposition from scratch is not needed and all the process can be done using the minreal command with the tolerance, which is linked to the threshold that is used to select the most important component in svd.
%
% Question 3.A
% The control canonical form is given by:
% A_c =\left\lbrack \begin{array}{cc}
% -7 & -12\\
% 1 & 0
% \end{array}\right\rbrack
% B_c =\left\lbrack \begin{array}{c}
% 1\\
% 0
% \end{array}\right\rbrack
% C_c =\left\lbrack \begin{array}{cc}
% 1 & a
% \end{array}\right\rbrack
% The observer canonical form is given by:
% A_o =\left\lbrack \begin{array}{cc}
% -7 & 1\\
% -12 & 0
% \end{array}\right\rbrack
% B_o =\left\lbrack \begin{array}{c}
% 1\\
% a
% \end{array}\right\rbrack
% C_c =\left\lbrack \begin{array}{cc}
% 1 & 0
% \end{array}\right\rbrack
%
% Question 3.B
% The control canonical form is controllable for all values of a. However, it is unob�servable for a = 3 and a = 4 (pole zero cancellation in the transfer function). The observer canonical form is observable for all values of a. However, it is uncontrollable for a = 3 and a = 4.
%
% Question 3.C
% From the previous observations we learn that it is in general not possible to decide observability/controllability from a transfer function, it is a function of the state of the system.
%
% Question 3.D
% The following statements are true:
%  – If a system is controllable, its state space representation can always be converted to control canonical form (via a nonsingular linear transformation).
% – If a system is observable, its state space representation can always be converted to observer canonical form (via a nonsingular linear transformation).
%  – Controllability and observability are invariant under nonsingular linear transformations.
% In the above example, it is not possible to transform the observer canonical form into the control canonical form for a = 3 or a = 4 (or vice versa) via a nonsingular transformation because the systems are not observable/controllable. Hence, these systems are not equivalent for those values of a. On the other hand, if you would investigate any other value of parameter a, you would easily find that there the linear nonsingular transformation between both realizations does exist.