H T I R T E E N
Digital Control Systems
SOLUTIONS TO CASE STUDIES CHALLENGES
Antenna Control: Transient Design via Gain
a. From the answer to the antenna control challenge in Chapter 5, the equivalent forward transfer function found by neglecting the dynamics of the power amplifier, replacing the pots with unity gain, and including the integration in the sample-and-hold is
Ge(s) =
But,
0.16K s (s + 1.32)
2
z-1 Thus, Ge(z) = 0.16K z Gz , or, (z+0.95696) Ge(z) = 7.659x10-4K (z-1) (z-0.87634) b. Draw the root locus and overlay it over the = 0.5 (i.e. 16.3% overshoot) curve.
.
Solutions to Case Studies Challenges 581
We find that the root locus crosses at approximately 0.93 j0.11 with 7.659x10-4K = 8.63x10-3. Hence, K = 11.268. c.
1 (7.659x10-4 K)(1.95696) lim(z - 1)Ge (z) = = 0.1366; T z 1 0.12366 1 e() = = 7.321 Kv Kv =
d. Program:
T=0.1; numf=0.16; denf=[1 1.32 0 0]; 'F(s)' F=tf(numf,denf) numc=conv([1 0],numf); denc=denf; %Input sampling time %Numerator of F(s) %Denominator of F(s) %Display label %Display F(s) %Differentiate F(s) to compensate %for c2dm which assumes series zoh %Denominator of continuous system %same as denominator of F(s) %Form continuous system, C(s) %Cancel common poles and zeros %Convert to z assuming zoh
C=tf(numc,denc); C=minreal(C,1e-10); D=c2d(C,T,'zoh'); 'F(z)' D=minreal(D,1e-10) %Cancel common poles and zeros and display rlocus(D) pos=(16.3); z=-log(pos/100)/sqrt(pi^2+[log(pos/100)]^2); zgrid(z,0) title(['Root Locus with ' , num2str(pos), ' Percent Overshoot Line']) [K,p]=rlocfind(D) %Allows input by selecting point on %graphic
Computer response:
ans = F(s) Transfer function: 0.16 -------------s^3 + 1.32 s^2 ans = F(z) Transfer function: 0.0007659 z + 0.000733 ---------------------z^2 - 1.876 z + 0.8763 Sampling time: 0.1 Select a point in the graphics window selected_point = 9.2969e-001 +1.0219e-001i K = 9.8808e+000
582 Chapter 13: Digital Control Systems
p = 9.3439e-001 +1.0250e-001i 9.3439e-001 -1.0250e-001i
Antenna Control: Digital Cascade Compensator Design
a. Let the compensator be KGc(s) and the plant be Gp (s ) =
0.16 . For 10% overshoot and a s( s + 1.32)
peak time of 1 second, = 0.591 and n = 3.895, which places the dominant poles at 2.303 j3.142. If we place the compensator zero at 1.32 to cancel the plant's pole, then the following geometry results.
j3.142 s-plane
-pc
X
-2.303
X
Hence, pc = 4.606. Thus, Gc (s )
=
K (s + 1.32) 0.16 K and Gc (s )G p ( s ) = . Using the ( s + 4.606) s(s + 4.606)
product of pole lengths to find the gain, 0.16K = (3.896)2, or K = 94.87. Hence,
Solutions to Case Studies Challenges 583
Gc (s ) =
is Gc (z )
94.87(s + 1.32) . Using a sampling interval of 0.01 s, the Tustin transformation of Gc(s) (s + 4.606)
= 93.35( z - 0.9869) 93.35z - 92.12 = . ( z - 0.955) z - 0.955
(z - 0.955)X(z) = (93.35z 92.12)E(z)
b. Cross multiplying, Solving for the highest power of z operating on X(z), zX(z) = (93.35z 92.12)E(z) + 0.955X(z) Solving for X(z), X(z) = (93.35 92.12z-1)E(z) + 0.955z-1X(z) Implementing this equation as a flowchart yields the following diagram
c.
Program:
's-plane lead design for Challenge - Lead Comp' clf %Clear graph on screen. 'Uncompensated System' %Display label. numg=0.16; %Generate numerator of G(s). deng=poly([0 -1.32]); %Generate denominator of G(s). 'G(s)' %Display label. G=tf(numg,deng); %Create G(s). Gzpk=zpk(G) %Display G(s). pos=input('Type desired percent overshoot '); %Input desired percent overshoot. z=-log(pos/100)/sqrt(pi^2+[log(pos/100)]^2); %Calculate damping ratio. Tp=input('Type Desired Peak Time '); %Input desired peak time. wn=pi/(Tp*sqrt(1-z^2)); %Evaluate desired natural frequency. b=input('Type Lead Compensator Zero, (s+b). b= '); %Input lead compensator zero. done=1; %Set loop flag. while done==1 %Start loop for trying lead %compensator pole. a=input('Enter a Test Lead Compensator Pole, (s+a). a = '); %Enter test lead compensator pole. numge=conv(numg,[1 b]); %Generate numerator of Gc(s)G(s). denge=conv([1 a],deng); %Generate denominator of Gc(s)G(s). Ge=tf(numge,denge); %Create Ge(s)=Gc(s)G(s). clf %Clear graph on screen.
584 Chapter 13: Digital Control Systems
rlocus(Ge)
%Plot compensated root locus with %test lead compensator pole. axis([-5 2 -8 8]); %Change axes ranges. sgrid(z,wn) %Overlay grid on lead-compensated %root locus. title(['Lead-Compensated Root Locus with ' , num2str(pos),... '% Overshoot Line, Lead Pole at ', num2str(-a),... ' and Required Wn']) %Add title to lead-compensated root %locus. done=input('Are you done? (y=0,n=1) '); %Set loop flag. end %End loop for trying compensator %pole. [K,p]=rlocfind(Ge); %Generate gain, K, and closed-loop %poles, p, for point selected %interactively on the root locus. 'Gc(s)' %Display label. Gc=K*tf([1 b],[1 a]) %Display lead compensator. 'Gc(s)G(s)' %Display label. Ge %Display Gc(s)G(s). 'Closed-loop poles = ' %Display label. p %Display lead-compensated system's %closed-loop poles. f=input('Give pole number that is operating point '); %Choose lead-compensated system %dominant pole. 'Summary of estimated specifications for selected point on lead' 'compensated root locus' %Display label. operatingpoint=p(f) %Display lead-compensated dominant %pole. gain=K %Display lead-compensated gain. estimated_settling_time=4/abs(real(p(f))) %Display lead-compensated settling %time. estimated_peak_time=pi/abs(imag(p(f))) %Display lead-compensated peak time. estimated_percent_overshoot=pos %Display lead-compensated percent %overshoot. estimated_damping_ratio=z %Display lead-compensated damping %ratio. estimated_natural_frequency=sqrt(real(p(f))^2+imag(p(f))^2) %Display lead-compensated natural %frequency. s=tf([1 0],1); %Create transfer function, "s". sGe=s*Ge; %Create sGe(s) to evaluate Kv. sGe=minreal(sGe); %Cancel common poles and zeros. Kv=dcgain(K*sGe) %Display lead-compensated Kv. ess=1/Kv %Display lead-compensated steady%state error for unit ramp input. 'T(s)' %Display label. T=feedback(K*Ge,1) %Create and display lead-compensated %T(s). 'Press any key to continue and obtain the lead-compensated step' 'response' %Display label pause step(T) %Plot step response for lead %compensated system. title(['Lead-Compensated System with ' ,num2str(pos),'% Overshoot']) %Add title to step response of PD %compensated system. pause 'z-plane conversion for Challenge - Lead Comp' clf %Clear graph. 'Gc(s) in polynomial form' %Print label. Gcs=Gc %Create Gc(s) in polynomial form. 'Gc(s) in polynomial form' %Print label.
Solutions to Case Studies Challenges 585
Gcszpk=zpk(Gcs) %Create Gc(s) in factored form. 'Gc(z) in polynomial form via Tustin Transformation' %Print label. Gcz=c2d(Gcs,1/100,'tustin') %Form Gc(z) via Tustin %transformation. 'Gc(z) in factored form via Tustin Transformation' %Print label. Gczzpk=zpk(Gcz) %Show Gc(z) in factored form. 'Gp(s) in polynomial form' %Print label. Gps=G %Create Gp(s) in polynomial form. 'Gp(s) in factored form' %Print label. Gpszpk=zpk(Gps) %Create Gp(s) in factored form. 'Gp(z) in polynomial form' %Print label. Gpz=c2d(Gps,1/100,'zoh') %Form Gp(z) via zoh transformation. 'Gp(z) in factored form' %Print label. Gpzzpk=zpk(Gpz) %Form Gp(z) in factored form. pole(Gpz) %Find poles. Gez=Gcz*Gpz; %Form Ge(z) = Gc(z)Gp(z). 'Ge(z) = Gc(z)Gp(z) in factored form' %Print label. Gezzpk=zpk(Gez) %Form Ge(z) in factored form. 'z-1' %Print label. zm1=tf([1 -1],1,1/100) %Form z-1. zm1Gez=minreal(zm1*Gez,.00001); '(z-1)Ge(z)' %Print label. zm1Gezzpk=zpk(zm1Gez) pole(zm1Gez) Kv=300*dcgain(zm1Gez) Tz=feedback(Gez,1) step(Tz) title('Closed-Loop Digital Step Response') %Add title to step response.
Computer response:
ans = s-plane lead design for Challenge - Lead Comp ans = Uncompensated System ans = G(s) Zero/pole/gain: 0.16 ---------s (s+1.32) Type desired percent overshoot 10 Type Desired Peak Time 1 Type Lead Compensator Zero, (s+b). b= 1.32 Enter a Test Lead Compensator Pole, (s+a). a = Are you done? (y=0,n=1) 0 Select a point in the graphics window selected_point = -2.3045 + 3.1056i ans =
4.606
586 Chapter 13: Digital Control Systems
Gc(s) Transfer function: 93.43 s + 123.3 --------------s + 4.606 ans = Gc(s)G(s) Transfer function: 0.16 s + 0.2112 -----------------------s^3 + 5.926 s^2 + 6.08 s ans = Closed-loop poles = p = -2.3030 + 3.1056i -2.3030 - 3.1056i -1.3200 Give pole number that is operating point ans = Summary of estimated specifications for selected point on lead ans = compensated root locus operatingpoint = -2.3030 + 3.1056i gain = 93.4281 estimated_settling_time = 1.7369 estimated_peak_time = 1.0116 estimated_percent_overshoot = 10 1
Solutions to Case Studies Challenges 587
estimated_damping_ratio = 0.5912 estimated_natural_frequency = 3.8663 Kv = 3.2454 ess = 0.3081 ans = T(s) Transfer function: 14.95 s + 19.73 --------------------------------s^3 + 5.926 s^2 + 21.03 s + 19.73 ans = Press any key to continue and obtain the lead-compensated step ans = response ans = z-plane conversion for Challenge - Lead Comp ans = Gc(s) in polynomial form Transfer function: 93.43 s + 123.3 --------------s + 4.606 ans = Gc(s) in polynomial form Zero/pole/gain: 93.4281 (s+1.32) ----------------
588 Chapter 13: Digital Control Systems
(s+4.606) ans = Gc(z) in polynomial form via Tustin Transformation Transfer function: 91.93 z - 90.72 --------------z - 0.955 Sampling time: 0.01 ans = Gc(z) in factored form via Tustin Transformation Zero/pole/gain: 91.9277 (z-0.9869) -----------------(z-0.955) Sampling time: 0.01 ans = Gp(s) in polynomial form Transfer function: 0.16 -----------s^2 + 1.32 s ans = Gp(s) in factored form Zero/pole/gain: 0.16 ---------s (s+1.32) ans = Gp(z) in polynomial form Transfer function: 7.965e-006 z + 7.93e-006 -----------------------z^2 - 1.987 z + 0.9869 Sampling time: 0.01 ans = Gp(z) in factored form Zero/pole/gain:
Solutions to Case Studies Challenges 589
7.9649e-006 (z+0.9956) ---------------------(z-1) (z-0.9869) Sampling time: 0.01 ans = 1.0000 0.9869 ans = Ge(z) = Gc(z)Gp(z) in factored form Zero/pole/gain: 0.0007322 (z+0.9956) (z-0.9869) ------------------------------(z-1) (z-0.9869) (z-0.955) Sampling time: 0.01 ans = z-1 Transfer function: z - 1 Sampling time: 0.01 ans = (z-1)Ge(z) Zero/pole/gain: 0.0007322 (z+0.9956) -------------------(z-0.955) Sampling time: 0.01 ans = 0.9550 Kv = 9.7362 Transfer function: 0.0007322 z^2 + 6.387e-006 z - 0.0007194 ---------------------------------------z^3 - 2.941 z^2 + 2.884 z - 0.9432 Sampling time: 0.01
590 Chapter 13: Digital Control Systems
Solutions to Case Studies Challenges 591
592 Chapter 13: Digital Control Systems
ANSWERS TO REVIEW QUESTIONS
1. (1) Supervisory functions external to the loop; (2) controller functions in the loop 2. (1) Control of multiple loops by the same hardware; (2) modifications made with software, not hardware; (3) more noise immunity (4) large gains usually not required 3. Quantization error; conversion time 4. An ideal sampler followed by a sample-and-hold 5. z = esT 6. The value of the time waveform only at the sampling instants 7. Partial fraction expansion; division to yield power series 8. Partial fraction 9. Division to yield power series 10. The input must be sampled; the output must be either sampled or thought of as sampled. 11. c(t) is c*(t) = c(kT), i.e. the output only at the sampling instants. 12. No; the waveform is only valid at the sampling instants. Instability may be apparent if one could only see between the sampling instants. The roots of the denominator of G(z) must be checked to see that they are within the unit circle. 13. A sample-and-hold must be present between the cascaded systems. 14. Inside the unit circle 15. Raible table; Jury's stability test 16. z=+1 17. There is no difference. 18. Map the point back to the s-plane. Since z = esT, s = (1/T) ln z. Thus, = (1/T) ln (Re z), and = (1/T) ln (Im z). 19. Determine the point on the s-plane and use z = esT. Thus, Re z = eT cos , and Im z = eT sin . 20. Use the techniques described in Chapters 9 and 11 and then convert the design to a digital compensator using the Tustin transformation. 21. Both compensators yield the same output at the sampling instants.
SOLUTIONS TO PROBLEMS
1. a. f(t) = e-at; f*(t) =
e-akT(t-kT) ; F*(s) =
k=0
e-akT e-kTs
= 1 + e-aT e-Ts + e-a2T e-2Ts + . . . Thus,
k=0
F(z) = 1 + e-aT z-1 + e-a2T z-2 + . . . = 1 + x-1 + x-2 + . . . where x = e-aT z-1.
Solutions to Problems 593 1 1 z = = . 1 - x-1 1 - e-aT z-1 z - e-aT b. f(t) = u(t); f*(t) = (t-kT) ; F*(s) = e-kT s = 1 + e-Ts + e-2Ts + . . . k=0 k=0 But, F(z) = Thus, F(z) = 1 + z-1 + z-2 + . . z z-1 . c. f(t) = t2 e-at; f*(t) = Since 1 = 1 + z-1 + z-2 + z-3, F(z) = 1 - z-1
z-k
k=0
=
1 = 1 - z-1
(kT)2e-akT(t-kT) ; F*(s) = T2 k2 e-akT e-kTs
k=0 k=0 = T2
= T2
k2 (e-(s+a)T)k
k=0
k2 xk
k=0
= T2(x + 4x2 + 9x3 + 16x4 + . . . ) , where x = e- (s+a)T.
Let s1 = x + 4x2 + 9x3 + 16x4 + . . . Thus, xs1 = x2 + 4x3 + 9x4 + 16x5 + . . . Let s2 = s1 - xs1 = x + 3x2 + 5x3 + 7x4 + . . .Thus, xs2 = x2 + 4x3 + 9x4 + 16x3 + . . . Let s3 = s2 - xs2 = x + 2x2 + 2x3 + 2x4 + . . . Thus xs3 = x2 + 2x3 + 2x4 + 2x3 + . . . Let s4 = s3 - xs3 = x + x2. Solving for s3, x + x2 s3 = 1 - x and s3 x + x2 s2 = 1 - x = (1 - x)2 and s2 x + x2 s1 = 1 - x = (1 - x)3 Thus F*(s) = T2 s1 = T2
-(s+a)T+ e-2(s+a)T) x + x2 2 (e = 3 =T (1 - e-(s+a)T)3 (1 - x)
T 2 [z -1 e -aT + z -2 e -2aT ] T 2 ze -aT [z + e - aT ] = - aT 3 -3 -aT 3 (z - e ) z (z - e )
d. f(t) = cos(kT); f*(t) =
cos(kT) (t
k =0
- kT) ; F*(s) =
cos(kT)e
k =0
- kTs
=
(e jkT + e - jkT )e -kTs 1 T (s - j ) - k T (s + j ) - k ) +(e ) = (e 2 2 k=0 k =0
But,
594 Chapter 13: Digital Control Systems
x-k
k=0 Thus, F*(s) =
=
1 . 1 - x-1
1 1 1 2 - e -Ts ( e jT + e- jT ) = 1 + -T (s - j ) -T (s+ j ) -T (s - j ) -T (s + j ) -T (s- j ) - T (s + j ) 2 1- e -e )+e e 2 1- e 1-e
=
1 2 - e - Ts (2cos(T ) 1 - z -1 cos(T ) = 2 1- e -Ts (e jT + e - jT ) + e -2Ts 1 - 2z -1 cos(T ) + z -2
Therefore,
F( z ) =
2. Program:
syms T a w n '(a)' 'f(kT)' f=exp(-a*n*T); pretty(f) 'F(z)' F=ztrans(f); pretty(F) '(b)' 'f(kT)' f=exp(-0*n*T); pretty(f) 'F(z)' F=ztrans(f); pretty(F) '(c)' 'f(kT)' f=(n*T)^2*exp(-a*n*T); pretty(f) 'F(z)' F=ztrans(f); pretty(F) '(d)' 'f(kT)' f=cos(w*n*T); pretty(f) 'F(z)' F=ztrans(f); pretty(F)
z (z - cos( T)) z - 2 z cos(T ) + 1
2
%Construct symbolic objects for %'T', 'a','w', and 'n'. %Display label. %Display label. %Define f(kT). %Pretty print f(kT) %Display label. %Find z-transform, F(z). %Pretty print F(z). %Display label. %Display label. %Define f(kT) %Pretty print f(kT) %Display label. %Find z-transform, F(z). %Pretty print F(z). %Display label. %Display label. %Define f(kT) %Pretty print f(kT) %Display label. %Find z-transform, F(z). %Pretty print F(z). %Display label. %Display label. %Define f(kT) %Pretty print f(kT) %Display label. %Find z-transform, F(z). %Pretty print F(z).
Computer response:
ans = (a) ans = f(kT) exp(-a n T)
Solutions to Problems 595
ans = F(z) z ------------------------/ z \ exp(-a T) |--------- - 1| \exp(-a T) / ans = (b) ans = f(kT) 1 ans =
F(z) z ----z - 1 ans = (c) ans = f(kT) 2 2 n T exp(-a n T) ans = F(z) 2 T z exp(-a T) (z + exp(-a T)) -----------------------------3 (z - exp(-a T)) ans = (d) ans = f(kT) cos(w n T) ans = F(z) (z - cos(w T)) z --------------------2 z - 2 z cos(w T) + 1
596 Chapter 13: Digital Control Systems
3. a.
F( z ) =
z( z + 3)( z + 5) ( z - 0.4)(z - 0.6)(z - 0.8)
275.5 504 F(z) 229.5 + - = z - 0.4 z - 0.6 z - 0.8 z
F(z) =
275.5z 229.5z 504z + - z - 0.4 z - 0.6 z - 0.8
k k k
f (kT) = 229.5(0.4) - 504(0.6) + 275.5(0.8) , k = 0, 1, 2, 3,. . .
b.
F( z) =
(z + 0.2)( z + 0.4) ( z - 0.1)(z - 0.5)(z - 0.9)
F(z) 1.778 4.6875 7.875 4.9653 + - + =- z - 0.1 z - 0.5 z - 0.9 z z
F(z) = -1.778 + 4.6875z 7.875 z 4.9653 z - + z - 0.1 z - 0.5 z - 0.9
k k k
f (kT) = 4.6875(0.1) - 7.875(0.5) + 4.9653(0.9) , k = 1 ,2, 3,. . .
c.
( z + 1)( z + 0.3)( z + 0.4) z ( z - 0.2)( z - 0.5)( z - 0.7) F ( z) ( z + 1)( z + 0.3)( z + 0.4) = 2 z z ( z - 0.2)( z - 0.5)( z - 0.7) 38.1633 72 60 26.1633 1.7143 = - + - - z - 0.7 z - 0.5 z - 0.2 z z2 38.1633z 72 z 60 z 1.7143 F ( z) = - + - 26.1633 - z z - 0.7 z - 0.5 z - 0.2 k k k F = 38.1633(0.7) - 72(0.5) + 60(0.2) for k = 2,3, 4,... =1 for k = 1 =0 for k = 0
F ( z) =
4. Program:
'(a)' syms z k F=vpa(z*(z+3)*(z+5)/((z-0.4)*(z-0.6)*(z-0.8)),4); pretty(F) f=vpa(iztrans(F),4); pretty(f) '(b)' syms z k F=vpa((z+0.2)*(z+0.4)/((z-0.1)*(z-0.5)*(z-0.9)),4); pretty(F)
Solutions to Problems 597
f=vpa(iztrans(F),4); pretty(f) '(c)' syms z k F=vpa((z+1)*(z+0.3)*(z+0.4)/(z*(z-0.2)*(z-0.5)*(z-0.7)),4); pretty(F) f=vpa(iztrans(F),4); pretty(f)
Computer response:
ans = (a) z (z + 3.) (z + 5.) ----------------------------------(z - .4000) (z - .6000) (z - .8000) n n n 275.5 .8000 - 504.0 .6000 + 229.5 .4000 ans = (b) (z + .2000) (z + .4000) ----------------------------------(z - .1000) (z - .5000) (z - .9000) n n n -1.778 charfcn[0](n) + 4.965 .9000 - 7.875 .5000 + 4.688 .1000 ans = (c) (z + 1.) (z + .3000) (z + .4000) ------------------------------------z (z - .2000) (z - .5000) (z - .7000) n n -1.714 charfcn[1](n) - 26.16 charfcn[0](n) + 38.16 .7000 - 72.00 .5000 n + 60.00 .2000
598 Chapter 13: Digital Control Systems
5. a.
By division Instant 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Value 1 9.8 31.6 46.88 53.4016 53.43488 49.64608 44.043776 37.90637056 31.95798733 26.5581568 21.84639857 17.83896791 14.48905384 11.72227881 9.456567702 7.612550239 6.118437551 4.911796342 3.939668009 3.15787423 2.529983782 2.026197867 1.622284879 1.298623886 1.039376712 0.831787937 0.665602292 0.532584999 0.4261299 0.34094106 k 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 By Formula Value 1 9.8 31.6 46.88 53.4016 53.43488 49.64608 44.043776 37.90637056 31.95798733 26.5581568 21.84639857 17.83896791 14.48905384 11.72227881 9.456567702 7.612550239 6.118437551 4.911796342 3.939668009 3.15787423 2.529983782 2.026197867 1.622284879 1.298623886 1.039376712 0.831787937 0.665602292 0.532584999 0.4261299 0.34094106
Solutions to Problems 599
b.
By division Instant 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 Value 1 2.1 2.64 2.766 2.6859 2.51571 2.313354 2.1066276 1.90826949 1.723594881 1.554311564 1.400418494 1.261145687 1.13541564 1.022066337 0.919955834 0.828008315 0.745231516 0.670720381 0.603654351 0.54329192 0.48896423 0.440068558 0.396062078 0.356456058 0.320810546 0.288729538 0.259856608 0.233870959 0.210483869 0.189435485 k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 By Formula Value 1.00002 2.100018 2.6400162 2.76601458 2.685913122 2.51572181 2.313364629 2.106637166 1.908278099 1.723602629 1.554318538 1.40042477 1.261151336 1.135420724 1.022070912 0.919959951 0.828012021 0.745234852 0.670723383 0.603657053 0.543294352 0.488966419 0.440070528 0.39606385 0.356457653 0.320811982 0.28873083 0.259857771 0.233872006 0.210484811 0.189436333
600 Chapter 13: Digital Control Systems
c.
6. a.
G (s ) =
(s + 4) 0.6667 0.3333 = + ( s + 2)( s + 5) s+2 s +5
G(z) = 0.6667z 0.3333z + -5T z - e -2T z - e
For T = 0.5 s,
G( z) =
b.
0.6667z 0.3333z z(z - 0.1774) + = z - 0.3679 z - 0.082085 ( z - 0.3679)( z - 0.082085) (s + 1)(s + 2) 0.1667 0.6667 1.5 = - + s (s + 3)(s + 4) s s +3 s +4
0.1667z 0.6667z 1.5z - -3T + z-1 z -e z - e-4T
G (s ) =
G(z) =
For T = 0.5 s,
G( z) =
0.1667z 0.6667 z 1.5 z z(z - 0.29675)(z - 0.8408) - + = z-1 z - 0.22313 z - 0.13534 ( z - 1)(z - 0.22313)( z - 0.13534)
Solutions to Problems 601
c.
G(s) =
20 1.25 1.25s + 3.57 1.25 1.25( s + 3) = - = - ( s + 3)( s + 6 s + 25) s + 3 s 2 + 6s + 25 s + 3 ( s + 3) 2 + 42 z z 2 - zae - aT cos T G ( z ) = -1.25 - 1.25 2 z - e - aT z - 2 ze - aT cos T + e -2 aT
2
For a = 3;
= 4; T = 0.5 , z z 2 + 0.0929 z G ( z ) = -1.25 - 1.25 2 z - 0.2231 z + 0.1857 z + 0.0498 z ( z + 0.2232) = 0.395 ( z - 0.2231)( z 2 + 0.1857 z + 0.0498)
15 0.1852 0.2083 s + 0.9978 = - + 0.02314 2 s ( s + 1)( s + 10s + 81) s s +1 s + 10s + 81
2
d.
G ( s) =
0.1852 0.2083 ( s + 5) - 0.5348 56 - + 0.02314 s s +1 ( s + 5) 2 + 56 z z z 2 - zae - aT cos T ze - aT sin T G ( z ) = 0.1852 - 0.2083 + 0.02314 2 - 0.0124 2 z -1 z - eT z - 2 ze - aT cos T + e -2 aT z - 2 ze- aT cos T + e -2 aT =
For a = 5;
=1; = 56; T = 0.5 ,
G ( z ) = 0.1852
z z z 2 + 0.0678 z z - 0.2083 + 0.02314 2 + 0.0005748 2 z -1 z - 0.6065 z + 0.1355 z + 0.006738 z + 0.1355 z + 0.006738
=
0.00004 z 4 + 0.05781z 3 + 0.02344 z 2 + 0.001946 z 0.05781z 3 + 0.02344 z 2 + 0.001946 z ( z - 1)( z - 0.6065) z 2 + 0.1355 z + 0.006738 ( z - 1)( z - 0.6065) z 2 + 0.1355 z + 0.006738
(
)
(
)
)
= 0.05781
z 3 + 0.4055 z 2 + 0.0337 z z ( z + 0.2888)( z + 0.1167) = 0.05781 2 ( z - 1)( z - 0.6065) z + 0.1355 z + 0.006738 ( z - 1)( z - 0.6065) z 2 + 0.1355 z + 0.006738
(
)
(
7. Program:
'(a)' syms s z n T Gs=(s+4)/((s+2)*(s+5)); 'G(s)' pretty(Gs) %'g(t)' gt=ilaplace(Gs); %pretty(gt) gnT=compose(gt,n*T); %'g(kT)' %pretty(gnT) Gz=ztrans(gnT); Gz=simplify(Gz); %Construct symbolic objects for %'s', 'z', 'n',and 'T'. %Form G(s). %Display label. %Pretty print G(s). %Display label. %Find g(t). %Pretty print g(t). %Find g(nT). %Display label. %Pretty print g(nT). %Find G(z). %Simplify G(z).
602 Chapter 13: Digital Control Systems
%'G(z)' %pretty(Gz) Gz=subs(Gz,T,0.5); Gz=vpa(simplify(Gz),6); Gz=vpa(factor(Gz),6); 'G(z) evaluated for T=0.5' pretty(Gz)
%Display label. %Pretty print G(z). %Let T = 0.5 in G(z). %Simplify G(z) and evaluate numerical %values to 6 places. %Factor G(z). %Display label. %Pretty print G(z) with numerical %values.
'(b)' Gs=(s+1)*(s+2)/(s*(s+3)*(s+4)); %Form G(s) = G(s). 'G(s)' %Display label. pretty(Gs) %Pretty print G(s). %'g(t)' %Display label. gt=ilaplace(Gs); %Find g(t). %pretty(gt) %Pretty print g(t). gnT=compose(gt,n*T); %Find g(nT). %'g(kT)' %Display label. %pretty(gnT) %Pretty print g(nT). Gz=ztrans(gnT); %Find G(z). Gz=simplify(Gz); %Simplify G(z). %'G(z)' %Display label. %pretty(Gz) %Pretty print G(z). Gz=subs(Gz,T,0.5); %Let T = 0.5 in G(z). Gz=vpa(simplify(Gz),6); %Simplify G(z) and evaluate numerical %values to 6 places. Gz=vpa(factor(Gz),6); %Factor G(z). 'G(z) evaluated for T=0.5' %Display label. pretty(Gz) %Pretty print G(z) with numerical %values. '(c)' Gs=20/((s+3)*(s^2+6*s+25)); %Form G(s) = G(s). 'G(s)' %Display label. pretty(Gs) %Pretty print G(s). %'g(t)' %Display label. gt=ilaplace(Gs); %Find g(t). %pretty(gt) %Pretty print g(t). gnT=compose(gt,n*T); %Find g(nT). %'g(kT)' %Display label. %pretty(gnT) %Pretty print g(nT). Gz=ztrans(gnT); %Find G(z). Gz=simplify(Gz); %Simplify G(z). %'G(z)' %Display label. %pretty(Gz) %Pretty print G(z). Gz=subs(Gz,T,0.5); %Let T = 0.5 in G(z). Gz=vpa(simplify(Gz),6); %Simplify G(z) and evaluate numerical %values to 6 places. Gz=vpa(factor(Gz),6); %Factor G(z). 'G(z) evaluated for T=0.5' %Display label. pretty(Gz) %Pretty print G(z) with numerical %values. '(d)' Gs=15/(s*(s+1)*(s^2+10*s+81)); %Form G(s) = G(s). 'G(s)' %Display label. pretty(Gs) %Pretty print G(s). %'g(t)' %Display label. gt=ilaplace(Gs); %Find g(t). %pretty(gt) %Pretty print g(t). gnT=compose(gt,n*T); %Find g(nT). %'g(kT)' %Display label. %pretty(gnT) %Pretty print g(nT). Gz=ztrans(gnT); %Find G(z). Gz=simplify(Gz); %Simplify G(z). %'G(z)' %Display label.
Solutions to Problems 603
%pretty(Gz) Gz=subs(Gz,T,0.5); Gz=vpa(simplify(Gz),6); Gz=vpa(factor(Gz),6); 'G(z) evaluated for T=0.5' pretty(Gz)
%Pretty print G(z). %Let T = 0.5 in G(z). %Simplify G(z) and evaluate numerical %values to 6 places. %Factor G(z). %Display label. %Pretty print G(z) with numerical %values.
Computer response:
ans = (a) ans = G(s) s + 4 --------------(s + 2) (s + 5) ans = G(z) evaluated for T=0.5 z (z - .177350) 1.00000 ---------------------------(z - .0820850) (z - .367880) ans = (b) ans = G(s) (s + 1) (s + 2) ----------------s (s + 3) (s + 4) ans = G(z) evaluated for T=0.5 z (z - .296742) (z - .840812) 1.00000 -----------------------------------(z - .135335) (z - .223130) (z - 1.) ans = ans = (c) ans = G(s) 20 ----------------------2 (s + 3) (s + 6 s + 25) ans = G(z) evaluated for T=0.5
604 Chapter 13: Digital Control Systems
(z + .223130) z .394980 ----------------------------------------2 (z - .223135) (z + .185705 z + .0497861) ans = (d) ans = G(s) 15 -------------------------2 s (s + 1) (s + 10 s + 81) ans = G(z) evaluated for T=0.5 (z + .289175) (z + .116364) z .0578297 -------------------------------------------------------2 (z - .606535) (z - .999995) (z + .135489 z + .00673794)
8. a.
Thus, Hence,
Letting T = 0.3,
b.
Solutions to Problems 605
Thus,
Hence,
Letting T = 0.3,
c.
Ge ( z ) = Ga ( z )G ( z )
where Ga ( z ) is the answer to part (a) and G ( z ) , the pulse transfer function for with a zero-order-hold will now be found:
1 in cascade s+3
Thus,
Hence,
606 Chapter 13: Digital Control Systems
Letting T = 0.3,
Thus,
Ge ( z ) = Ga ( z )G ( z ) = 0.19464
9.
6.482 z + 3.964 ( z - 1)(4.482 z - 1)(2.46 z - 1)
a. Add phantom samplers at the input, output, and feedback path after H(s). Push G2 (s) and its input sampler to the right past the pickoff point. Add a phantom sampler after G1(s). Hence,
R(s)
+ -
C(s) G (s)
1
G (s)
2
H (s) G (s)
2
G1(z)G2(z) From this block diagram, T(z) = 1 + G (z)HG (z) .
1 2
b. Add phantom samplers to the input, output, and the output of H(s). Push G1(s)G2(s) and its input sampler to the right past the pickoff point. Add a phantom sampler at the output.
R(s)
+
C(s)
G (s) G (s)
1 2
H(s)
G (s) G (s) 1
2
Solutions to Problems 607 G1G2(z) From this block diagram, T(z) = 1 + G G (z)H(z) .
1 2
10. Add phantom samplers after G1(s), G3(s), G4(s), H1(s), and H2(s).
Push G1(s) and its sampler to the left past the summing junction. Also, push G4(s) and its input sampler to the right past the pickoff point. The resulting block diagram is,
Converting to z transforms,
608 Chapter 13: Digital Control Systems
G2G3 ( z ) 1 * (1 + G2G3 ( z ) H1 ( z )) (1 + H 2G4 ( z )) C ( s) = RG1 ( z )G4 ( z ) 1 G2G3 ( z ) 1 + * H G G ( z) (1 + G2G3 ( z ) H1 ( z )) (1 + H 2G4 ( z )) 3 1 4
=
11. Push gain of 2 to the left past the summing junction and add phantom samplers as shown.
RG1 ( z )G4 ( z )G2G3 ( z ) (1 + G2G3 ( z ) H1 ( z ))(1 + H 2G4 ( z )) + G2G3 ( z ) H 3G1G4 ( z )
Push the z.o.h. output.
1 and to the right past the pickoff point. Also, add a phantom sampler at the s(s+1)
Add phantom samplers after the gain of 2 at the input and in the feedback. Also, represent the z.o.h. as Laplace transforms.
Solutions to Problems 609
Push the last block to the right past the pickoff point and get,
Find the z transform for each transfer function.
G1 ( s ) = 2
transforms into
G1 ( z ) = 2 .
H1 ( s ) = (1 - e - sT )
transforms into
2 2 2 2 = (1 - e - sT ) 2 - + s ( s + 1) s s + 1 s
2
z -1 Tz z z Tz - Te -T + ze -T - z - e-T + 1 H1 ( z ) = 2 -2 +2 =2 z ( z - 1) 2 z -1 z - e -T ( z - 1)( z - e -T )
H 2 ( s ) = (1 - e - sT )
transforms into
2 2 2 2 2 = (1 - e - sT ) - - 2 + 3 s ( s + 1) s s s +1 s
3
H 2 ( z) =
z -1 2z 2z 2Tz T 2 z ( z + 1) - - + z z - 1 z - e -T ( z - 1) 2 ( z - 1)3
610 Chapter 13: Digital Control Systems
=
(T 2 - 2e-T + 2 - 2T ) z 2 + (4e-T - 4 + 2Te-T + 2T + T 2 - T 2 e-T ) z + (2 - 2e-T - 2Te -T - T 2 e-T ) ( z - 1)2 ( z - e-T )
G2 ( s) = (1 - e - sT )
transforms into
1 s ( s + 1)
3
1 H 2 ( z) 2
Thus, the closed-loop transfer function is
1 T ( z ) = G1 ( z )G2 ( z ) 1 + H1 ( z ) + H 2 ( z )
12. z-1 G(z) = z Using Eq. (13.49) T (1 - e-T) (T-1+e-T)z + (1-e-T-Te-T) G(z) = z - 1 -T = z-e (z-1)(z-e-T) But, G(z) (T-1+e-T) z + (1-e-T-Te-T) T(z) = 1 + G(z) = 2 z + (T-2) z + (1-Te-T) The roots of the denominator are inside the unit circle for 0<T<3.923. 13. Program:
numg1=10*[1 7]; deng1=poly([-1 -3 -4 -5]); G1=tf(numg1,deng1); for T=5:-.01:0; Gz=c2d(G1,T,'zoh'); Tz=feedback(Gz,1); r=pole(Tz); rm=max(abs(r)); if rm<=1; break; end; end; T r rm
z s2(s1+ 1) .
Computer response:
T = 3.3600 r = -0.9990 -0.0461 -0.0001
Solutions to Problems 611
-0.0000 rm = 0.9990 >> T = 3.3600 r = -0.9990 -0.0461 -0.0001 -0.0000 rm = 0.9990
14.
612 Chapter 13: Digital Control Systems
15. a.
First, check to see that the system is stable.
Solutions to Problems 613 Since the closed-loop poles are inside the unit circle, the system is stable. Next, evaluate the static error constants and the steady-state error.
K p = lim G(z) = 0.5
z 1
e * () =
e () =
*
1 2 = 1+ K p 3
1 = Kv
1 = Ka
Kv =
Ka =
1 lim(z - 1)G(z) = 0 T z 1
1 2 2 lim(z - 1) G(z) = 0 T z1
e () =
*
b.
K p = lim G(z) =
z 1
e* () =
*
1 =0 1 + Kp
1 = 0.1 Kv
1 = Ka
Kv =
Ka =
1 lim(z - 1)G(z) = 10 T z 1
1 2 2 lim(z - 1) G(z) = 0 z1 T
e () =
e () =
*
614 Chapter 13: Digital Control Systems
c.
K p = lim G(z) = 2.03
z 1
e * () =
e () =
*
1 = 0.33 1+ K p
1 = Kv
1 = Ka
Kv =
Ka =
d.
1 lim(z - 1)G(z) = 0 T z 1
1 2 2 lim(z - 1) G(z) = 0 T z1
e () =
*
K p = lim G(z) =
z 1
e* () =
1 =0 1 + Kp
1 1 * lim(z - 1)G(z) = 10 e () = = 0.1 T z 1 Kv 1 1 2 * Ka = 2 lim(z - 1) G(z) = 0 e () = = z1 T Ka Kv =
16. Program:
T=0.1; numgz=[0.04406 -0.03624 -0.03284 0.02857]; dengz=[1 -3.394 +4.29 -2.393 +0.4966]; 'G(z)' Gz=tf(numgz,dengz,0.1) 'Zeros of G(z)' zeros=roots(numgz) 'Poles of G(z)'
Solutions to Problems 615
poles=roots(dengz) %Check stability Tz=feedback(Gz,1); 'Closed-Loop Poles' r=pole(Tz) M=abs(r) pause 'Find Kp' Gz=minreal(Gz,.00001); Kp=dcgain(Gz) 'Find Kv' factorkv=tf([1 -1],[1 0],0.1); Gzkv=factorkv*Gz; Gzkv=minreal(Gzkv,.00001);
%Makes transfer function %proper and yields same Kv %Cancel common poles and %zeros
Kv=(1/T)*dcgain(Gzkv) 'Find Ka' factorka=tf([1 -2 1],[1 0 0],0.1);%Makes transfer function %proper and yields same Ka Gzka=factorka*Gz; Gzka=minreal(Gzka,.00001); Ka=(1/T)^2*dcgain(Gzka) %Cancel common poles and %zeros
Computer response:
ans = G(z) Transfer function: 0.04406 z^3 - 0.03624 z^2 - 0.03284 z + 0.02857 ----------------------------------------------z^4 - 3.394 z^3 + 4.29 z^2 - 2.393 z + 0.4966 Sampling time: 0.1 ans = Zeros of G(z) zeros = -0.8753 0.8489 + 0.1419i 0.8489 - 0.1419i ans = Poles of G(z) poles = 1.0392 0.8496 + 0.0839i 0.8496 - 0.0839i 0.6557 ans =
616 Chapter 13: Digital Control Systems
Closed-Loop Poles r = 0.9176 0.9176 0.7573 0.7573 M = 0.9332 0.9332 0.7765 0.7765 ans = Find Kp Kp = -8.8750 ans = Find Kv Kv = 0 ans = Find Ka Ka = 0 + + 0.1699i 0.1699i 0.1716i 0.1716i
17. First find G(z)
Solutions to Problems 617 Next, plot the root locus.
18.
9.68 Step Sum Zero-Order Hold Gain
1 2 +4s+3 s Transfer Fcn Scope
618 Chapter 13: Digital Control Systems
19. Program:
numgz=0.13*[1 1]; dengz=poly([1 0.74]); Gz=tf(numgz,dengz,0.1) Gzpkz=zpk(Gz) Tz=feedback(Gz,1) Ltiview
Computer response:
Transfer function: 0.13 z + 0.13 ------------------z^2 - 1.74 z + 0.74 Sampling time: 0.1
Zero/pole/gain: 0.13 (z+1) -------------(z-1) (z-0.74) Sampling time: 0.1 Transfer function: 0.13 z + 0.13 ------------------z^2 - 1.61 z + 0.87 Sampling time: 0.1
Solutions to Problems 619
20. Program:
%Digitize G1(s) preceded by a sample and hold numg1=1; deng1=poly([-1 -3]); 'G1(s)' G1s=tf(numg1,deng1) 'G(z)' Gz=c2d(G1s,0.1,'zoh') %Input transient response specifications Po=input('Type %OS '); %Determine damping ratio z=(-log(Po/100))/(sqrt(pi^2+log(Po/100)^2)); %Plot root locus rlocus(Gz) zgrid(z,0) title(['Root Locus']) [K,p]=rlocfind(Gz) %Allows input by selecting point on graphic pause 'T(z)' Tz=feedback(K*Gz,1) step(Tz)
Computer response:
ans = G1(s) Transfer function:
620 Chapter 13: Digital Control Systems
1 ------------s^2 + 4 s + 3 ans = G(z) Transfer function: 0.004384 z + 0.003837 ---------------------z^2 - 1.646 z + 0.6703 Sampling time: 0.1 Type %OS 16.3 Select a point in the graphics window selected_point = 0.8016 + 0.2553i K = 9.7200 p = 0.8015 + 0.2553i 0.8015 - 0.2553i ans = T(z) Transfer function: 0.04262 z + 0.0373 ---------------------z^2 - 1.603 z + 0.7076 Sampling time: 0.1
Solutions to Problems 621
622 Chapter 13: Digital Control Systems
21.
13.12
25.46 deg 9 deg
Solutions to Problems 623
22. z-1 1 z-1 First find G(z). G(z) = K z z{ 2 }=K z s (s+1)(s+3) z{ - 1
18
1 +1 1 -4 1 +1 1 } s +3 2 s +1 9 s 3 s 2
1 z 4z 1z - 1 z z -1 30 9 2 18 + + - For T=0.1, G(z) = K z - 0.74082 z - 0.90484 z - 1 z - 1 2 z z + 0.24204 z + 3.3828 = 0.00015103 K . Plotting the root locus and overlaying z - 1 z - 0.74082 z - 0.90484
the 20% overshoot curve, we select the point of intersection and calculate the gain: 0.00015103K = 1.789. Thus, K = 11845.33. Finding the intersection with the unit circle yields 0.00015103K = 9.85. Thus, 0 < K < 65218.83 for stability.
23. z-1 (s+2) z-1 First find G(z). G(z) = K z z{ 2 } = K z z{ 1 1 + 1 1 - 5 1 + 2 1 } = 18 s (s+1)(s+3) s +3 2 s +1 9 s 3 s 2
1 z 1 z 5 z 2 z 18 2 - 9 + + 3 For T=1, G(z) = K z - 1 z - 0.049787 z - 0.36788 z - 1 z - 1 z
2
= 0.29782K
z - 0.13774 z + 0.55935 . Plotting the root locus and overlaying the Ts = 15 z - 1 z - 0.049787 z - 0.36788
second circle, we select the point of intersection (0.4 + j0.63) and calculate the gain: 0.29782K = 1.6881. Thus, K = 5.668. Finding the intersection with the unit circle yields 0.29782K = 4.4. Thus, 0 < K < 14.77 for stability.
624 Chapter 13: Digital Control Systems
Root Locus 1
0.5
0
-0.5
-1
-1
-0.5
0 Real Axis
0.5
1
24. Substituting Eq. (13.88) into Gc(s) and letting T = 0.01 yields
(z - 0.9959)(z - 0.5625) 1180z 2 - 1839z + 661.1 = 1180 Gc (z) = 2 (z + 1)(z - 1) z -1
25. Since %OS = 10%, = %OS - ln ( 100 ) 4 = 0.591. Since Ts = = 2 seconds, n %OS 2 + ln2 ( 100 )
n = 3.383 rad/s. Hence, the location of the closed-loop poles must be 2 j2.729. The summation of angles from open-loop poles to 2 j2.729 is 192.99o. Therefore, the design point is not on the root locus. A compensator whose angular contribution is 192.99o - 180o = 12.99o is required. Assume a compensator zero at -5 canceling the pole at -5. Adding the compensator zero at 5 to the plant's poles yields 150.690 at to 2 j2.729. Thus, the compensator's pole must contribute 1800 150.69 = 29.310. The following geometry results.
j2.729
29.31 -2
Thus,
2.729 = tan(29.310 ) pc - 2
Solutions to Problems 625
Hence, pc = 6.86.Adding the compensator pole and zero to the system poles, the gain at the design point is found to be 124.3. Summarizing the results: Eq. (13.88) into Gc(s) and letting T = 0.01 yields
Gc (s) =
124.3(s + 5) . Substituting (s + 6.86)
Gc (z) =
26. Program:
123.2z - 117.2 123.2(z - 0.9512) = z - 0.9337 (z - 0.9337)
'Design of digital lead compensation' clf %Clear graph on screen. 'Uncompensated System' %Display label. numg=1; %Generate numerator of G(s). deng=poly([0 -5 -8]); %Generate denominator of G(s). 'G(s)' %Display label. G=tf(numg,deng) %Create and display G(s). pos=input('Type desired percent overshoot '); %Input desired percent overshoot. z=-log(pos/100)/sqrt(pi^2+[log(pos/100)]^2); %Calculate damping ratio. rlocus(G) %Plot uncompensated root locus. sgrid(z,0) %Overlay desired percent overshoot %line. title(['Uncompensated Root Locus with ' , num2str(pos),... '% Overshoot Line']) %Title uncompensated root locus. [K,p]=rlocfind(G); %Generate gain, K, and closed-loop %poles, p, for point selected %interactively on the root locus. 'Closed-loop poles = ' %Display label. p %Display closed-loop poles. f=input('Give pole number that is operating point '); %Choose uncompensated system %dominant pole. 'Summary of estimated specifications for selected point on' 'uncompensated root locus' %Display label. operatingpoint=p(f) %Display uncompensated dominant %pole. gain=K %Display uncompensated gain. estimated_settling_time=4/abs(real(p(f))) %Display uncompensated settling %time. estimated_peak_time=pi/abs(imag(p(f))) %Display uncompensated peak time. estimated_percent_overshoot=pos %Display uncompensated percent %overshoot. estimated_damping_ratio=z %Display uncompensated damping %ratio. estimated_natural_frequency=sqrt(real(p(f))^2+imag(p(f))^2) %Display uncompensated natural %frequency. numkv=conv([1 0],numg); %Set up numerator to evaluate Kv. denkv=deng; %Set up denominator to evaluate Kv. sG=tf(numkv,denkv); %Create sG(s). sG=minreal(sG); %Cancel common poles and zeros. Kv=dcgain(K*sG) %Display uncompensated Kv. ess=1/Kv %Display uncompensated steady-state %error for unit ramp input. 'T(s)' %Display label. T=feedback(K*G,1) %Create and display T(s). step(T) %Plot step response of uncompensated %system. title(['Uncompensated System with ' ,num2str(pos),'% Overshoot']) %Add title to uncompensated step %response.
626 Chapter 13: Digital Control Systems
`Press any key to go to lead compensation' %Display label. pause Ts=input('Type Desired Settling Time '); %Input desired settling time. b=input('Type Lead Compensator Zero, (s+b). b= '); %Input lead compensator zero. done=1; %Set loop flag. while done==1 %Start loop for trying lead %compensator pole. a=input('Enter a Test Lead Compensator Pole, (s+a). a = '); %Enter test lead compensator pole. numge=conv(numg,[1 b]); %Generate numerator of Gc(s)G(s). denge=conv([1 a],deng); %Generate denominator of Gc(s)G(s). Ge=tf(numge,denge); %Create Ge(s)=Gc(s)G(s). wn=4/(Ts*z); %Evaluate desired natural frequency. clf %Clear graph on screen. rlocus(Ge) %Plot compensated root locus with %test lead compensator pole. axis([-10,10,-10,10]) %Change lead-compensated root locus %axes. sgrid(z,wn) %Overlay grid on lead-compensated %root locus. title(['Lead-Compensated Root Locus with ' , num2str(pos),... '% Overshoot Line, Lead Pole at ', num2str(-a),... ' and Required Wn']) %Add title to lead-compensated root %locus. done=input('Are you done? (y=0,n=1) '); %Set loop flag. end %End loop for trying compensator %pole. [K,p]=rlocfind(Ge); %Generate gain, K, and closed-loop %poles, p, for point selected %interactively on the root locus. 'Gc(s)' %Display label. Gc=tf([1 b],[1 a]) %Display lead compensator. 'Gc(s)G(s)' %Display label. Ge %Display Gc(s)G(s). 'Closed-loop poles = ' %Display label. p %Display lead-compensated system's %closed-loop poles. f=input('Give pole number that is operating point '); %Choose lead-compensated system %dominant pole. 'Summary of estimated specifications for selected point on lead' 'compensated root locus' %Display label. operatingpoint=p(f) %Display lead-compensated dominant %pole. gain=K %Display lead-compensated gain. estimated_settling_time=4/abs(real(p(f))) %Display lead-compensated settling %time. estimated_peak_time=pi/abs(imag(p(f))) %Display lead-compensated peak time. estimated_percent_overshoot=pos %Display lead-compensated percent %overshoot. estimated_damping_ratio=z %Display lead-compensated damping %ratio. estimated_natural_frequency=sqrt(real(p(f))^2+imag(p(f))^2) %Display lead-compensated natural %frequency. s=tf([1 0],1); %Create transfer function, "s". sGe=s*Ge; %Create sGe(s) to evaluate Kv. sGe=minreal(sGe); %Cancel common poles and zeros. Kv=dcgain(K*sGe) %Display lead-compensated Kv. ess=1/Kv %Display lead-compensated steady%state error for unit ramp input.
Solutions to Problems 627
%Display label. %Create and display lead-compensated %T(s). 'Press any key to continue and obtain the lead-compensated step' 'response' %Display label pause step(T) %Plot step response for lead %compensated system. title(['Lead-Compensated System with ' ,num2str(pos),'% Overshoot']) %Add title to step response of PD %compensated system. pause 'Digital design' %Print label. T=0.01 %Define sampling interval. clf %Clear graph. 'Gc(s) in polynomial form' %Print label. Gcs=K*Gc %Create Gc(s) in polynomial form. 'Gc(s) in polynomial form' %Print label. Gcszpk=zpk(Gcs) %Create Gc(s) in factored form. 'Gc(z) in polynomial form via Tustin Transformation' %Print label. Gcz=c2d(Gcs,T,'tustin') %Form Gc(z) via Tustin transformation. 'Gc(z) in factored form via Tustin Transformation' %Print label. Gczzpk=zpk(Gcz) %Show Gc(z) in factored form. 'Gp(s) in polynomial form' %Print label. Gps=G %Create Gp(s) in polynomial form. 'Gp(s) in factored form' %Print label. Gpszpk=zpk(Gps) %Create Gp(s) in factored form. 'Gp(z) in polynomial form' %Print label. Gpz=c2d(Gps,T,'zoh') %Form Gp(z) via zoh transformation. 'Gp(z) in factored form' %Print label. Gpzzpk=zpk(Gpz) %Form Gp(z) in faactored form. pole(Gpz) %Find poles of Gp(z). Gez=Gcz*Gpz; %Form Ge(z) = Gc(z)Gp(z). 'Ge(z) = Gc(z)Gp(z) in factored form' %Print label. Gezzpk=zpk(Gez) %Form Ge(z) in factored form. 'z-1' %Print label. zm1=tf([1 -1],1,T) %Form z-1. zm1Gez=minreal(zm1*Gez,.00001); %Cancel common factors. '(z-1)Ge(z)' %Print label. zm1Gezzpk=zpk(zm1Gez) %Form & display (z-1)Ge(z) in %factored form. pole(zm1Gez) %Find poles of (z-1)Ge(z). Kv=10*dcgain(zm1Gez) %Find Kv. Tz=feedback(Gez,1) %Find closed-loop %transfer function, T(z) step(Tz) %Find step reponse. title('Closed-Loop Digital Step Response') %Add title to step response.
'T(s)' T=feedback(K*Ge,1)
Computer response:
ans = Design of digital lead compensation ans = Uncompensated System ans = G(s)
628 Chapter 13: Digital Control Systems
Transfer function: 1 ------------------s^3 + 13 s^2 + 40 s Type desired percent overshoot 10 Select a point in the graphics window selected_point = -1.6435 + 2.2437i ans = Closed-loop poles = p = -9.6740 -1.6630 + 2.2492i -1.6630 - 2.2492i Give pole number that is operating point ans = Summary of estimated specifications for selected point on ans = uncompensated root locus operatingpoint = -1.6630 + 2.2492i gain = 75.6925 estimated_settling_time = 2.4053 estimated_peak_time = 1.3968 estimated_percent_overshoot = 10 estimated_damping_ratio = 0.5912 estimated_natural_frequency = 2
Solutions to Problems 629
2.7972 Kv = 1.8923 ess = 0.5285 ans = T(s) Transfer function: 75.69 --------------------------s^3 + 13 s^2 + 40 s + 75.69 ans = Press any key to go to lead compensation Type Desired Settling Time 2 Type Lead Compensator Zero, (s+b). b= 5 Enter a Test Lead Compensator Pole, (s+a). a = Are you done? (y=0,n=1) 0 Select a point in the graphics window selected_point = -1.9709 + 2.6692i ans = Gc(s) Transfer function: s + 5 ------s + 6.8 ans = Gc(s)G(s) Transfer function: s + 5 ---------------------------------s^4 + 19.8 s^3 + 128.4 s^2 + 272 s ans = Closed-loop poles =
6.8
630 Chapter 13: Digital Control Systems
p = -10.7971 -5.0000 -2.0015 + 2.6785i -2.0015 - 2.6785i Give pole number that is operating point ans = Summary of estimated specifications for selected point on lead ans = compensated root locus operatingpoint = -2.0015 + 2.6785i gain = 120.7142 estimated_settling_time = 1.9985 estimated_peak_time = 1.1729 estimated_percent_overshoot = 10 estimated_damping_ratio = 0.5912 estimated_natural_frequency = 3.3437 Kv = 2.2190 ess = 0.4507 ans = T(s) 3
Solutions to Problems 631
Transfer function: 120.7 s + 603.6 -------------------------------------------s^4 + 19.8 s^3 + 128.4 s^2 + 392.7 s + 603.6 ans = Press any key to continue and obtain the lead-compensated step ans = response ans = Digital design T = 0.0100 ans = Gc(s) in polynomial form Transfer function: 120.7 s + 603.6 --------------s + 6.8 ans = Gc(s) in polynomial form Zero/pole/gain: 120.7142 (s+5) -------------(s+6.8) ans = Gc(z) in polynomial form via Tustin Transformation Transfer function: 119.7 z - 113.8 --------------z - 0.9342 Sampling time: 0.01 ans = Gc(z) in factored form via Tustin Transformation
632 Chapter 13: Digital Control Systems
Zero/pole/gain: 119.6635 (z-0.9512) ------------------(z-0.9342) Sampling time: 0.01 ans = Gp(s) in polynomial form Transfer function: 1 ------------------s^3 + 13 s^2 + 40 s ans = Gp(s) in factored form Zero/pole/gain: 1 ------------s (s+8) (s+5) ans = Gp(z) in polynomial form Transfer function: 1.614e-007 z^2 + 6.249e-007 z + 1.512e-007 -----------------------------------------z^3 - 2.874 z^2 + 2.752 z - 0.8781 Sampling time: 0.01 ans = Gp(z) in factored form Zero/pole/gain: 1.6136e-007 (z+3.613) (z+0.2593) -------------------------------(z-1) (z-0.9512) (z-0.9231) Sampling time: 0.01 ans = 1.0000 0.9512 0.9231 ans = Ge(z) = Gc(z)Gp(z) in factored form Zero/pole/gain
Solutions to Problems 633
1.9308e-005 (z+3.613) (z-0.9512) (z+0.2593) ------------------------------------------(z-1) (z-0.9512) (z-0.9342) (z-0.9231) Sampling time: 0.01 ans = z-1 Transfer function: z - 1 Sampling time: 0.01 ans = (z-1)Ge(z) Zero/pole/gain: 1.9308e-005 (z+3.613) (z+0.2593) -------------------------------(z-0.9342) (z-0.9231) Sampling time: 0.01 ans = 0.9342 0.9231 Kv = 0.2219 Transfer function: 1.931e-005 z^3 + 5.641e-005 z^2 - 5.303e-005 z - 1.721e-005 ----------------------------------------------------------z^4 - 3.809 z^3 + 5.438 z^2 - 3.45 z + 0.8203 Sampling time: 0.01
634 Chapter 13: Digital Control Systems
Solutions to Problems 635
+
636 Chapter 13: Digital Control Systems
SOLUTIONS TO DESIGN PROBLEMS
27. a. Push negative sign from vehicle dynamics to the left past the summing junction. The computer will be the area inside the large box with the inputs and outputs shown sampled. G(s) is the combined rudder actuator and vehicle dynamics. Also, the yaw rate sensor is shown equivalently before the integrator with unity feedback.
0.25(s+0.437) where G(s) = (s+2)(s+1.29)(s+0.193) .
Solutions to Design Problems 637
b. Add a phantom sampler at the output and push G(s) along with its sample and hold to the right past the pickoff point.
G(s) Move the outer-loop sampler to the output of s and write the z transforms of the transfer functions.
where 0.25(s+0.437) G1(s) = (1 - e-Ts) s(s+2)(s+1.29)(s+0.193) and 0.25(s+0.437) . G2(s) = (1 - e-Ts) 2 s (s+2)(s+1.29)(s+0.193) Now find the z transforms of G1(s) and G2(s). For G1(z).
For G2(z):
638 Chapter 13: Digital Control Systems
Now find the closed-loop transfer function. First find the equivalent forward transfer function.
Thus, Ge(z) T(z) = 1 + G (z) e Substituting values,
c. Using Ge(z), plot the root locus and see where it crosses the unit circle.
Solutions to Design Problems 639
The root locus crosses the unit circle when 3.8642x10-5K = 5.797x10-4, or K = 15. 28. a. First find G(z). z-1 G(z) = K z z{
2
s
1 } s + 7 s + 1220
2
7 s + 3.5 - 34.4 1207.8 - 7 1 + 1220 1 ) } z-1 = K z z{6.7186x10-7( 2 2 s s s + 3.5 + 1207.8
For T = 0.1,
z z z z + 0.66582 z z-1 + 7.8472 + 122 -7 = K z {6.7186x107(7 2 2 z -1 z -1 z + 1.3316 z + 0.49659 z + 1.3316 z + 0.49659
2
2
)}
G(z) = K 7.940510 -5
z + 0.63582 + 0.49355 i z + 0.63582 - 0.49355 i z - 1 z + 0.66582 + 0.2308 i z + 0.66582 - 0.2308 i
640 Chapter 13: Digital Control Systems
b.
c. The root locus intersects the unit circle at -1 with a gain, 7.9405x10-5K = 10866, or
0 < K < 136.84x106.
d. Program:
%Digitize G1(s) preceded by a sample and hold numg1=1; deng1=[1 7 1220 0]; 'G1(s)' G1s=tf(numg1,deng1) 'G(z)' Gz=c2d(G1s,0.1,'zoh') [numgz,dengz]=tfdata(Gz,'v'); 'Zeros of G(z)' roots(numgz) 'Poles of G(z)' roots(dengz) %Plot root locus rlocus(Gz) title(['Root Locus']) [K,p]=rlocfind(Gz)
Computer response:
ans = G1(s)
Solutions to Design Problems 641
Transfer function: 1 -------------------s^3 + 7 s^2 + 1220 s ans = G(z) Transfer function: 7.947e-005 z^2 + 0.0001008 z + 5.15e-005 ---------------------------------------z^3 + 0.3316 z^2 - 0.8351 z - 0.4966 Sampling time: 0.1 ans = Zeros of G(z) ans = -0.6345 + 0.4955i -0.6345 - 0.4955i ans = Poles of G(z) ans = 1.0000 -0.6658 + 0.2308i -0.6658 - 0.2308i Select a point in the graphics window selected_point = -0.9977 K = 1.0885e+004 p = -0.9977 -0.0995 + 0.2330i -0.0995 - 0.2330i
See part (b) for root locus plot.
29.
z-1 20000 z-1 a. First find G(z). G(z) = K z z{ 2 } = K z z{2 1 - 2 1 + 200 12 } s (s+100) s s + 100 s z z -1 -2 z +2 z +2 For T = 0.01, G(z) = K 2 z -1 z - 0.36788 z z -1
642 Chapter 13: Digital Control Systems
z + 0.71828 . z - 1 z - 0.36788
= 0.73576 K
b. Plotting the root locus. Finding the intersection with the unit circle yields 0.73576K = 1.178. Thus, 0
< K < 1.601 for stability.
c. Using the root locus, we find the intersection with the 15% overshoot curve ( = 0.517) at
0.5955 + j0.3747 with 0.73576K = 0.24. Thus K = 0.326.
d. Program:
%Digitize G1(s) preceded by a sample and hold numg1=20000; deng1=[1 100 0]; 'G1(s)' G1s=tf(numg1,deng1) 'G(z)' Gz=c2d(G1s,0.01,'zoh') [numgz,dengz]=tfdata(Gz,'v'); 'Zeros of G(z)' roots(numgz) 'Poles of G(z)' roots(dengz) %Input transient response specifications Po=input('Type %OS '); %Determine damping ratio z=(-log(Po/100))/(sqrt(pi^2+log(Po/100)^2)) %Plot root locus rlocus(Gz) zgrid(z,0) title(['Root Locus']) [K,p]=rlocfind(Gz) %Allows input by selecting point on graphic.
Computer response:
ans = G1(s)
Solutions to Design Problems 643
Transfer function: 20000 ----------s^2 + 100 s ans = G(z) Transfer function: 0.7358 z + 0.5285 ---------------------z^2 - 1.368 z + 0.3679 Sampling time: 0.01 ans = Zeros of G(z) ans = -0.7183 ans = Poles of G(z) ans = 1.0000 0.3679 Type %OS z = 0.5169 Select a point in the graphics window selected_point = 0.5949 + 0.3888i K = 0.2509 p = 0.5917 + 0.3878i 0.5917 - 0.3878 15
644 Chapter 13: Digital Control Systems
30.
G PID (s) =
Substituting Eq. (13.88) with T = 1/3 second,
0.5857(s + 0.19)(s + 0.01) . s
Gc (z) =
3.632z2 - 7.028z + 3.397 3.632(z - 0.9967)(z - 0.9386) = 2 (z + 1)(z - 1) z -1
Drawing the flow diagram yields
Solutions to Design Problems 645
3.632
-7.028
0
3.397 T = 1/3 second
31.
-1
a. From Chapter 9, the plant without the pots and unity gain power amplifier is
64.88 (s+53.85) Gp(s) = (s2 + 8.119s + 376.3) (s^2 + 15.47s + 9283)
The PID controller and notch filter with gain adjusted for replacement of pots (i.e. divided by 100) was 26.82 (s+24.1) (s+0.1) (s^2 + 16.s + 9200) --------------------------------------------------------------s (s+60)^2
Gc(s) =
Thus, Ge(s) = Gp(s)Gc(s) is 1740.0816 (s+53.85)(s2 + 16s + 9200)(s+24.09)(s+0.1) Get(s) = s (s2 + 8.119s + 376.3) (s^2 + 15.47s + 9283)(s+60)2
A Bode magnitude plot of Ge(s) shows c = 36.375 rad/s. Thus, the maximum T should be in the range 0.15/c to 0.5/c or 4.1237e-03 to 1.3746e-02. Let us select T = 0.001. Performing a Tustin transformation on Gc(s) yields 5.166e04 z4 - 2.041e05 z3 + 3.029e05 z2 - 2.001e05 z + 4.963e04 Gc(z) = -----------------------------------------------------------------------------------------Z4 - 1.883 z3 - 0.1131 z2 + 1.883 z - 0.8869
b. Drawing the flowchart
646 Chapter 13: Digital Control Systems
5.166e4
-2.041e5
+ ++
-1.883 - --
3.029e5
-0.1131
-2.001e5
1.883
4.963e4
-0.8869
T = 0.001
c. Program:
syms s 'Compensator from Chapter 9' T=.001 Gc=26.82*(s^2+16*s+9200)*(s+24.09)*(s+.1)/(s*((s+60)^2)); Gc=vpa(Gc,4); [numgc,dengc]=numden(Gc); numgc=sym2poly(numgc); dengc=sym2poly(dengc); Gc=tf(numgc,dengc); 'Gc(s)' Gczpk=zpk(Gc) 'Gc(z)' Gcz=c2d(Gc,T,'tustin') 'Gc(z)' Gczzpk=zpk(Gcz) 'Plant from Chapter 9'
Solutions to Design Problems 647
Gp=64.88*(s+53.85)/[(s^2+15.47*s+9283)*(s^2+8.119*s+376.3)]; Gp=vpa(Gp,4); [numgp,dengp]=numden(Gp); numgp=sym2poly(numgp); dengp=sym2poly(dengp); 'Gp(s)' Gp=tf(numgp,dengp) 'Gp(s)' Gpzpk=zpk(Gp) 'Gp(z)' Gpz=c2d(Gp,T,'zoh') 'Gez=Gcz*Gpz' Gez=Gcz*Gpz Tz=feedback(Gez,1); t=0:T:1; step(Tz,t) pause t=0:T:50; step(Tz,t)
Computer response:
ans = Compensator from Chapter 9 T = 0.0010 ans = Gc(s) Zero/pole/gain: 26.82 (s+24.09) (s+0.1) (s^2 + 16s + 9198) ------------------------------------------s (s+60)^2 ans = Gc(z) Transfer function: 5.17e004 z^4 - 2.043e005 z^3 + 3.031e005 z^2 - 2.002e005 z + 4.966e004 ---------------------------------------------------------------------z^4 - 1.883 z^3 - 0.1131 z^2 + 1.883 z - 0.8869 Sampling time: 0.001 ans = Gc(z) Zero/pole/gain: 51699.4442 (z-1) (z-0.9762) (z^2 - 1.975z + 0.9842) ---------------------------------------------------(z+1) (z-1) (z-0.9417)^2 Sampling time: 0.001 ans =
648 Chapter 13: Digital Control Systems
Plant from Chapter 9 ans = Gp(s) Transfer function: 64.88 s + 3494 ---------------------------------------------------s^4 + 23.59 s^3 + 9785 s^2 + 8.119e004 s + 3.493e006 ans = Gp(s) Zero/pole/gain: 64.88 (s+53.85) ---------------------------------------------(s^2 + 8.119s + 376.3) (s^2 + 15.47s + 9283) ans = Gp(z) Transfer function: 1.089e-008 z^3 + 3.355e-008 z^2 - 3.051e-008 z - 1.048e-008 ----------------------------------------------------------z^4 - 3.967 z^3 + 5.911 z^2 - 3.92 z + 0.9767 Sampling time: 0.001 ans = Gez=Gcz*Gpz Transfer function: 0.000563 z^7 - 0.0004901 z^6 - 0.005129 z^5 + 0.01368 z^4 - 0.01328 z^3 + 0.004599 z^2 + 0.0005822 z - 0.0005203 -------------------------------------------------------------------------z^8 - 5.85 z^7 + 13.27 z^6 - 12.72 z^5 - 0.6664 z^4 + 13.25 z^3 - 12.74 z^2 + 5.317 z - 0.8662 Sampling time: 0.001
Solutions to Design Problems 649