% Problem 8.11
disp('Problem 8.11')
clear
K=4;
n=[1];
d=[1 2 4];

% Transfer function of the system
disp('(a) Uncompensated system:')
printsys(K*n,d,'s')

% Open-loop poles
disp('open-loop poles:')
p=roots(d)

subplot(221)
rlocus(n,d)
axis([-2 0.1 -3 3])
grid
title('(a) Uncompensated system')

% Operating point sd
OS=0.1; ts=3;
xi=sqrt(log(OS)^2/(log(OS)^2+pi^2)); xi=xi+0.005
om_n=-log(0.05*sqrt(1-xi^2))/(xi*ts)

sd=-xi*om_n + om_n*sqrt(1-xi^2)*i

% Controller phase contribution
th_p1=angle(sd-p(1))
th_p2=angle(sd-p(2))
alpha_c=(pi+th_p1+th_p2)/2

% Zero of PD compensator
if alpha_c > pi/2
   %zc=om_n*(xi*tan(pi-alpha_c)+sqrt(1-xi^2))/tan(pi-alpha_c);
else
   zc=om_n*(xi*tan(alpha_c)+sqrt(1-xi^2))/tan(alpha_c)
end
nc=conv([1 zc],[1 zc])

disp('Closed-loop uncompensated system:')
[ncl,dcl]=cloop(K*n,d,-1)
printsys(ncl,dcl,'s')

disp('Closed-loop poles of uncompensated system:')
roots(dcl)

disp('System with a phase-lead compensator:')
n_c=conv(nc,n);
d_c=d;

% Static gain of compensated system
dp1=abs(sd-p(1))
dp2=abs(sd-p(2))
dzc=abs(sd+zc)
disp('Static gain:')
Kc=dp1*dp2/(dzc*dzc)
subplot(222)
rlocus(n_c,d_c)
axis([-20 0.1 -10 10])
grid
title('(b) Compensated system')

disp('Closed-loop compensated system:')
[n_ccl,d_ccl]=cloop(Kc*n_c,d_c,-1);
printsys(n_ccl,d_ccl,'s')

disp('Poles of closed-loop compensated system:')
roots(d_ccl)

subplot(212)
t=0:0.1:5;
[y,x]=step(ncl,dcl,t);
hold
[yc,xc]=step(n_ccl,d_ccl,t);
plot(t,y,'--',t,yc,'-')
grid
xlabel('time [s]')
ylabel('y,yc')

% Find the settling time from the graph
yss=yc(length(yc));
delta=0.05;
j=length(yc);
while abs((yc(j)-yss)/yss) < delta,
   j=j-1;
end
disp('5%-settling time [s] =')
j*0.1

% Find the overshoot from the graph
j=1;
ymax=yc(j);
j=j+1;
while ymax < yc(j),
      ymax=yc(j);
      j=j+1;
end
disp('Maximum percentage overshoot=')
(ymax-yss)*100/yss