% Problem 8.17
disp('Problem 8.17 - Double PD controller')
clear
n=[0 0 0 0 0 1];
d1=conv([1 0.2],[1 0.5]);
d2=conv([1 10],conv([1 14.28],[1 25]));
d=conv(d1,d2);
disp('Original system (K=1)')
printsys(n,d,'s')

% (a) Root locus
figure(1)
subplot(221)
rlocus(n,d)
axis([-30 5 -10 10])
title('(a) Uncompensated system')
grid
subplot(223)
rlocus(n,d)
axis([-1 0.1 -1 1])
grid

% Double PD controller
nc=conv([1 1.2079],[1 1.2079]);
dc=[0 0 1];

num_c=conv(n,nc);
den_c=conv(d,dc);
subplot(222)
rlocus(num_c,den_c)
axis([-30 5 -10 10])

title('(b) Compensated system')
grid
subplot(224)
rlocus(num_c,den_c)
axis([-2 0.1 -1 1])
grid

% Operating point
xi=0.7071;
disp('Pick up a point from the root locus by using the left mouse button!')
[sdx,sdy]=ginput(1);
om_n=-sdx*sqrt(2)
disp(' ')
disp('Operating point:')
sd=sdx+sdy*i

% 4.17(b)
% Static gain
dp1=abs(sd+0.2);
dp2=abs(sd+0.5);
dp3=abs(sd+10);
dp4=abs(sd+14.28);
dp5=abs(sd+25);
dzc=abs(sd+1.2079);
disp('Static gain:')
K=dp1*dp2*dp3*dp4*dp5/dzc^2

% Steady state errors
% Uncompensated system
Kp=K*n(length(n))/d(length(d))
disp('Steady state error of the uncompensated system:')
ess=1/(1+Kp)

% Compensated system
Kpc=K*num_c(length(num_c))/den_c(length(den_c))
disp('Steady state error of the compensated system:')
essc=1/(1+Kpc)

% Closed-loop poles
[n_cl,d_cl]=cloop(K*n,d,-1);
[nc_cl,dc_cl]=cloop(K*num_c,den_c,-1);
disp('Closed-loop eigenvalues - compensated system:')
pcomp=roots(dc_cl)

% Step responses
t=0:0.05:25;
[y,x]=step(n_cl,d_cl,t);
[yc,xc]=step(nc_cl,dc_cl,t);

% Settling time ts
% Uncompensated system
j=length(t); ysteady=y(length(t)); delta=abs(ysteady-y(j))/ysteady;
while delta<0.05,
  j=j-1;
  delta=abs(ysteady-y(j))/ysteady;
end
disp('ts [s] - uncompensated system')
t(j+1)

% Compensated system
j=length(t); ycsteady=yc(length(t)); delta=abs(ycsteady-yc(j))/ycsteady;
while delta<0.05,
  j=j-1;
  delta=abs(ycsteady-yc(j))/ycsteady;
end
disp('ts [s] - compensated system')
t(j+1)

% Maximum percentage overshoot
% Uncompensated system
j=1; k=j+1;
while (y(k)>y(j) | y(k)<ysteady),
  j=k; k=k+1;
end
disp('Overshoot [%] of the uncompensated system:')
OS=(y(j)-ysteady)*100/ysteady

% Compensated system
j=1; k=j+1;
while (yc(k)>yc(j) | yc(k)<ycsteady),
  j=k; k=k+1;
end
disp('Overshoot [%] of the compensated system:')
OS=(yc(j)-ycsteady)*100/ycsteady

figure(2)
plot(t,y,'--',t,yc,'-')
xlabel('time [s]')
ylabel('y, yc')
grid
