Basic Root Locus

Topic #2 16.30/31 Feedback Control Systems Basic Root Locus • Basic aircraft control concepts • Basic control approaches Fall 2010 16.30/31 2–2 Aircraft Longitudinal Control • Consider the short period approximate model of an 747 aircraft. ẋsp = Aspxsp + Bspδe where δe is the elevator input, and � � � � U0 w Zw /m xsp = , Asp = −1 −1 (Mw + Mẇ Zw /m) Iyy (Mq + Mẇ U0) Iyy q � � Zδe /m Bsp = −1 Iyy (Mδe + Mw˙ Zδe /m) • Add that θ̇ = q, so sθ = q • Take the output as θ, input is δe, then form the transfer function1 � 1 q(s) 1� θ(s) = = 0 1 (sI − Asp)−1Bsp δe(s) s δe(s) s • For the 747 (40Kft, M = 0.8) this reduces to: θ(s) 1.1569s + 0.3435 =− 2 ≡ Gθδe (s) δe(s) s(s + 0.7410s + 0.9272) so that the dominant roots have a frequency of approximately 1 rad/sec and damping of about 0.4 Pole−Zero Map 1 1 0.68 0.54 0.42 0.3 0.2 0.09 0.8 0.8 0.6 0.84 0.6 0.4 0.4 Imaginary Axis 0.95 0.2 0.2 0 −0.2 0.2 0.95 −0.4 −0.6 0.4 0.6 0.84 −0.8 −1 −1 0.8 0.68 −0.9 0.54 −0.8 −0.7 −0.6 0.42 −0.5 −0.4 Real Axis 0.3 0.2 −0.3 −0.2 0.09 −0.1 10 0.1 Fig. 1: Note - this is the Pole-zero map for Gqδe 1 Much more on how to do this part later Fall 2010 16.30/31 2–3 • Basic problem is that there are vast quantities of empirical data to show that pilots do not like the flying qualities of an aircraft with this combination of frequency and damping • What is preferred? Undamped natural frequency ωs rad/sec 7 6 5 POOR ACCEPTABLE 4 3 SATISFACTORY 2 UNACCEPTABLE 1 0 0.1 0.2 0.4 0.6 0.8 1 2 4 Damping ratio ςs Image by MIT OpenCourseWare. Fig. 2: “Thumb Print” criterion • This criterion was developed in 1950’s, and more recent data is pro­ vided in MILSPEC8785C • Based on this plot, a good target: frequency ≈ 3 rad/sec and damping of about ≈ 0.6 • Problem is that the short period dynamics are no where near these numbers, so we must modify them. • Could do it by redesigning the aircraft, but it is a bit late for that. . . Fall 2010 16.30/31 2–4 First Short Period Autopilot • First attempt to control the vehicle response: measure θ and feed it back to the elevator command δe. • Unfortunately the actuator is slow, so there is an apparent lag in the response that we must model δec 4 s+4 δea Gθδe (s) θ θc kθ − • Dynamics: δea is the actual elevator deflection, δec is the actuator command created by our controller 4 θ = Gθδe (s)δea; δea = H(s)δec; H(s) = s+4 The control is just basic proportional feedback δec = −kθ (θ − θc) which gives that θ = −Gθδe (s)H(s)kθ (θ − θc) or that θ(s) Gθδe (s)H(s)kθ = θc(s) 1 + Gθδe (s)H(s)kθ • Looks good, but how do we analyze what is going on? • Need to be able to predict where the poles are going as a function of kθ ⇒ Root Locus Fall 2010 16.30/31 2–5 Root Locus Basics r e − Gc(s) u Gp(s) y • Assume that the plant transfer function is of the form �npz Np (s − zpi) Gp = Kp = Kp �ni pp Dp i (s − ppi ) and the controller transfer function is �ncz Nc (s − zci) Gc(s) = Kc = Kc �nicp Dc i (s − pci ) 2 • Assume that npp > npz and ncp > ncz • Signals are: u y r e control commands output/measurements reference input response error • Unity feedback form. We could add the controller Gc in the feedback path without changing the pole locations. • Will discuss performance and add disturbances later, but for now just focus on the pole locations 2 Errata: Added n values for the number of poles and zeros Fall 2010 16.30/31 2–6 • Basic questions: • Analysis: Given Nc and Dc, where do the closed loop poles go as a function of Kc? • Synthesis: Given Kp, Np and Dp, how should we chose Kc, Nc, Dc to put the closed loop poles in the desired locations? • Block diagram analysis: Since y = GpGce and e = r − y, then easy to show that y Gc Gp = ≡ Gcl (s) r 1 + Gc Gp where Kc Kp N c N p DcDp + KcKpNcNp is the closed loop transfer function Gcl (s) = • Denominator called characteristic equation φc(s) and the roots of φc(s) = 0 are called the closed-loop poles (CLP). • The CLP are clearly functions of Kc for a given Kp, Np, Dp, Nc, Dc ⇒ a “locus of roots” [Evans, 1948] Fall 2010 16.30/31 2–7 Root Locus Analysis • General root locus is hard to determine by hand and requires Matlab tools such as rlocus(num,den) to obtain full result, but we can get some important insights by developing a short set of plotting rules. • Full rules in FPE, page 279 (4th edition). • Basic questions: 1. What points are on the root locus? 2. Where does the root locus start? 3. Where does the root locus end? 4. When/where is the locus on the real line? 5. Given that s0 is found to be on the locus, what gain is need for that to become the closed-loop pole location? 6. What are the departure and arrival angles? 7. Where are the multiple points on the locus? • Question #1: is point s0 on the root locus? Assume that Nc and Dc are known, let Nc Np Ld = and K = KcKp Dc Dp ⇒ φc(s) = 1 + KLd(s) = 0 So values of s for which Ld(s) = −1/K, with K real are on the RL. • For K positive, s0 is on the root locus if �Ld(s0) = 180◦ ± l · 360◦, l = 0, 1, . . . • If K negative, s0 is on the root locus if [0◦ locus] �Ld(s0) = 0◦ ± l · 360◦, l = 0, 1, . . . These are known as the phase conditions. Fall 2010 16.30/31 2–8 • Question #2: Where does the root locus start? φc = 1 + K Nc N p =0 DcDp ⇒ DcDp + KNcNp = 0 So if K → 0, then locus starts at solutions of DcDp = 0 which are the poles of the plant and compensator. • Question #3: Where does the root locus end? Already shown that for s0 to be on the locus, must have 1 Ld(s0) = − K So if K → ∞, the poles must satisfy: Ld = N c Np =0 DcDp • There are several possibilities: 1. Poles are located at values of s for which NcNp = 0, which are the zeros of the plant and the compensator 2. If Loop Ld(s) has more poles than zeros • As |s| → ∞, |Ld(s)| → 0, but we must ensure that the phase condition is still satisfied. Fall 2010 16.30/31 2–9 • More details as K → ∞: • Assume there are n zeros and p poles of Ld(s) • Then for large |s|, Ld(s) ≈ 1 (s − α)p−n • So the root locus degenerates to: 1 1+ =0 (s − α)p−n • So n poles head to the zeros of Ld(s) • Remaining p − n poles head to |s| = ∞ along asymptotes defined by the radial lines 180◦ + 360◦ · (l − 1) φl = p−n l = 1, 2, . . . so that the number of asymptotes is governed by the number of poles compared to the number of zeros (relative degree). • If zi are the zeros if Ld and pj are the poles, then the centroid of the asymptotes is given by: α= p � pj − n � p−n zi Fall 2010 16.30/31 2–10 • Example: L(s) = s−4 Im ✖ Re • Number of asymptotes and α? • Example G(s) = s+1 s2(s + 4) Im ✖ Re ✖ • Number of asymptotes and α? • Example G(s) = s−1 s2(s − 4) Im ✖ • Number of asymptotes and α? ✖ Re Fall 2010 16.30/31 2–11 • Question #4: When/where is the locus on the real line? • Locus points on the real line are to the left of an odd number of real axis poles and zeros [K positive]. • Follows from the phase condition and the fact that the phase contribution of the complex poles/zeros cancels out • Question #5: Given that s0 is found to be on the locus, what gain is needed for that to become the closed-loop pole location? • Need � � � Dp(s0)Dc(s0) � 1 � K≡ = �� |Ld(s0)| Np(s0)Nc(s0) � • Since K = KpKc, sign of Kc depends on sign of Kp ∗ e.g., assume that �Ld(s0) = 180◦, then need Kc and Kp to be same sign so that K > 0 Fall 2010 16.30/31 2–12 Root Locus Examples Im Re ✖ Fig. 3: Basic Im ✖ Re ✖ Fig. 4: Two poles Im ✖ Re ✖ Fig. 5: Add zero • Examples similar to control design process: add compensator dynam­ ics to modify root locus and then chose gain to place CLP at desired location on the locus. Fall 2010 16.30/31 2–13 Im ✖ Re ✖ ✖ Fig. 6: Three poles Im ✖ Re ✖ ✖ Fig. 7: Add a zero again Im ✖ ✖ Re Fig. 8: Complex Case Im ✖ ✖ ✖ Fig. 9: Very Complex Case Re Fall 2010 16.30/31 2–14 Performance Issues • Interested in knowing how well our closed loop system can track var­ ious inputs • Steps, ramps, parabolas • Both transient and steady state • For perfect steady state tracking want error to approach zero lim e(t) = 0 t→∞ • Can determine this using the closed-loop transfer function and the final value theorem lim e(t) = lim se(s) t→∞ s→0 • So for a step input r(t) = 1(t) → r(s) = 1/s y(s) Gc(s)Gp(s) = r(s) 1 + Gc(s)Gp(s) e(s) 1 = r(s) 1 + Gc(s)Gp(s) y(s) = Gc(s)Gp(s) e(s) so in the case of a step input, we have r(s) 1/s = 1 + Gc(s)Gp(s) 1 + Gc(s)Gp(s) 1/s 1 ⇒ lim se(s) = lim s = ≡ e(∞) s→0 s→0 1 + Gc (s)Gp (s) 1 + Gc(0)Gp(0) e(s) = • So the steady state error to a step is given by 1 ess = 1 + Gc(0)Gp(0) • To make the error small, we need to make one (or both) of Gc(0), Gp(0) very large Fall 2010 16.30/31 2–15 • Clearly if Gp(s) has a free integrator (or two) so that it resembles 1 sn (s+α)m with n ≥ 1, then lim Gp(s) → ∞ s→0 ⇒ ess → 0 • Can continue this discussion by looking at various input types (step, ramp, parabola) with systems that have a different number of free integrators (type), but the summary is this: step ramp parabola type 0 1 1 + Kp ∞ ∞ type 1 0 1 Kv ∞ type 2 0 0 1 Ka where Kp = lim Gc(s)Gp(s) Position Error Constant Kv = lim sGc(s)Gp(s) Velocity Error Constant s→0 s→0 Ka = lim s2Gc(s)Gp(s) s→0 Acceleration Error Constant which are a good simple way to keep track of how well your system is doing in terms of steady state tracking performance. Fall 2010 16.30/31 2–16 Dynamic Compensation • For a given plant, can draw a root locus versus K. But if desired pole locations are not on that locus, then need to modify it using dynamic compensation. • Basic root locus plots give us an indication of the effect of adding compensator dynamics. But need to know what to add to place the poles where we want them. • New questions: • What type of compensation is required? • How do we determine where to put the additional dynamics? • There are three classic types of controllers u = Gc(s)e 1. Proportional feedback: Gc ≡ Kg a gain, so that Nc = Dc = 1 • Same case we have been looking at. 2. Integral feedback: u(t) = Ki � t 0 e(τ )dτ ⇒ Gc(s) = • Used to reduce/eliminate steady-state error Ki s • If e(τ ) is approximately constant, then u(t) will grow to be very large and thus hopefully correct the error. Fall 2010 16.30/31 2–17 • Consider error response of Gp(s) = 1/(s + a)(s + b) (a > 0, b > 0) to a step, r(t) = 1(t) → r(s) = 1/s where e 1 = = S(s) r 1 + Gc G p → e(s) = r(s) (1 + GcGp) • where S(s) is the Sensitivity Transfer Function for the closed-loop system • To analyze error, use FVT limt→∞ e(t) = lims→0 se(s) so that with proportional control, �s� 1 1 lim ess = lim = K t→∞ s→0 s 1 + Kg Gp (s) 1 + abg so can make ess small, but only with a very large Kg • With integral control, lims→0 Gc(s) = ∞, so ess → 0 • Integral control improves the steady state, but this is at the expense of the transient response ∗ Typically gets worse because the system is less well damped Fall 2010 16.30/31 2–18 1 , add integral feedback to (s + a)(s + b) improve the steady state response. Example #1: G(s) = Im ✖✖ ✖ Re Fig. 10: RL after adding integral FB • Increasing Ki to increase speed of the response pushes the poles towards the imaginary axis → more oscillatory response. Now combine proportional and integral (PI) feedback: K2 K1 s + K 2 = s s which introduces a pole at the origin and zero at s = −K2/K1 Gc = K1 + • PI solves many of the problems with just integral control Im ✖✖ ✖ Re Fig. 11: RL with proportional and integral FB Fall 2010 16.30/31 2–19 3. Derivative Feedback: u = Kdė so that Gc(s) = Kds • Does not help with the steady state • Provides feedback on the rate of change of e(t) so that the control can anticipate future errors. Example # 2: G(s) = with Gc(s) = Kds 1 , (a > 0, b > 0) (s − a)(s − b) Im ✖ ✖ Re Fig. 12: RL with derivative FB • Derivative feedback is very useful for pulling the root locus into the LHP - increases damping and more stable response. Typically used in combination with proportional feedback to form proportional-derivative feedback PD Gc(s) = K1 + K2s which moves the zero from the origin. • Unfortunately pure PD is not realizable in the lab as pure differenti­ ation of a measured signal is typically a bad idea • Typically use band-limited differentiation instead, by rolling-off the PD control with a high-frequency pole (or two). Fall 2010 16.30/31 2–20 Controller Synthesis • First determine where the poles should be located • Will proportional feedback do the job? • What types of dynamics need to be added? Use main building block GB (s) = Kc (s + z) (s + p) • Looks like various controllers, depending how Kc, p, and z picked • If pick z > p, with p small, then Im GB (s) ≈ Kc ✖ (s + z) s Re which is essentially a PI compensator, called a lag. • If pick p � z, then at low frequency, the impact of p/(s + p) is small, so Im ✖ Re GB (s) ≈ Kc(s + z) which is essentially PD compensator, called a lead. • Various algorithms exist to design the components of the lead and lag compensators Fall 2010 16.30/31 2–21 Classic Root Locus Approach • Consider a simple system Gp = s−2 for which we want the closed loop poles to be at −1 ± 2j • Will proportional control be sufficient? no • So use compensator with 1 pole. Gc = K (s + z) (s + p) So there are 3 CLP. • To determine how to pick the p, z, and k, we must use the phase and magnitude conditions of the RL • To proceed, evaluate the phase of the loop s+z Ld(s) = (s + p)s2 at s0 = −1 + 2j. Since we want s0 to be on the new locus, we know that �Ld(s0) = 180◦ ± 360◦l Im (-1,2j) • ✖ −p α β −z ✖ φ0 Re Fig. 13: Phase Condition • As shown in the figure, there are four terms in �Ld(s0) – the two poles at the origin contribute 117◦ each • Given the assumed location of the compensator pole/zero, can work out their contribution as well Fall 2010 16.30/31 2–22 • Geometry for the real zero: tan α = 2 p−1 2 z−1 and for the real pole: tan β = • Since we expect the zero to be closer to the origin, put it first on the negative real line, and then assume that p = γz, where typically 5 ≤ γ ≤ 10 is a good ratio. • So the phase condition gives: −2(117◦) + α − β = 180◦ � � � � 2 2 − arctan = 53◦ arctan z−1 10z − 1 but recall that tan(A − B) = so tan(A) − tan(B) 1 + tan(A) tan(B) ( z−2 1 ) − ( 10z2−1 ) 1 + ( z−2 1 )( 10z2−1 ) = 1.33 which give z = 2.2253, p = 22.2531, kc = 45.5062 1 2 3 4 5 6 7 8 9 10 11 12 13 % RL design using angles clear all target = −1+2*j; phi origin= 180−atan(imag(target)/−real(target))*180/pi; syms z M; ratio=10; phi z=(imag(target)/(z+real(target))); phi p=(imag(target)/(ratio*z+real(target))); M=(phi z−phi p)/(1+phi z * phi p); test=solve(M−tan(pi/180*(2*phi origin−180))); Z=eval(test(1)); P=ratio*Z; K=1/abs((target+Z)/(targetˆ2*(target+P))); [Z P K] Fall 2010 16.30/31 2–23 Pole Placement • Another option for simple systems is called pole placement. • Know that the desired characteristic equation is φd(s) = (s2 + 2s + 5)(s + α) = 0 • Actual closed loop poles solve: φc(s) = 1 + GpGc = 0 → s2(s + p) + K(s + z) = 0 → s3 + s2p + Ks + Kz = 0 • Clearly need to pull the poles at the origin into the LHP, so need a lead compensator → Rule of thumb:3 take p = (5 − 10)z. • Compare the characteristic equations: φc(s) = s2 + 10zs2 + Ks + Kz = 0 φd(s) = (s2 + 2s + 5)(s + α) = s3 + s2(α + 2) + s(2α + 5) + 5α = 0 gives s2 α + 2=10z s 2α + 5=K s0 5α=zK solve for α, z, K 25 5z ; α= 5 − 2z 5 − 2z → z = 2.23, α = 20.25, K = 45.5 K= 3 Errata: changed the rule of the pole zero ratio. Fall 2010 16.30/31 2–24 Fig. 14: CLP with pole placement Code: Pole Placement 1 2 3 4 5 6 7 8 9 10 11 % % Fall 2009 % close all figure(1);clf set(gcf,'DefaultLineLineWidth',2) set(gcf,'DefaultlineMarkerSize',10) set(gcf,'DefaultlineMarkerFace','b') clear all;%close all; set(0, 'DefaultAxesFontSize', 14, 'DefaultAxesFontWeight','demi') set(0, 'DefaultTextFontSize', 14, 'DefaultTextFontWeight','demi') 12 13 14 15 16 17 18 19 20 21 22 %Example: G(s)=1/2ˆ2 %Design Gc(s) to put the clp poles at −1 + 2j z=roots([−20 49 −10]);z=max(z),k=25/(5−2*z),alpha=5*z/(5−2*z), num=1;den=[1 0 0]; knum=k*[1 z];kden=[1 10*z]; rlocus(conv(num,knum),conv(den,kden)); hold;plot(−alpha+eps*j,'d');plot([−1+2*j,−1−2*j],'d');hold off r=rlocus(conv(num,knum),conv(den,kden),1)' axis([−25 5 −15 15]) print −dpng −r300 rl pp.png Fall 2010 16.30/31 2–25 Observations • In a root locus design it is easy to see the pole locations, and thus we can relatively easily identify the dominant time response • Caveat is that near pole/zero cancelation complicates the process of determining which set of poles will dominate • Some of the performance specifications are given in the frequency response, and it is difficult to determine those (and the corresponding system error gains) in the RL plot • Easy for low-order systems, very difficult / time consuming for higher order ones • As we will see, extremely difficult to identify the robustness margins using a RL plot • A good approach for a fast/rough initial design • Matlab tool called sisotool provides a great interface for designing and analyzing controllers