% Script for exercises in EGR 390 class #36, Chapter 13
% J Cardell, April 24, 2009

% 1: Test observability and controllability of following system
A = [0 1 0 0; 0 0 1 0; 0 0 0 1; -1 -4 0 -3];
B = [0; 0; 0; 1];
C = [1 0 0 0];

[M, D] = eig(A);
Bn = inv(M)*B
Cn = C*M

% Controllability Matrix
P = [B A*B A*A*B A*A*A*B]
Prank = rank(P)
Pdet = det(P)

% Observability
% Since (AB)' = B'A', Q is calculated as below, though writing the
% expression for Q in-line text is Q' = [C' A'C' (A^2)'C' ... ]
Q = [ C; C*A; C*A*A; C*A*A*A]
Qrank = rank(Q)
Qdet = det(Q)

% 2: Test observability and controllability of following system
A = [0 1; 1 1.5];
B = [1; 2];
C = [1 0];

[M, D] = eig(A);
Bn = inv(M)*B
Cn = C*M

P = [B A*B]
Prank = rank(P)
Pdet = det(P)

Q = [ C; C*A ]
Qrank = rank(Q)
Qdet = det(Q)

% 3: Find B matrix for #2, B = [1; q];  so find 'q' to make the system
% controllable.  This can be done (a) by constructing P and finding values
% of 'q' that make det(P)=0 - then any OTHER values will make system
% controllable;  or (b) construct Bn = inv(M)*B and find values of 'q' that
% make any row of Bn = 0.  Then, any OTHER value of q will make system
% controllable

A = [0 1; 1 1.5];
% B = [1; q];

[M, D] = eig(A);

% 4:  Find feedback gain matrix K that will place poles at -1 and -10
% To calculate by hand, find det([sI - A + BK]) = 0.  Plug in s1 = p1 = -1,
% and s2 = p2 = -10, write out the polynomial and solve for K1 and K2
A = [0 1; -2 -3];
B = [1; 1];

K = place(A, B, [-1 -10])

% 5:  Pole placement for mass-spring system
% See mass_spring_script.m script

% System matrix A, for C = 0.5 ('c' is the damping constant 'b')
% A = [0 1; -k/m -c/m]
% k = 2
% m = 1
% 
% A =      0    1.0000
%    -2.0000   -0.0500
% 
% Eigenvalues for C = 0.5
% D = -0.2500 + 1.3919i        0          
%          0            -0.2500 - 1.3919i

% For selecting pole placement, note that Re{s} = 1/tau so from the
% knowledge (chapter 7 of EGR 220) that at one time constant the
% amplitude is decreased to 0.368 of max value, we find that for c =
% 0.5, the max is around 0.55, (or > 0.6 if extrapolating to axis), and
% 36.8% of this is ~0.20, which occurs around 4 sec which = 1/s = 1/.25
%
% Period = 4.5s = 1/f = 1/{omega/(2*pi)} = 1/{1.3919/(2*pi)}