16.333 Lecture # 10 State Space Control ? Basic state space control approaches ? ? ? ? ? ? Fall 2004 16.333 9–1 State Space Basics ? State space models are of the form x˙(t) = Ax(t) + Bu(t) y(t) = Cx(t) + Du(t) with associated transfer function G(s) = C(sI ? A) ?1 B + D Note: must form symbolic inverse of matrix (sI ? A), which is hard. ? Time response: Homogeneous part x˙ = Ax, x(0) known – Take Laplace transform X(s) = (sI ? A) ?1 x(0) ? x(t) = L ?1 (sI ? A) ?1 x(0) I A A 2 – But can show (sI ? A) ?1 = + s 2 + s 3 + . . . s 1 so L ?1 (sI ? A) ?1 = I + At + 2! (At) 2 + . . . = e At – Gives x(t) = e At x(0) where e At is Matrix Exponential 1 3 Calculate in MATLAB R? using expm.m and not exp.m ? Time response: Forced Solution – Matrix case x˙ = Ax + Bu where x is an n-vector and u is a m-vector. Cam show t x(t) = e At x(0) + e A(t?τ ) Bu(τ )dτ 0 t y(t) = Ce At x(0) + Ce A(t?τ ) Bu(τ )dτ + Du(t) 0 – Ce At x(0) is the initial response – Ce A(t) B is the impulse response of the system. 1 MATLAB R? is a trademark of the Mathworks Inc. ? ? ? Fall 2004 16.333 9–2 Dynamic Interpretation ? Since A = T ΛT ?1 , then ? ?? ?? ? T e At | | λ 1 t . . . ? w . . 1 ? ? e = T e Λt T ?1 = ? v 1 ?? ?? .··· v n | | λ n t T e ? w n ? where we have written ? ? T ? w . 1 ? ?. T ?1 = ? . T ? w n ? which is a column of rows. ? Multiply this expression out and we get that n At λ i t T e = e v i w i i=1 ? Assume A diagonalizable, then x˙ = Ax, x(0) given, has solution x(t) = e At x(0) = T e Λt T ?1 x(0) n = e λ i t v i {w i T x(0)} i=1 n λ i t = e v i β i i=1 ? State solution is a linear combination of the system modes v i e λ i e λ i t – Determines the nature of the time response v i – Determines extent to which each state contributes to that mode β i – Determines extent to which the initial condition excites the mode ? ? ? ? Fall 2004 16.333 9–3 ? Note that the v i give the relative sizing of the response of each part of the state vector to the response. 1 v 1 (t) = e ?t mode 1 0 0.5 v 2 (t) = e ?3t mode 2 0.5 ? Clearly e λ i t gives the time modulation – λ i real – growing/decaying exponential response – λ i complex – growing/decaying exponential damped sinusoidal ? Bottom line: The locations of the eigenvalues determine the pole locations for the system, thus: – They determine the stability and/or performance & transient be- havior of the system. – It is their locations that we will want to modify with the controllers. Fall 2004 16.333 9–4 Full-state Feedback Controller ? Assume that the single-input system dynamics are given by x˙ = Ax + Bu y = Cx so that D = 0. – The multi-actuator case is quite a bit more complicated as we would have many extra degrees of freedom. ? Recall that the system poles are given by the eigenvalues of A. – Want to use the input u(t) to modify the eigenvalues of A to change the system dynamics. r u y ? OO // A, B, C // x  K Assume a full-state feedback of the form:? u = r ? Kx where r is some reference input and the gain K is R 1×n – If r = 0, we call this controller a regulator ? ? ? ? ? ? ? ? Fall 2004 16.333 9–5 ? Find the closed-loop dynamics: x˙ = Ax + B(r ? Kx) = (A ? BK)x + Br = A cl x + Br y = Cx ? Objective: Pick K so that A cl has the desired properties, e.g., – A unstable, want A cl stable – Put 2 poles at ?2± 2j ? Note that there are n parameters in K and n eigenvalues in A, so it looks promising, but what can we achieve? ? Example #1: Consider: 1 1 x˙ = 1 2 1 x + u 0 – Then det(sI ? A) = (s ? 1)(s ? 2)? 1 = s 2 s + 1 = 0 ? 3 so the system is unstable. – De?ne u = k 1 k 2 x = ?Kx, then ? ? ? ? ? 1 1 1 ? ? 1? k 1 1? k 2 A cl = A?BK = k 1 k 2 = 1 2 ? 0 1 2 – So then we have that det(sI ? A cl ) = s 2 + (k 1 ? 3)s + (1? 2k 1 + k 2 ) = 0 ? ? ? ? ? ? ? ? ? ? ? ? ? Fall 2004 16.333 9–6 – Thus, by choosing k 1 and k 2 , we can put λ i (A cl ) anywhere in the complex plane (assuming complex conjugate pairs of poles). ? To put the poles at s = ?5, ? 6, compare the desired characteristic equation (s + 5)(s + 6) = s 2 + 11s + 30 = 0 with the closed-loop one 2 s + (k 1 ? 3)x + (1? 2k 1 + k 2 ) = 0 to conclude that k 1 ? 3 = 11 k 1 = 14 1? 2k 1 + k 2 = 30 k 2 = 57 so that K = 14 57 , which is called Pole Placement. ? Of course, it is not always this easy, as the issue of controllability must be addressed. ? Example #2: Consider this system: 1 1 x˙ = 0 2 1 x + u 0 with the same control approach 1 1 1 ? ? 1? k 1 1? k 2 A cl = A ? BK = k 1 k 2 = 0 2 ? 0 0 2 so that det(sI ? A cl ) = (s ? 1 + k 1 )(s ? 2) = 0 The feedback control can modify the pole at s = 1, but it cannot move the pole at s = 2. ? This system cannot be stabilized with full-state feedback control. ? What is the reason for this problem? ? ? ? ? ? Fall 2004 16.333 9–7 – It is associated with loss of controllability of the e 2t mode. Basic test for controllability: rank M c = n? ? ? ? ? ? B AB 1 1 = 1 0 1 0 M c = 0 2 So that rank M c 1 = < 2. ? Must assume that the pair (A, B) are controllable. ? ? ? ? ? ? ? ? ?? ?? Fall 2004 16.333 9–8 Ackermann’s Formula ? The previous outlined a design procedure and showed how to do it by hand for second-order systems. – Extends to higher order (controllable) systems, but tedious. ? Ackermann’s Formula gives us a method of doing this entire design process is one easy step. K = 0 ... 0 1 M ?1 Φ d (A) c – M c = B AB ... A n?1 B – Φ d (s) is the characteristic equation for the closed-loop poles, which we then evaluate for s = A. – It is explicit that the system must be controllable because we are inverting the controllability matrix. ? Revisit Example #1: Φ d (s) = s 2 + 11s + 30 ? ? ? ? ? ? ? 1 1 1 1 1 ? 1 = B AB = =M c 0 1 2 0 0 1 So ? ? ? 1 1 ? ?1 ? ? ? 2 ? ? ? 1 1 1 1 K = 0 1 ? + 11 + 30I ? 0 1 1 2 1 2 ? ? 43 14 ? ? = 0 1 = 14 57 14 57 ? Automated in Matlab: place.m & acker.m (see polyvalm.m too) ? ? ? ? Fall 2004 16.333 9–9 ? Origins? For simplicity, consider a third-order system (case #2), but this extends to any order. ? ? ? ? ?a 1 ?a 2 ?a 3 1 A = ? 1 0 0 ? B = ? 0 ? C = b 1 b 2 b 3 0 1 0 0 – This form is useful because the characteristic equation for the 2 system is obvious ? det(sI ? A) = s 3 + a 1 s + a 2 s + a 3 = 0 Can show that ? ? ? ? ? ?a 1 ?a 2 ?a 3 1 A cl = A ? BK = ? 1 0 0 ? ? ? 0 ? k 1 k 2 k 3 0 1 0 0 ? ? ?a 1 ? k 1 ?a 2 ? k 2 ?a 3 ? k 3 = ? 1 0 0 ? 0 1 0 so that the characteristic equation for the system is still obvious: 2 Φ cl (s) = det(sI ? A cl ) = s 3 + (a 1 + k 1 )s + (a 2 + k 2 )s + (a 3 + k 3 ) = 0 ? We then compare this with the desired characteristic equation devel- oped from the desired closed-loop pole locations: 2 Φ d (s) = s 3 + (α 1 )s + (α 2 )s + (α 3 ) = 0 to get that ? a 1 + k 1 = α 1 ? k 1 = α . 1 ? a 1 . . . . . ? a n + k n = α n k n = α n ? a n ? Pole placement is a very powerful tool and we will be using it for most of our state space work. ? ? ? ? ? ? Fall 2004 16.333 9–10 Aircraft State Space Control ? Can now design a full state feedback controller for the dynamics: x˙ sp = A sp x sp + B sp δ e with desired poles being at ω n = 3 and ζ = 0.6 ? s = ?1.8 ± 2.4i φ d (s) = s 2 + 3.6s + 9 Ksp=place(Asp,Bsp,[roots([1 2*0.6*3 3^2])’]) w ? Design controller u = ?0.0264 ?2.3463 q ? With full model, could arrange it so phugoid poles remain in the same place, just move the ones associated with the short period mode s = ?1.8 ± 2.4i, ?0.0033 ± 0.0672i ev=eig(A); % damp short period, but leave the phugoid where it is Plist=[roots([1 2*.6*3 3^2])’ ev([3 4],1)’]; K1=place(A,B(:,1),Plist) ? ? ? ? ? ? u w ? ? ? ? ? u = 0.0026 ?0.0265 ?2.3428 0.0363 q θ ? ? ? ?? ? ? ? ? ? Fall 2004 16.333 9–11 ? Can also add the lag dynamics to short period model with θ included x˙ sp = A ? sp x sp + ? δ a δ a = 4 δ c B sp e ; e s + 4 e → x˙ δ = ?4x δ + 4δ e c , δ a = x δ e x˙ sp = ? A ? sp B sp x sp 0 + δ c e ? x˙ δ 0 ?4 x δ 4 ? Add s = ?3 to desired pole list Plist=[roots([1 2*.6*3 3^2])’,-.25,-3]; At2=[Asp2 Bsp2(:,1);zeros(1,3) -4];Bt2=[zeros(3,1);4]; Kt=place(At2,Bt2,Plist); step(ss(At2-Bt2*Kt2,Bt2,[0 0 1 0],0),35) ? ? ? ? ? ? w q θ x δ ? ? ? ? u = 0.0011 ?3.4617 ?4.9124 0.5273 0 5 10 15 20 25 30 35?0.25 ?0.2 ?0.15 ?0.1 ?0.05 0 Step Response Time (sec) q rads ? No problem working with larger systems with state space tools ? Main control issue is ?nding “good” locations for closed-loop poles Fall 2004 16.333 9–12 Estimators/Observers Problem: So far we have assumed that we have full access to the ? state x(t) when we designed our controllers. – Most often all of this information is not available. ? Usually can only feedback information that is developed from the sensors measurements. – Could try “output feedback” ? u = Kx ? u = Ky – Same as the proportional feedback we looked at at the beginning of the root locus work. – This type of control is very di?cult to design in general. ? Alternative approach: Develop a replica of the dynamic system that provides an “estimate” of the system states based on the mea- sured output of the system. ? New plan: 1. Develop estimate of x(t) that will be called x?(t). 2. Then switch from u = ?Kx(t) to u = ?Kx?(t). ? Two key questions: – How do we ?nd x?(t)? – Will this new plan work? ? Fall 2004 16.333 9–13 Estimation Schemes ? Assume that the system model is of the form: x˙ = Ax + Bu , x(0) unknown y = Cx where 1. A, B, and C are known. 2. u(t) is known 3. Measurable outputs are y(t) from C = I ? Goal: Develop a dynamic system whose state x?(t) = x(t) for all time t ≥ 0. Two primary approaches: – Open-loop. – Closed-loop. Fall 2004 16.333 9–14 Open-loop Estimator ? Given that we know the plant matrices and the inputs, we can just perform a simulation that runs in parallel with the system x? ˙ (t) = Ax? + Bu(t) ? Then x?(t) ≡ x(t) ? t provided that x?(0) = x(0) ? Major Problem: We do not know x(0) System A,B,C ? x(t) y(t) // u(t) // // Observer A,B,C ? ?x(t) ?y(t) // ? Analysis of this case. Start with: x˙(t) = Ax + Bu(t) x? ˙ (t) = Ax? + Bu(t) ? De?ne the estimation error: x?(t) = x(t) ? x?(t). – Now want x?(t) = 0 ? t. – But is this realistic? ? Fall 2004 16.333 9–15 ? Subtract to get: d (x ?x?) = A(x ?x?) ? x? ˙ (t) = Ax? dt which has the solution x?(t) = e At x?(0) – Gives the estimation error in terms of the initial error. ? Does this guarantee that x? = 0 ? t? Or even that x? → 0 as t →∞? (which is a more realistic goal). – Response is ?ne if x?(0) = 0. But what if x?(0) = 0? ? If A stable, then x? → 0 as t →∞, but the dynamics of the estima- tion error are completely determined by the open-loop dynamics of the system (eigenvalues of A). – Could be very slow. – No obvious way to modify the estimation error dynamics. ? Open-loop estimation does not seem to be a very good idea. Fall 2004 16.333 9–16 Closed-loop Estimator ? Obvious way to ?x the problem is to use the additional information available: – How well does the estimated output match the measured output? Compare: y = Cx with y? = Cx? – Then form y? = y ? y? ≡ Cx? System A,B,C → x(t) y(t) // +   u(t) // // L Observer A,B,C → ?x(t) ?y(t) // ? OO ? Approach: Feedback y? to improve our estimate of the state. Basic form of the estimator is: x? ˙ (t) = Ax?(t) + Bu(t) + Ly?(t) y?(t) = Cx?(t) where L is a user selectable gain matrix. ? Analysis: ˙ x? = x˙ ? x? ˙ = [Ax + Bu] [Ax? + Bu + L(y ? y?)]? = A(x ? x?) ? L(Cx ? Cx?) = Ax? ? LCx? = (A ? LC)?x Fall 2004 16.333 9–17 ? So the closed-loop estimation error dynamics are now x? ˙ = (A ?LC)?x with solution x?(t) = e (A?LC)t x?(0) ? Bottom line: Can select the gain L to attempt to improve the convergence of the estimation error (and/or speed it up). – But now must worry about observability of the system model. ? Note the similarity: – Regulator Problem: pick K for A ?BK 3 Choose K ∈R 1×n (SISO) such that the closed-loop poles det(sI ?A + BK) = Φ c (s) are in the desired locations. – Estimator Problem: pick L for A ?LC 3 Choose L ∈R n×1 (SISO) such that the closed-loop poles det(sI ?A + LC) = Φ o (s) are in the desired locations. ? These problems are obviously very similar – in fact they are called dual problems. Fall 2004 16.333 9–18 Estimation Gain Selection ? For regulation, were concerned with controllability of (A, B) For a controllable system we can place the eigenvalues of A ? BK arbitrarily. ? For estimation, were concerned with observability of pair (A, C). For an observable system we can place the eigenvalues of A ? LC arbitrarily. ? Test using the observability matrix: ?? rank M o ? rank ? ? ? ? ? ? C CA CA 2 . . . CA n?1 ? ? ? ? ? ? = n ? The procedure for selecting L is very similar to that used for the regulator design process. Fall 2004 16.333 9–19 ? One approach: – Note that the poles of (A ? LC) and (A ? LC) T are identical. – Also we have that (A ? LC) T = A T ? C T L T – So designing L T for this transposed system looks like a standard regulator problem (A ? BK) where A A T ? B C T ? K L T ? So we can use K e = acker(A T , C T , P ) , L ≡ K T e ? Note that the estimator equivalent of Ackermann’s formula is that ? ? L = Φ e (s)M ?1 o ? ? ? ? 0 . . . 0 1 ? ? ? ? ? ? ? ? ? ? ? ? Fall 2004 16.333 9–20 Simple Estimator Example ? Simple system ? ? ? ? ? ? ?1 1.5 1 .5 A = , B = , x(0) = ?0 1 ?2 0 ?1 C = 1 0 , D = 0 – Assume that the initial conditions are not well known. – System stable, but λ max (A) = ?0.18 – Test observability: C rank CA 1 0 = rank ?1 1.5 ? Use open and closed-loop estimators. Since the initial conditions are 0 not well known, use x?(0) = 0 ? Open-loop estimator: x? ˙ = Ax? + Bu y? = Cx? ? Closed-loop estimator: x? ˙ = Ax? + Bu + Ly? = Ax? + Bu + L(y ? y?) = (A ? LC)?x + Bu + Ly y? = Cx? – Which is a dynamic system with poles given by λ i (A ? LC) and which takes the measured plant outputs as an input and generates an estimate of x. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Fall 2004 16.333 9–21 ? Typically simulate both systems together for simplicity ? Open-loop case: x˙ = Ax + Bu y = Cx x? ˙ = Ax? + Bu y? = Cx? ? ? ?0.5 ?1 x˙ A 0 x ? ? ? ? ? ? ? ? B x(0) = + u , = x?(0) ? x? ˙ 0 A x? B y ? C 0 x = y? 0 C x? ? Closed-loop case: x˙ = Ax + Bu x? ˙ = (A ? LC)?x + Bu + LCx x˙ A 0 x B = u B ? x? ˙ LC A ? LC x? + ? Example uses a strong u(t) to shake things up 0 0 Fall 2004 16.333 9–22 0 0.5 1 1.5 2 2.5 3 3.5 4?1 ?0.5 0 0.5 1 Open loop estimator states time x1 x2 0 0.5 1 1.5 2 2.5 3 3.5 4?1 ?0.5 0 0.5 1 estimation error time Figure 1: Open-loop estimator. Estimation error converges to zero, but very slowly. 0 0.5 1 1.5 2 2.5 3 3.5 4?1 ?0.5 0 0.5 1 Closed?loop estimator states time x1 x2 0 0.5 1 1.5 2 2.5 3 3.5 4?1 ?0.5 0 0.5 1 estimation error time Figure 2: Closed-loop estimator. Convergence looks much better. ? Fall 2004 16.333 9–23 Aircraft Estimation Example ? Take Short period model and assume that we can measure q. Can we estimate the motion associated with the short period mode? x˙ sp = A sp x sp ? + B sp u y = 0 1 x sp – Take x sp (0) = [?0.5; ?0.05] T ? System stable, so could use an open loop estimator ? For closed-loop estimator, put desired poles at ?3, ?4 ? For the various dynamics models as before Csp=[0 1]; % sense q Ke=place(Asp’,Csp’,[-3 -4]);Le=Ke’; 0 2 4 6 8 10?20 ?15 ?10 ?5 0 5 10 x1 Open?loop estimator time 0 2 4 6 8 10 ?0.1 ?0.05 0 0.05 0.1 x1 time 0 2 4 6 8 10?8 ?6 ?4 ?2 0 2 4 x1 error time 0 2 4 6 8 10 ?0.05 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 0.02 x2 error time 0 2 4 6 8 10?20 ?15 ?10 ?5 0 5 10 15 x1 Closed?loop estimator time 0 2 4 6 8 10 ?0.1 ?0.05 0 0.05 0.1 x1 time 0 2 4 6 8 10?20 ?15 ?10 ?5 0 x1 error time 0 2 4 6 8 10 ?0.05 ?0.04 ?0.03 ?0.02 ?0.01 0 0.01 x2 error time Figure 3: Closed-loop estimator. Convergence looks much better. ? As expected, the OL estimator does not do well, but the closed-loop one converges nicely Fall 2004 16.333 9–24 Where to put the Estimator Poles? ? Location heuristics for poles still apply – use Bessel, ITAE, ... – Main di?erence: probably want to make the estimator faster than you intend to make the regulator – should enhance the control, which is based on x?(t). – ROT: Factor of 2–3 in the time constant ζω n associated with the regulator poles. ? Note: When designing a regulator, were concerned with “band- width” of the control getting too high ? often results in control commands that saturate the actuators and/or change rapidly. Di?erent concerns for the estimator: ? – Loop closed inside computer, so saturation not a problem. – However, the measurements y are often “noisy”, and we need to be careful how we use them to develop our state estimates. ? High bandwidth estimators tend to accentuate the e?ect of sens- ing noise in the estimate. – State estimates tend to “track” the measurements, which are ?uc- tuating randomly due to the noise. ? Low bandwidth estimators have lower gains and tend to rely more heavily on the plant model – Essentially an open-loop estimator – tends to ignore the measure- ments and just uses the plant model. Fall 2004 16.333 9–25 ? Can also develop an optimal estimator for this type of system. – Which is apparently what Kalman did one evening in 1958 while taking the train from Princeton to Baltimore... – Balances e?ect of the various types of random noise in the system on the estimator: x˙ = Ax + Bu + B w w y = Cx + v where: 3 w: “process noise” – models uncertainty in the system model. 3 v: “sensor noise” – models uncertainty in the measurements. Final Thoughts ? Note that the feedback gain L in the estimator only stabilizes the estimation error. – If the system is unstable, then the state estimates will also go to ∞, with zero error from the actual states. ? Estimation is an important concept of its own. – Not always just “part of the control system” – Critical issue for guidance and navigation system ? More complete discussion requires that we study stochastic processes and optimization theory. ? Estimation is all about which do you trust more: your mea- surements or your model. ? ? ? ? ? ? ? ? Fall 2004 16.333 9–26 Combined Regulator and Estimator ? As advertised, we can change the previous control u = ?Kx to the new control u = ?Kx? (same K). We now have x˙ = Ax + Bu y = Cx x? ˙ = Ax? + Bu + L(y ? y?) y? = Cx? with closed-loop dynamics ? ? ? ?? ? x˙ A ?BK x x? ˙ = LC A ? BK ? LC x? ? x˙ cl = A cl x cl ? Not obvious that this system will even be stable: λ i (A cl ) < 0? ? To analyze, introduce x? = x ? x?, and the similarity transform I 0 T = = T ?1 I ?I x x ? Rewrite the dynamics in terms of the state = T x? x? A cl ? T ?1 A cl T ≡ A cl and when you work through the math, you get A cl = A ? BK BK !!! 0 A ? LC Fall 2004 16.333 9–27 ? Absolutely key points: 1. λ i (A cl ) ≡ λ i (A cl ) why? 2. A cl is block upper triangular, so can ?nd poles by inspection: det(sI ? A cl ) = det(sI ? (A ? BK)) · det(sI ? (A ? LC)) The closed-loop poles of the system consist of the union of the regulator and estimator poles ? So we can design the estimator and regulator separately with con?- dence that combination of the two will work VERY well. ? Compensator is a combination of the estimator and regulator. x? ˙ = Ax? + Bu + L(y ? y?) = (A ? BK ? LC)?x + Ly u = ?Kx? ? x˙ c = A c x c + B c y u = ?C c x c – Keep track of this minus sign. We need one in the feedback path, but we can move it around to suit our needs. Fall 2004 16.333 9–28 ? Let G c (s) be the compensator transfer function where u = ?C c (sI ? A c ) ?1 B c = ?G c (s) y = ?K(sI ? (A ? BK ? LC)) ?1 L so by my de?nition, ? u = ?G c y ≡ G c (?y) ? Reason for making the de?nition is that when we implement the controller, we often do not just feedback ?y(t), but instead have to include a reference command r(t) – Use servo approach and feed back e(t) = r(t)? y(t) instead r // e // G c (s) u // G(s) y // ? OO – So now u = G c e = G c (r ? y). – And if r = 0, then we still have u = G c (?y) ? Important points: – Closed-loop system will be stable, but the compensator dynamics need not be. – Often very simple and useful to provide classical interpretations of the compensator dynamics G c (s). ? ? ? ? Fall 2004 16.333 9–29 ? Mechanics of closing the loop G(s) : ˙x = Ax + Bu y = Cx G c (s) : ˙x c = A c x c + B c e u = C c x c and e = r ? y, u = G c e, y = Gu. ? Loop dynamics L = GG c ? y = L(s)e x˙ = Ax + Bu = Ax + BC c x c x˙ c = A c x c + B c e ? ? ? ?? ? ? ? x˙ A BC c x 0 = + e x˙ c 0 A c x c B c ? ? x y = C 0 x c ? Now form the closed-loop dynamics by inserting e = r ? y ? ? ? ?? ? ? ?? ? ?? x˙ A BC c x 0 ? ? x = x˙ c 0 A c x c + B c r ? C 0 x c ? ?? ? ? ? A BC c x 0 = + r ?B c C A c x c B c ? ? x y = C 0 x c ? ? ? ? ? ? ? ? Fall 2004 16.333 9–30 Performance Issue ? Often ?nd with state space controllers that the DC gain of the closed loop system is not 1. So y = r in steady state. ? Relatively simple ?x is to modify the original controller with scalar N u = r ? Kx ? u = N r ? Kx ? Closed-loop system on page 5 becomes x˙ = Ax + B(N r ? Kx) = A cl x + BN r G cl (s) = C(sI?A cl ) ?1 BN y = Cx – Analyze steady state step response ? y ss = G cl (0)r step G cl (0) = C(?A cl ) ?1 BN – And pick N so that G cl (0) = 1 ? N = 1 (C(?A cl ) ?1 B) ? A bit more complicated with a combined estimator and regulator – One simple way (not the best) of achieving a similar goal is to add N to r and force G cl (0) = 1 – Now the closed-loop dynamics on page 29 become: ? ? ? ? ? x˙ x = A cl + B cl N r x˙ c 1 N = → (C cl (?A cl ) ?1 B cl ) x c ? ? ? ? x y = C cl x c ? Note that this ?xes the steady state tracking error problems, but in my experience can create strange transients (often NMP). ? ? ? ? ? ? Fall 2004 16.333 9–31 Example: Compensator Design 1 x˙ = Ax + Bu G(s) = s 2 + s + 1 ? y = Cx where ? 0 1 ? ? 0 ? ? ? A = B = C = 1 0 ?1 ?1 1 ? Regulator: Want regulator poles to have a time constant of τ c = 1/(ζω n ) = 0.25 sec ? λ(A ? BK r ) = ?4 ± 4j which can be found using place or acker K_r=acker(a,b,[-4+4*j;-4-4*j]); to give K r = [ 31 7 ] ? Estimator: want the estimator poles to be faster, so use τ e = 1/(ζω n ) = 0.1 sec. Use real poles, ? λ(A ? L e C) = ?10 L_e=acker(a’,c’,[-10 -10]’)’; 19 which gives L e = 80 ? Form compensator G c (s) ac=a-b*K_r-L_e*c;bc=L_e;cc=K_r;dc=0; A c = ?19 1 19 ? ? B c = C c = 31 7 ?112 80?8 (s + 2.5553) u G c (s) = 1149 = s 2 + 27s + 264 e Low frequency zero, with higher frequency poles (like a lead) Fall 2004 16.333 9–32 10?1 100 101 102 100 101 102 Freq (rad/sec) Mag Plant G Compensator Gc 10?1 100 101 102 ?200 ?150 ?100 ?50 0 50 Freq (rad/sec) Phase (deg) Plant G Compensator Gc Figure 4: The compensator does indeed look like a high frequency lead (ampli?cation from 2–16 rad/sec). Plant pretty simple looking. Fall 2004 16.333 9–33 10?1 100 101 102 100 101 102 Freq (rad/sec) Mag Loop L 10?1 100 101 102 ?250 ?200 ?150 ?100 ?50 0 Freq (rad/sec) Phase (deg) Figure 5: The loop transfer function L = G c G shows a slope change around ω c = 5 rad/sec due to the e?ect of the compensator. Signi?cant gain and phase margins. Fall 2004 16.333 9–34 ?150 ?100 ?50 0 50 Magnitude (dB) 10?2 10?1 100 101 102 103 ?270 ?225 ?180 ?135 ?90 ?45 0 Phase (deg) Bode Diagram Gm = 14.3 dB (at 14.9 rad/sec) , Pm = 45.8 deg (at 4.84 rad/sec) Frequency (rad/sec) Figure 6: Quite signi?cant gain and phase margins. Fall 2004 16.333 9–35 ?40 ?35 ?30 ?25 ?20 ?15 ?10 ?5 0 5 10 ?25 ?20 ?15 ?10 ?5 0 5 10 15 20 25 Root Locus Real Axis Imaginary Axis Figure 7: Freeze the compensator poles and zeros and draw a root locus versus an ? additional plant gain α, G(s) ? G(s) = (s 2 + α s+1) . Note location of the closed-loop poles!! Fall 2004 16.333 9–36 10?1 100 101 102 10?2 10?1 100 101 102 Freq (rad/sec) Mag Factor of N=1.0899 applied to Closed loop Plant G closed?loop Gcl Figure 8: Closed-loop transfer – system bandwidth has increased substantially. Fall 2004 16.333 9–37 Estimator Design (est1.m) 1 clear all 2 close all 3 figure(1);clf 4 set(gcf,’DefaultLineLineWidth’,2) 5 set(gcf,’DefaultlineMarkerSize’,10) 6 figure(2);clf 7 set(gcf,’DefaultLineLineWidth’,2) 8 set(gcf,’DefaultlineMarkerSize’,10) 9 10 load b747 % get A B Asp Bsp 11 Csp=[0 1]; % sense q 12 13 Ke=place(Asp’,Csp’,[-3 -4]);Le=Ke’; 14 15 xo=[-.5;-.05]; % start somewhere 16 17 t=[0:.01:10];N=floor(.15*length(t)); 18 % hit on the system with an input 19 %u=0;u=[ones(15,1);-ones(15,1);ones(15,1)/2;-ones(15,1)/2;zeros(41,1)]/5; 20 u=0;u=[ones(N,1);-ones(N,1);ones(N,1)/2;-ones(N,1)/2]/20; 21 u(length(t))=0; 22 23 [y,x]=lsim(Asp,Bsp,Csp,0,u,t,xo); 24 plot(t,y) 25 26 % closed-loop estimator 27 % hook both up so that we can simulate them at the same time 28 % bigger state = state of the system then state of the estimator 29 A_cl=[Asp zeros(size(Asp));Le*Csp Asp-Le*Csp]; 30 B_cl=[Bsp;Bsp]; 31 C_cl=[Csp zeros(size(Csp));zeros(size(Csp)) Csp]; 32 D_cl=zeros(2,1); 33 34 % note that we start the estimators at zero, since that is 35 % our current best guess of what is going on (i.e. we have no clue :-) ) 36 % 37 [y_cl,x_cl]=lsim(A_cl,B_cl,C_cl,D_cl,u,t,[xo;0;0]); 38 figure(1) 39 subplot(221) 40 plot(t,x_cl(:,[1]),t,x_cl(:,[3]),’--’) 41 ylabel(’x1’);title(’Closed-loop estimator’);xlabel(’time’);grid 42 subplot(222) 43 plot(t,x_cl(:,[2]),t,x_cl(:,[4]),’--’) 44 ylabel(’x1’);xlabel(’time’);grid 45 subplot(223) 46 plot(t,x_cl(:,[1])-x_cl(:,[3])) 47 ylabel(’x1 error’);xlabel(’time’);grid 48 subplot(224) 49 plot(t,x_cl(:,[2])-x_cl(:,[4])) 50 ylabel(’x2 error’);xlabel(’time’);grid 51 print -depsc spest_cl.eps 52 jpdf(’spest_cl’) 53 54 % open-loop estimator 55 % hook both up so that we can simulate them at the same time 56 % bigger state = state of the system then state of the estimator 57 A_ol=[Asp zeros(size(Asp));zeros(size(Asp)) Asp]; 58 B_ol=[Bsp;Bsp]; 59 C_ol=[Csp zeros(size(Csp));zeros(size(Csp)) Csp]; 60 D_ol=zeros(2,1); 61 62 [y_ol,x_ol]=lsim(A_ol,B_ol,C_ol,D_ol,u,t,[xo;0;0]); 63 figure(2) 64 subplot(221) 65 plot(t,x_ol(:,[1]),t,x_ol(:,[3]),’--’) 66 ylabel(’x1’);title(’Open-loop estimator’);xlabel(’time’);grid Fall 2004 16.333 9–38 67 subplot(222) 68 plot(t,x_ol(:,[2]),t,x_ol(:,[4]),’--’) 69 ylabel(’x1’);xlabel(’time’);grid 70 subplot(223) 71 plot(t,x_ol(:,[1])-x_ol(:,[3])) 72 ylabel(’x1 error’);xlabel(’time’);grid 73 subplot(224) 74 plot(t,x_ol(:,[2])-x_ol(:,[4])) 75 ylabel(’x2 error’);xlabel(’time’);grid 76 print -depsc spest_ol.eps 77 jpdf(’spest_ol’) 78 Fall 2004 16.333 9–39 Regular/Estimator Design (reg est.m) 1 % Combined estimator/regulator design for a simple system 2 % G= 1/(s^2+s+1) 3 % 4 % Jonathan How 5 % Fall 2004 6 % 7 close all;clear all 8 for ii=1:5 9 figure(ii);clf;set(gcf,’DefaultLineLineWidth’,2);set(gcf,’DefaultlineMarkerSize’,10) 10 end 11 12 a=[0 1;-1 -1];b=[0 1]’;c=[1 0];d=0; 13 k=acker(a,b,[-4+4*j;-4-4*j]); 14 l=acker(a’,c’,[-10 -10]’)’; 15 % 16 % For state space for G_c(s) 17 % 18 ac=a-b*k-l*c;bc=l;cc=k;dc=0; 19 20 G=ss(a,b,c,d); 21 Gc=ss(ac,bc,cc,dc); 22 23 f=logspace(-1,2,400); 24 g=freqresp(G,f*j);g=squeeze(g); 25 gc=freqresp(Gc,f*j);gc=squeeze(gc); 26 27 figure(1);clf 28 subplot(211) 29 loglog(f,abs(g),f,abs(gc),’--’);axis([.1 1e2 .2 1e2]) 30 xlabel(’Freq (rad/sec)’);ylabel(’Mag’) 31 legend(’Plant G’,’Compensator Gc’);grid 32 subplot(212) 33 semilogx(f,180/pi*angle(g),f,180/pi*angle(gc),’--’); 34 axis([.1 1e2 -200 50]) 35 xlabel(’Freq (rad/sec)’);ylabel(’Phase (deg)’);grid 36 legend(’Plant G’,’Compensator Gc’) 37 38 L=g.*gc; 39 40 figure(2);clf 41 subplot(211) 42 loglog(f,abs(L),[.1 1e2],[1 1]);axis([.1 1e2 .2 1e2]) 43 xlabel(’Freq (rad/sec)’);ylabel(’Mag’) 44 legend(’Loop L’); 45 grid 46 subplot(212) 47 semilogx(f,180/pi*phase(L.’),[.1 1e2],-180*[1 1]); 48 axis([.1 1e2 -290 0]) 49 xlabel(’Freq (rad/sec)’);ylabel(’Phase (deg)’);grid 50 % 51 % loop dynamics L = G Gc 52 % 53 al=[a b*cc;zeros(2) ac]; 54 bl=[zeros(2,1);bc]; 55 cl=[c zeros(1,2)]; 56 dl=0; 57 figure(3) 58 rlocus(al,bl,cl,dl) 59 % 60 % closed-loop dynamics 61 % unity gain wrapped around loop L 62 % 63 acl=al-bl*cl;bcl=bl;ccl=cl;dcl=d; 64 65 N=inv(ccl*inv(-acl)*bcl) 66 Fall 2004 16.333 9–40 67 hold on;plot(eig(acl),’d’);hold off 68 grid 69 % 70 % closed-loop freq response 71 % 72 Gcl=ss(acl,bcl*N,ccl,dcl); 73 gcl=freqresp(Gcl,f*j);gcl=squeeze(gcl); 74 75 figure(4);clf 76 loglog(f,abs(g),f,abs(gcl),’--’); 77 axis([.1 1e2 .01 1e2]) 78 xlabel(’Freq (rad/sec)’);ylabel(’Mag’) 79 legend(’Plant G’,’closed-loop G_{cl}’);grid 80 title([’Factor of N=’,num2str(N),’ applied to Closed loop’]) 81 82 figure(5);clf 83 margin(al,bl,cl,dl) 84 85 figure(1);orient tall;print -depsc reg_est1.eps 86 jpdf(’reg_est1’) 87 figure(2);orient tall;print -depsc reg_est2.eps 88 jpdf(’reg_est2’) 89 figure(3);print -depsc reg_est3.eps 90 jpdf(’reg_est3’) 91 figure(4);print -depsc reg_est4.eps 92 jpdf(’reg_est4’) 93 figure(5);print -depsc reg_est5.eps 94 jpdf(’reg_est5’)