%% System Poles and Zeros
% Investigation of how poles and zeros affect the time response of a system
% Adapted from: http://seit.unsw.adfa.edu.au/staff/sites/valu/TEACHING/CT1/effects.pdf
%
% BJ Furman
%
% 27SEP2014
%% System Transfer Function example
% 
syms s t
Ns = 5*s^2 + 10*s;  % Numerator polynomial in s
Ds = s^3 + 5*s^2+11*s+15;  % Denominator polynomial in s
Gs = Ns/Ds;  % The transfer function
den_resp=ilaplace(1/Ds); % The inverse Laplace transform of the denominator gives the time response for the natural motion (homogeneous solution)
imp_resp = matlabFunction(den_resp) % converts the symbolic time function into a numeric function
fplot(imp_resp,[0 10]) % plots the impulse response of the natural motion from t=0 to t=10
grid
hold    % hold the plot for plotting in
tot_resp_t = ilaplace(Gs) % find the total solution (natural + forced)
tot_resp = matlabFunction(tot_resp_t) % converts the symbolic time function into a numeric function
fplot(tot_resp,[0 10])  % plots the total response to an impulse
xlabel('Time, s'); ylabel('Amplitude'); title('Step Response')
hold    % release the plot
%% Using tf() to Generate the Transfer Function
%
sys = tf(1,[1 5 11 15]) % builds the transfer function using the coefficients of s in numerator, denominator from highest order of s to lowest moving from left to right
step(sys)   % plot the step response
xlabel('Time, s'); ylabel('Amplitude'); title('Step Response Using tf()')
grid
%% First Order System With Addition of a Pole on the Real Axis
% System: Y(s)/R(s) = p/[(s+1)*(s+p)] = 1/[(s+p)*((1/p)*s + 1)]
%
% Additional pole added at s = -p tends to slow the system down
pvals=[.5 1 10 inf];    % values of p to generate additional pole on negative real axis
t=[0:.05:5]; % vector of times for plotting
ypvals=zeros(length(t),length(pvals));   % establish a variable to hold the amplitude values of the step response for the cases of different poles
for i=4:-1:1    % iterate through the additional pole values starting with the last value
    p=pvals(i);     % grab the ith additional pole value and hold onto it in the variable p
    G=tf(1,conv([1 1],[1/p 1]));    % form the transfer function by multiplying the factored terms
    [y,x]=step(G,t);    % generate the step response
    ypvals(:,i) = y;    % store the result. Each column of ypvals corresponds to a step response for pvals(i)
end
plot(t,ypvals(:,1),'-',t,ypvals(:,2),'--',t,ypvals(:,3),'-.',t,ypvals(:,4),'.')
xlabel('Time, s'); ylabel('Amplitude'); title('Response of a 1st Order System With the Addition of a Pole on the Real Axis in the Left-Half Plane','FontSize',9)
legend('p=\infty','p=10','p=1','p=0.5','Location','Best');
grid
%% Effect of a Zero on a Dominantly First-Order System
% System: Y(s)/R(s) = [(1/z)*s + 1]/[(s+1)*((1/10)*s+1)]
%
% Note if z=inf, then the system looks like what we had previously with an
% extra pole added at s=-10 to the first order system dominated by the pole
% at s=-1. The additional pole at s=-10 resulted in a slight slow down, but not much, because the
% addtional pole was relatively far to the left from s=-1 down the negative real
% axis
%
% Consider five cases:
%
%  z >> 10 (which means a zero far to the left of the pole at s=-10)
%  z = 10  (which means the zero sits on top of the pole at s=-10, effectively 'canceling' the extra pole
%  1 < z < 10  (which means the zero is between the pole at s=-10 and the pole at s=-1
%  z = 1  (which means the zero sits on top of the pole at s=-1, effectively canceling the pole at s=-1)
%  z < 1  (which means the zero is between the imaginary axis and the dominant pole at s=-1)
%
% The addition of a zero on the negative real axis tends to speed up the
% response, and will lead to overshoot if it becomes dominant (closest to the origin)

zvals=[.2 .5 1 2 10 inf];    % values of z to generate added zero on negative real axis
t=[0:.05:5]; % vector of times for plotting
ynzvals=zeros(length(t),length(zvals));   % establish and pre-allocate a variable to hold the amplitude values of the step response for the cases of different zeros
for i=6:-1:1    % iterate through the additional zero values starting with the last value
    z=zvals(i);     % grab the ith additional z value and hold onto it in the variable z
    G=tf([1/z 1],conv([1 1],[1/10 1]));    % form the transfer function by multiplying the factored terms using conv()
    [y,x]=step(G,t);    % generate the step response
    ynzvals(:,i)=y;    % store the result. Each column of ynzvals corresponds to a step response for zvals(i)
end
plot(t,ynzvals(:,1),'--',t,ynzvals(:,2),'-',t,ynzvals(:,3),'-.',t,ynzvals(:,4),'-o',t,ynzvals(:,5),'-x',t,ynzvals(:,6),'-*')
xlabel('Time, s'); ylabel('Amplitude'); title('Response of a 1st Order System With the Addition of a Zero on the Real Axis in the Left-Half Plane','FontSize',9)
legend('z=0.2','z=0.5','z=1','z=2','z=10','z=\infty','Location','Best');
grid
%% Effect of a Right Half-Plane Zero
% System: Y(s)/R(s) = [(1/z)*s + 1]/[(s+1)*((1/10)*s+1)]
%
% If z < 0, then the added zero will be on the real axis in the right half-plane
%
% Look at dy(t)/dt at t=0  (=10/z). If z<0, then y(t) will start off in the
% negative direction. Like the previous case, the closer the zero
% approaches the origin, the more overshoot will occur (but in the
% *negative*
% direction).

zpvals=[-0.2 -0.5 -2 -inf];    % values of z to generate additional zero on positive real axis 
t=[0:.05:5]; % vector of times for plotting
ypzvals=zeros(length(t),length(zpvals));   % establish and pre-allocate a variable to hold the amplitude values of the step response for the cases of different zeros
for i=4:-1:1    % iterate through the additional pole values starting with the last value
    z=zpvals(i);     % grab the ith additional zero value and hold onto it in the variable z
    G=tf([1/z 1],conv([1 1],[1/10 1]));    % form the transfer function by multiplying the factored terms
    [y,x]=step(G,t);    % generate the step response
    ypzvals(:,i)=y;    % store the result. Each column of ypzvals corresponds to a step response for pvals(i)
end
plot(t,ypzvals(:,1),'-',t,ypzvals(:,2),'--',t,ypzvals(:,3),'-.',t,ypzvals(:,4),'-x')
xlabel('Time, s'); ylabel('Amplitude'); title('Response of a 1st Order System With the Addition of a Pole on the Real Axis in the Right-Half Plane','FontSize',9)
legend('z=-0.2','z=-0.5','z=-2','z=\infty','Location','Best');
grid
%% Step Response for a 2nd Order System as a Function of Damping Ratio
% System: Y(s)/R(s) = wn^2/[s^2 + 2*zeta*wn + wn^2]
%
% Consider the step response for several values of zeta, the damping ratio

wn = 1;
zeta_crit = sqrt(2)/2;
zeta = [ 0 0.1 0.3 0.5 zeta_crit 1];
t=[0:.05:25]; % vector of times for plotting
y2ndvals =zeros(length(t),length(zeta));   % establish and pre-allocate a variable to hold the amplitude values of the step response for the cases of different zeta values
for i=1:length(zeta)
    G(i) = tf(wn^2,[1 2*zeta(i)*wn wn^2]);
    [y,x] = step(G(i),t);
    y2ndvals(:,i) = y;
end
subplot(2,1,1)
plot(x,y2ndvals)
xlabel('Time, s'); ylabel('Amplitude'); title('Response of a 2nd Order System For Several Damping Ratios','FontSize',9)
legend('\zeta=0','\zeta=0.1','\zeta=0.3','\zeta=0.5',['\zeta=',num2str(zeta_crit,3)],'\zeta=1');
grid
subplot(2,1,2)
pzmap(G)
axis([-1.05 0.05 -1.1 1.1])
sgrid
%% Step Response for a 2nd Order System as Function of Natural Frequency (Constant Damping Ratio)
% System: Y(s)/R(s) = wn^2/[s^2 + 2*zeta*wn + wn^2]
%
% Consider the step response for several values of wn, but holding zeta, the damping ratio constant at 0.707

zeta = 0.707;
wn = [ 0.5 1 2];
t=[0:.05:25]; % vector of times for plotting
y2ndvals =zeros(length(t),length(wn));   % establish and pre-allocate a variable to hold the amplitude values of the step response for the cases of different zeta values
for i=1:length(wn)
    G(i) = tf(wn(i)^2,[1 2*zeta*wn(i) wn(i)^2]);
    [y,x] = step(G(i),t);
    y2ndvals(:,i) = y;
end
subplot(2,1,1)
plot(x,y2ndvals)
xlabel('Time, s'); ylabel('Amplitude'); title('Response of a 2nd Order System For Several Natural Frequencies at \zeta=0.707','FontSize',9)
legend(['\omega','n=0.5'],['\omega','n=1'],['\omega','n=2']);
grid
subplot(2,1,2)
pzmap(G)    % pole-zero map function
axis([-2 0.05 -2 2])
sgrid
%% Step Response Information Using stepinfo()
% We can get lots of information for a given transfer function using the
% stepinfo() command
close all
zeta = 0.5;
wn = 1;
G = tf(wn^2,[1 2*zeta*wn wn^2]);    % form the transfer function
S = stepinfo(G)
step(G)
xlabel('Time, s'); ylabel('Amplitude'); title(['Response of a 2nd Order System for ','\omega','n = 1, ','\zeta = 0.5'])
grid
%% Effects of Additional Left Half-Plane Zero on the Performance of a 2nd Order System
% 
zvals=[.5 1 2 10 inf];    % values of z to generate added zero on negative real axis
t=[0:.05:5]; % vector of times for plotting
wn = 1;
zeta = sqrt(2)/2;
ynzvals=zeros(length(t),length(zvals));   % establish and pre-allocate a variable to hold the amplitude values of the step response for the cases of different zeros
for i=length(zvals):-1:1    % iterate through the additional zero values starting with the last value
    z=zvals(i);     % grab the ith additional z value and hold onto it in the variable z
    G=tf([wn^2/z, wn^2],[1, 2*zeta*wn, wn^2]);    % form the transfer function by multiplying the factored terms using conv()
    [y,x]=step(G,t);    % generate the step response
    ynzvals(:,i)=y;    % store the result. Each column of ynzvals corresponds to a step response for zvals(i)
end
plot(t,ynzvals(:,1),'-*',t,ynzvals(:,2),'-o',t,ynzvals(:,3),'-.',t,ynzvals(:,4),'--',t,ynzvals(:,5),'-')
xlabel('Time, s'); ylabel('Amplitude'); title('Response of a 1st Order System With the Addition of a Zero on the Real Axis in the Left-Half Plane','FontSize',8)
legend('z=0.5','z=1','z=2','z=10','z=\infty','Location','Best');
grid
%% Effects of Additional Right Half-Plane Zero on the Performance of a 2nd Order System
% 
zvals = [-0.5 -1 -2 -10 -inf];    % values of z to generate added zero on negative real axis
t=[0:.05:5]; % vector of times for plotting
wn = 1;
zeta = sqrt(2)/2;
ynzvals=zeros(length(t),length(zvals));   % establish and pre-allocate a variable to hold the amplitude values of the step response for the cases of different zeros
for i=length(zvals):-1:1    % iterate through the additional zero values starting with the last value
    z=zvals(i);     % grab the ith additional z value and hold onto it in the variable z
    G=tf([wn^2/z, wn^2],[1, 2*zeta*wn, wn^2]);    % form the transfer function by multiplying the factored terms using conv()
    [y,x]=step(G,t);    % generate the step response
    ynzvals(:,i)=y;    % store the result. Each column of ynzvals corresponds to a step response for zvals(i)
end
plot(t,ynzvals(:,1),'-*',t,ynzvals(:,2),'-o',t,ynzvals(:,3),'-.',t,ynzvals(:,4),'--',t,ynzvals(:,5),'-')
xlabel('Time, s'); ylabel('Amplitude'); title('Response of a 1st Order System With the Addition of a Zero on the Real Axis in the Right-Half Plane','FontSize',8)
legend('z=0.5','z=1','z=2','z=10','z=\infty','Location','Best');
grid
%% The Effect of an Additional Forward-path Pole
% Open Loop System: Y(s)/R(s) = wn^2/[s*(s + 2*zeta*wn)(1 + s/p)]
% 
% Closed Loop System: Y(s)/R(s) = wn^2/[(1/p)*s^3 + (1+2*zeta*wn/p)*s^2 + 2*zeta*wn*s + wn^2]
%
% When p is large, the added pole has little effect on the system step
% response, because the system stays dominantly 2nd order. As p moves
% closer to the origin, the overshoot increases. If the pole gets too close
% to the origin, the system can go unstable.
%% The Effect of an Additional Closed Loop Pole
% System Y(s)/R(s) = wn^2/[(s/p + 1)*(s^2 + 2*zeta*wn + wn^2)]
pv=[0.1 .5 2 10 inf];
t=[0:.3:15];
wn=1;
zeta=sqrt(2)/2;
yv=[];
for i=5:-1:1
p=pv(i);
G=tf(wn^2,conv([1/p 1],[1 2*zeta*wn wn^2]));
[y,x]=step(G,t);
yv=[yv,y];
end
plot(t,yv(:,1),'-',t,yv(:,2),'--',t,yv(:,3),'-.',t,yv(:,4),'-o',t,yv(:,5),'-*');
grid;
legend('p=\infty','p=10','p=2','p=0.5','p=0.1');
%%