Solution 6.8.3.8 The MATLAB dialogue EDU>g = zpk([-2],[0 -4 -30],1) Zero/pole/gain: (s+2) -------------- s (s+4) (s+30) EDU>rlocus(g) EDU>print -deps rl6838a.eps draws and saves the root locus is shown in Figure 1. The root locus gives us an idea of where the line of constant damping ratio will intersect the root locus. The MATLAB program -30 -25 -20 -15 -10 -5 0 5 -8 -6 -4 -2 0 2 4 6 8 Real Axis Imag Axis Figure 1: Root locus z=2 p1 =0 p2 = 4 p3 = 30 zeta = 1/sqrt(2) 1 22 22.1 22.2 22.3 22.4 22.5 22.6 22.7 22.8 22.9 23 -181 -180.5 -180 -179.5 -179 -178.5 -178 -177.5 Figure 2: Angle versus real part of s omegan =linspace(22,23,200);; s=-zeta*omegan + j* omegan *sqrt(1-zeta^2);; ang = (angle(s + z)-angle(s + p1)-angle(s+p2)-angle(s+p3))*180/pi;; plot(omegan,ang) print -deps 6838a.eps creates the plot shown in Figure 2. Wesee that the root locus crosses the rayofconstantdamping ratio  =1=sqrt2ats = ;16:125+j16:125Indeed, at s = ;16:125+ j16:125, where the net angle is 180:00  . Then the gain to place closed loop poles at s = ;16:125+ j16:125 is K = jsjjs +4js +30jj s=;16:125+j16:125 = 456:55: The natural frequency of the closed loop poles is ! n = p 16:125 2 +16:125 2 =22:804: The MATLAB program z=2 p1 =0 p2 = 4 p3 = 30 zeta = 1/sqrt(2) 2 omegan =linspace(22,23,200);; s=-zeta*omegan + j* omegan *sqrt(1-zeta^2);; ang = (angle(s + z)-angle(s + p1)-angle(s+p2)-angle(s+p3))*180/pi;; plot(omegan,ang) print -deps 6838a.eps s=-zeta*omegan(161) + j*omegan(161)*sqrt(1-zeta^2) K=(abs(s + p1)*abs(s+p2)*abs(s+p3))/abs(s+z) omegana = omegan(161) sc = conj(s) tn2 = zpk([],[s sc],omegana^2) g=zpk([-z],[-p1 -p2 -p3],K) tc = feedback(g,1) T=linspace(0,5,200);; [Yn2,T] = step(tn2,T);; [Ytc,T] = step(tc,T);; plot(T,Yn2,'k-',T,Ytc,'k--') print -deps sr6838.eps t =0:0.01:1;; u=t;; [Yrtc,t] = lsim(tc,u,t);; [Yrn2,t] = lsim(tn2,u,t);; plot(t,Yrn2,'k-',t,Yrtc,'k--') print -deps rr6838.eps Prints and plots the step and ramp responses of T c (t)andT N2 (t) shown in Figure 3 and 4, respectively. The partial fraction expansion for the step response is C(s)= 1 s ; 0:136 s +1:756 + 0:6162 6 2:3477 s +16:125; j16:125 + 0:6162 6 ;2:3477 s +16:125+ j16:125 Then c(t)=[1; 0:136e ;1:756t +1:2323e ;16:125t cos(16:125t+2:3477)]1(t): The partial fraction expansion for the ramp response is C(s)= ;0:1315 s + 1 s 2 + 0:07745 s +1:756 + 0:0270 6 ;0:0085 s +16:125; j16:125 + 0:0270 6 0:0085 s +16:125+ j16:125 Then c(t)=[t;0:1315+0:077456e ;1:756t +0:0540e ;16:125t cos(16:125t;0:0085)]1(t): The residues were found using the MATLAB program 3 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 T N2 T c Figure 3: Comparison of step responses 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 T N2 T c Figure 4: Comparison of ramp responses 4 z=2 p1 =0 p2 = 4 p3 = 30 zeta = 1/sqrt(2) omegan =linspace(22,23,200);; s=-zeta*omegan + j* omegan *sqrt(1-zeta^2);; ang = (angle(s + z)-angle(s + p1)-angle(s+p2)-angle(s+p3))*180/pi;; plot(omegan,ang) print -deps 6838a.eps s=-zeta*omegan(161) + j*omegan(161)*sqrt(1-zeta^2) K=(abs(s + p1)*abs(s+p2)*abs(s+p3))/abs(s+z) omegana = omegan(161) sc = conj(s) tn2 = zpk([],[s sc],omegana^2) g=zpk([-z],[-p1 -p2 -p3],K) tc = feedback(g,1) T=linspace(0,5,200);; [Yn2,T] = step(tn2,T);; [Ytc,T] = step(tc,T);; plot(T,Yn2,'k-',T,Ytc,'k--') print -deps sr6838.eps t =0:0.01:1;; u=t;; [Yrtc,t] = lsim(tc,u,t);; [Yrn2,t] = lsim(tn2,u,t);; plot(t,Yrn2,'k-',t,Yrtc,'k--') print -deps rr6838.eps v0 = K*[1 z] v1 = [1 0] v2 = [1 -sc] v3 = [1 conj(-sc)] v4 = [1 1.756] B= K*[1 z] A=conv(v1,v2) A=conv(A,v3) A=conv(A,v4) [R,P,K1] = residue(B,A) M=R(2) absm =abs(M) abs2m = 2*abs(M) angm =angle(M) t=0 dt = 0.025 kount = 1 5 while t < 5 c(kount) = R(4)+ R(3)*exp(-1.756*t) + abs2m*exp(-16.125*t) *cos(16.1235*t + angm);; time(kount) = t;; t=t+dt;; kount = kount + 1;; end plot(time,c) B= K*[1 z] A=conv(v1,v1) A=conv(A,v2) A=conv(A,v3) A=conv(A,v4) [R,P,K1] = residue(B,A) M=R(2) absm =abs(M) abs2m = 2*abs(M) angm =angle(M) t=0 dt = 0.025 kount = 1 while t < 5 c(kount) = R(4)+ R(5)*t + R(3)*exp(-1.756*t) + abs2m*exp(-16.125*t) *cos(16.1235*t + angm);; time(kount) = t;; t=t+dt;; kount = kount + 1;; end plot(time,c) 6