clear;
close all;
echo on;
m = 2;
c = 0.1;
k = 1;
who
s = tf('s');
G1 = 1/(m*s*s + c*s + k); % definition of transfer function
A = [0, 1; -k/m, -c/m];
B = [0; 1/m];
C = [1, 0];
D = 0;
G2 = ss(A, B, C, D); % definition of state-space system
bode(G1, G2); % equivalent?
figure(2);
bode(G1, 'b', G2, 'r--'); % for better appearance
G3 = tf(G2); % TF and SS can be converted to each other
G3 % same to G1?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [Exercise 1]
% Confirm the following by drawing Bode diagram of 2nd-order systems
% - peak gain becomes higher as zeta becomes small
% - frequency at peak gain is omega_n
% - phase-lag is up to 180 degree
% - rate of gain attenuation is 40dB/dec at high frequency
zeta = 0.5;
omega_n = 1;
T1 = omega_n^2/(s^2 + 2*zeta*omega_n*s + omega_n^2)
zeta = 0.05;
omega_n = 10;
T2 = omega_n^2/(s^2 + 2*zeta*omega_n*s + omega_n^2)
figure(3);
bode(T1, 'b', T2, 'g');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T3 = 1/(s^3+1.5*s^2+1.5*s+1); % more complex plant
[p, z] = pzmap(T3) % stable?
roots(T3.den{1}) % stable?
get(T3) % show property of T3
T3.num{1} % numerator can be also extracted
figure(4);
nyquist(T3); % closed-loop is stable?
figure(5);
bode(T3); % GM and PM? relation to Nyquist plot?
Gcl = feedback(T3, 1);
help feedback % see on-line help
doc feedback % you can read detail documents in English
pole(Gcl) % Stability check can be also done by closed-loop poles
Gcl2 = T3/(1+T3) % Caution! This introduces redundant states
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [Exercise 2]
K = 1.5 % consider the static controller
% Answer if the closed-loop is stable or not.
L = K*T3; % loop transfer function
Gcl2 = feedback(L, 1);
pole(Gcl2)
figure(4);
hold on
nyquist(L, 'r'); % unstable
figure(5);
hold on;
bode(L, 'r'); % unstable
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% robustness against gain variation
pause
figure(6);
nyquist(T3);
hold on;
for K=0.5:0.1:1.5
K
L = T3*K; % loop transfer function
nyquist(L);
end
% relationship between poles and impulse response
pause
omega_n = 2*pi*2; % 2Hz
zeta = -0.05;
H1 = omega_n^2/(s^2 + 2*zeta*omega_n*s + omega_n^2)
zeta = 0;
H2 = omega_n^2/(s^2 + 2*zeta*omega_n*s + omega_n^2)
zeta = 0.05;
H3 = omega_n^2/(s^2 + 2*zeta*omega_n*s + omega_n^2)
figure(7);
pzmap(H1, H2, H3)
figure(8);
impulse(H1, H2, H3, 5);
figure(9);
step(H1, H2, H3, 5);
% Simulink
mod1 % open a model file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [Exercise 3]
% Confirm the following by executing simulation
% 1. By default setting of SS(initial state is 0), the outputs is consistent with TF.
% 2. By setting the initial state as [1; 0], then the output begin with y(0) = 1.
% 3. By setting the initial state as [5; 0], then the output begin with y(0) = 5.
% 4. By setting the initial state as [0; 1], then the output begin with y(0) = 0 but the rate of change is faster
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%