%
% Armando A. Rodriguez
% Copyright 1998
%

%
% EEE480: Feedback Systems
% FINAL EXAM: Problems 1-4
% Spring 1998
%

clear    % clear all variables from the workspace
clf      % deletes all children of the current figure with visible handles.

%
%*******************************************************************************
%

%
% Problem #1: Laplace Transforms, Second Order Systems In Standard Form, Time Constants, 
%             Overshoot, Time-to-Peak, Steady State
%
%             Plant:         P(s) = 1/(s+2)
%             Compensator:   K(s) = k/s
%
% Given the above plant and compensator, the closed loop characteristic 
% equation is given by:
%
%             Phi_cl = s^2 + 2s + k                    (*)
%
% and the closed loop poles are given by:
%
%             s = -1 +/- j SQRT( k - 1 )
%
% The desired time constant specification of tau = 1 second requires that the real part of the
% closed loop poles be at -1. The previous equation shows that this is satisfied for any k 
% greater than or equal to 1.
%
% The desired time-to-peak specification of t_p = pi/2 = pi/w_d seconds requires that the 
% imaginary part of the closed loop poles be at w_d = 2.
%
% Given the above, it follows that the stated closed loop design specifications require that 
% the closed loop poles be at:
%           
%             s = -1 +/- j 2. 
%
% This implies that the desired closed loop charactersistic equation be:
% 
%             Phi_cl = s^2 + 2s + 5.                 (**)
%
% Comparing (*) and (**), it follows that the required compensator design parameter k must be: 
%             
%             k = 5.
%
% Now lets use MATLAB to verify that the closed loop response to a step reference command 
% does in fact exhbit a time constant of tau = 1 second (i.e. settling time of 5tau = 5 seconds)
% and a time-to-peak of  
%                        tp = pi/2 = 1.5708 seconds.
%
num = 5;
den = [1 2 5];
sys = tf(num,den);
t= [0:0.05:5];
y = step(sys,t);
plot(t,y)
title('Response to Step Reference Command (Sp98 Final, Problem 1)')
xlabel('Time (seconds)')
ylabel('Output y')
grid
pause
%return


%
% Problem 1(b): Use MATLAB to compute coefficients in partial fraction expansion
%
s = 0;
A = 5/(s^2+2*s+5)                  % Coefficient associated with step command r.
                                   % A = T_ry (0) = 1.
                                   % That is, A is equal to the dc gain of the transfer 
                                   % function from r to y.

s = j*100;
B = 0.5*-5/(s^2+2*s+5)             % Coefficient associated with sinusoidal noise n.
                                   % B = 0.5 T_ny(j100).
magB = abs(B)                      % Magnitude of coefficient B.
angleB = angle(B)*180/pi           % Angle of coefficient B.                                   




s = -1+j*2;
C = 5*1e4/(s*(s^2+1e4)*(s+1+j*2))  % Coefficient associated with excitation of system's
                                   % natural dynamics by r and n.
magC = abs(C)                      % Magnitude of coefficient C.
angleC = angle(C)*180/pi           % Angle of coefficient C.                                   

%
% Problem 1(b): Use MATLAB to plot closed loop time response y to  a step reference
% command r = 1 and a sinusoidal sensor noise input n = cos100t.
%
t= [0:0.0005:5];
rmn = ones(1,size(t)*[0 1]') - cos(100*t);                     % From r minus n
plot(t,rmn)
title('Input to T_{ry}: r - n = 1 - cos 100t (Sp98 Final, Problem 1b)')
xlabel('Time (seconds)')
ylabel('Output y')
grid
pause
%return

y = lsim(sys,rmn,t);
plot(t,y)
title('Response to r = 1 and n = cos 100t (Sp98 Final, Problem 1b)')
xlabel('Time (seconds)')
ylabel('Output y')
grid
pause
%return

%
% Note: The effect of the noise n on the output y is quite small.
%       This is because n = cos 100t is a very high frequency signal and is
%       thus attenuated significantly by the transfer function from  n to y.
%       The sensor noise n is attenuated by a factor of 0.0005 (-66.0206 dB).
%
%       The stady state ouput y is now visualized.

t= [0:0.0005:15];
rmn = ones(1,size(t)*[0 1]') - cos(100*t);
y = lsim(sys,rmn,t);
plot( t(1,20000:30000), y(20000:30000, 1) )
title('Steady State Response to r = 1 and n = cos 100t (Sp98 Final, Problem 1b,c)')
xlabel('Time (seconds)')
ylabel('Output y')
grid
pause
%return

%
% Problem 1d: Compute Overshoot to Step Reference Command
%
zeta = 1/sqrt(5)                              % zeta = 0.4472
%overshoot = exp(-zeta*pi/sqrt(1 - zeta^2))    % Overshoot = 0.2079


%
%*******************************************************************************
%

%
% Problem #2: Steady State Analysis of a Feedback System
%
t= [0:0.001:15];
r = 1 - 0.01/(5e-4)*sin(100*t + pi/4);
yr = lsim(tf([5],[1 2 5]), r, t);       % Determine component of y due to r
d = sin(100*t - pi/4);
yd = lsim(tf([1 0],[1 2 5]), d, t);     % Determine component of y due to d
y = yr + yd;
plot (t, y)
title('Repsonse to r = 1 - 0.01/(5e-4) sin(100*t + pi/4) and  d = sin(100*t - pi/4) (Sp98 Final, Problem 2)')
xlabel('Time (seconds)')
ylabel('Output y')
grid
pause
%return


%
%*******************************************************************************
%


%
% Problem #3: Root Locus, Bode, Crossovers, Stability Margins
%

%
% Approximation #1
%
clf
num1 = 2;
den1 = conv([1 0],[1 3 2]);
[gm1,pm1,wg1,wp1] = margin(tf(num1,den1))
zeros = roots(num1)
poles = roots(den1)
k = [0:0.01:1e1];
rlocus(tf(num1,den1),k)
grid
title('Root Locus Using Nominal 3rd Order Approximation: L = 2/s(s+1)(s+2)')
pause
w = logspace(-2,8,300);
bode(tf(num1,den1),'-.', w)
title('Frequency Response Using Nominal 3rd Order Approximation: L = 2/s(s+1)(s+2)')
pause
%return


%
% Approximation #2
%
clf
num2 = num1*[1 sqrt(2)*4000 (4000)^2]/(4000)^2;
den2 = den1;
[gm2,pm2,wg2,wp2] = margin(tf(num2,den2))
format long
zeros = roots(num2)
poles = roots(den2)
format short
k = logspace(-2,12,1000);
rlocus(tf(num2,den2),k)
grid
title('Root Locus Using Nominal 3rd Order Approximation with 2nd Order Numerator (Zoom Out)')
pause
bode(tf(num1,den1), '-.',  tf(num2,den2), '--', w)
title('Frequency Response Using Nominal 3rd Order Approximation with 2nd Order Numerator')
pause
%
% The above yields:
%                     Gain Crossover Frequency = 0.7494 rad/sec       Phase Margin       = 32.6283 degrees
%                     Phase Crossover Frequency  = 1.4150 rad/sec     Upward Gain Margin = 3.0032      
%                        
%return

%
% Approximation #3
%
clf
num3 = num2*(5*1e6)^2;
den3 = conv(den2,[1 5*1e6]);
den3 = conv(den3,[1 5*1e6]);
[gm3,pm3,wg3,wp3] = margin(tf(num3,den3))
format long
zeros = roots(num3)
poles = roots(den3)
format short
k = logspace(-2,14,1000);
rlocus(tf(num3,den3),k)
grid
title('Root Locus: Zoom Out')
pause
bode(tf(num1,den1), '-.', tf(num2,den2),'--', tf(num3,den3),'-', w)
title('Frequency Response')
pause
%return


%
%*******************************************************************************
%

%
% Problem #4: Design, Root Locus, Bode, Crossovers, Stability Margins
%

%
% Plnat Dynamics: P = 1/(s-1) [1000/(s+1000)]^2
%
%
% To satisfy the zero steady state command following and disturbance rejection specification, 
% we choose
% 
%             K_1 (s) = k/s.
%
% To satisfy the stability and transient specifications, we choose
%
%             K_2 (s) = (s+a)/a.
%

%
% Given the above, the nominal closed loop charactersitic polynomial - obtained by neglecting 
% high frequency plant dynamics - is given by:
%
%                               Phi_cl = s(s-1)+ (k/a) (s+a)
%                                      = s^2 + (k/a  - 1)s + k
%
% The transient specifications call for a time constant 
%
%                               tau = 1 second 
% and a time-to-peak 
%                               t_p = pi/2 seconds.
%
% These require that the nominal closed loop poles be at
%
%                               s = -1 =/- j 2
%
% This, in turn, requires a nominal closed loop characteristic polynomial of
%
%                               Phi_cl = s^2 + 2s + 5
% This, however, implies that
%
%                               k = 5
% and
%                               a = 5/3.

%
% Given the above, the final open loop transfer function is given by:
%
%                    L(s) = 3 [(s+5/3)/(s(s-1))] [1000/(s+1000)]^2
%

%
% Approximation #1: L_1(s) = 3 [(s+5/3)/(s(s-1))] - a low frequency approximation to L(s)
%
clf
num1 = 3*[1 5/3];
den1 = [1 -1 0];
format long
zeros = roots(num1)
poles = roots(den1)
format short
[gm1,pm1,wp1,wg1] = margin(tf(num1,den1))
%
% The above margin command yields the following:
%
% Phase Crossover Frequency = 1.2910 rad/sec    Downward Gain Margin = 0.3333  
% Gain  Crossover Frequency = 3.2254 rad/sec    Phase Margin = 45.4475 degrees
%

%
% To estimate the downward gain margin and the low frequency phase 
% crossover frequency, one uses the low frequency approximation fot L:
%
%   L_1 (s) = 3 [(s+5/3)/(s(s-1))].
%
% This yields the nominal closed loop characteristic equation:
%
%                   Phi_cl = s(s-1) + 3k (s + 5/3)
%                          = s^2 + (3k - 1)s + 5k
%
% From this one obtains that
%                              k_1 = Downward Gain Margin = 1/3
%                              Phase crossover Frequency = sqrt(5/3) 
%                                                        = 1.2910 rad/sec
%
% (k_1, wp1) approximately specifiy the low frequency imaginary crossover 
% on the root locus for L(s).
%
k = logspace(-2,4,1000);
rlocus(tf(num1,den1),k)        % Plot approximate root locus for L 
                               % using low frequency approximation
grid
title('Root Locus: Using Second Order Low Frequency Approximation (Sp98 Final, Problem 4)')
pause
%return
w = logspace(-2,5,300);
bode(tf(num1,den1), '-.', w)
title('Frequency Response: Second Order Low Frequency Approximation (Sp98 Final, Problem 4)')
pause



%
% Approximation #2: L_2(s) = (s) = 3 [(s+5/3)/(s(s-1))] [1000/(s+1000)]^2
%
num2 = 1000^2*num1;
den2 = conv(den1,[1 2000 1e6]);
[gm2,pm2,wp2,wg2] = margin(tf(num2,den2))
%
% The above margin command yields the following:
%
% Phase Crossover Frequency = 1.2945 rad/sec    Downward Gain Margin = 0.3344  
% Gain  Crossover Frequency = 3.2254 rad/sec    Phase Margin = 45.0776 degrees
%
format long
zeros = roots(num2)
poles = roots(den2)
format short
k = logspace(-2,4,1000);
rlocus(tf(num2,den2),k)
grid
title('Root Locus: Zoomed Out (Sp98 Final, Problem 4)')
pause
%return

w = logspace(-2,5,300);
bode(tf(num1,den1), '-.', tf(num2,den2), '-', w)
title('Frequency Response (Sp98 Final, Problem 4)')
pause

return

%
%*******************************************************************************
%