Anti-Windup

Topic #23 16.30/31 Feedback Control Systems Analysis of Nonlinear Systems • Anti-windup • Notes developed in part from 16.30 Estimation and Con­ trol of Aerospace Systems, Lecture 21: Lyapunov Stability Theory by Prof. Frazzoli Fall 2010 16.30/31 23–2 Nonlinear System Analysis • Example: Car cruise control1 F Fg mg • Equation of motion in the direction parallel to the road surface: m dv = Feng + Faero + Ffrict + Fg . dt where 1 Faero = − ρCdAv · |v|, 2 Fg = −mg sin(θ), Ffrict = −mgCr cos(θ)sgn(v). 200 200 180 180 Torque T [Nm] Torque T [Nm] • Engine model 160 140 120 100 0 160 140 120 200 400 600 Angular velocity � [rad/s] n=5 100 0 n=1 n=2 n=3 20 40 Velocity v [m/s] � • Engine torque (at full throttle): Tω = Tm 1 − β � ω ωm n=4 60 −1 where ω = nr v = αnv, n is gear ratio, and r wheel radius. • The engine driving force can hence be written as Feng = αnT (αnv)u, 1 The example is taken from ˚ Astr¨ om and Murray: Feedback Systems, 2008 November 21, 2010 0 ≤ u ≤ 1. �2 � , Fall 2010 16.30/31 23–3 Jacobian Linearization • Any (feasible) speed corresponds to an equilibrium point. • Choose a reference speed vref > 0, and solve for dv/dt = 0 with respect to u, assuming a horizontal road (θ = 0). 1 0 = αnT (αnv̄)ū − ρCdAv̄ 2 − mgCr 2 i.e., + mgCr . αnT (αnv̄) 1 2 2 ρCd Av̄ ū = • Linearized system (ξ = v − v̄, η = u − ū): � � � � 1 d 1 ∂T (αnv) � ū − ρC ξ= αn αnT (αnv̄) η Av̄ ξ + d dt m ∂v �v̄ m � �� � � �� � Bdyn Adyn • Example: numerical values • Let us use the following numerical values (all units in SI): Tm = 190, β = 0.4, ωm = 420, α5 = 10, Cr = 0.01, m = 1500, g = 9.81, ρ = 1.2, CdA = 0.79. • For v̄ = 25 (90 km/h, or 55 mph), we get ū = 0.2497. • The linearization yields: Adyn = −0.0134, ⇒ November 21, 2010 G(s) = Bdyn = 1.1837 1.1837 s + 0.0134 Fall 2010 16.30/31 23–4 Cruise control design • A proportional controller would stabilize the closed-loop system. • Assume we want to maintain the commanded speed (cruise con­ trol): we need to add an integrator. • A PI controller will work, e.g., C(s) = 1.5 s+1 s • Linear control/model response to hill at specified angle between 5 and 25 sec November 21, 2010 Fall 2010 16.30/31 23–5 Nonlinear Simulation • Check with BOTH linear AND nonlinear simulation Signal 1 Signal Builder -KGain 2 vbar sin Trigonometric Function ubar -K- f(u) Fg engine Gain 4 To Workspace Clock Faero 1 1 -C­ Ffrict F del V Constant 3 s control -KGain 1.5s+1.5 0 Car1 f(u) Saturation 1 s v Integrator Add Product v -Kvbar 1.5s+1.5 Gain 3 Bdyn s s-Adyn control 1 Transfer Fcn Scope ubar November 21, 2010 Fall 2010 16.30/31 23–6 Review • (Jacobian) linearization: • Find the desired equilibrium condition (state and control). • Linearize the non-linear model around the equilibrium. • Control design: • Design a linear compensator for the linear model. • If the linear system is closed-loop stable, so will be the nonlinear system—in a neighborhood of the equilibrium. • Check in a (nonlinear) simulation the robustness of your design with respect to “typical” deviations. November 21, 2010 Fall 2010 16.30/31 23–7 Effects of the saturation • What if the slope is a little steeper (say 4.85 degrees)? • What is wrong? • Systems experiencing Integrator wind-up • Once the input saturates, the integral of the error keeps increasing. • When the error decreases, the large integral value prevents the controller from resuming “normal operations” quickly (the integral error must decrease first!) – so the response is delayed • Idea: once the input saturates, stop integrating the error (can’t do much about it anyway!) November 21, 2010 Fall 2010 16.30/31 23–8 Anti-windup Logic • One option is the following logic for the integral gain: � KI if the input does not saturate; KI = 0 if the input saturates • Another option is the following: • Compare the actual input and the commanded input. • If they are the same, the saturation is not in effect. • Otherwise, reduce the integral error by a constant times the dif­ ference. � uint = Ki (e + Kaw (sat(u) − u)) dt • With this choice, under saturation, the actuator/commanded differ­ ence is fedback to the integrator to so that eact = sat(u) − u tends to zero • Implies that the controller output is kept close to saturation limit. • The controller output will then change as soon as the error changes sign, thus avoiding the integral windup. • If there is no saturation, the anti-windup scheme has no effect. November 21, 2010 Fall 2010 16.30/31 23–9 Anti-windup Scheme Signal 1 Gain 2 Signal Builder vbar sin -K- Trigonometric Function f(u) -K- Fg Car2 engine Gain 4 f(u) To Workspace Clock ubar Faero Constant 1.5 -C1 0 -KGain 1 s v Integrator Constant 1 Transfer Fcn 2 1.5 Constant 3 F Ffrict Saturation Add Product s Transfer Fcn 3 -KGain 1 -K1.5s+1.5 vbar Gain 3 Bdyn s s-Adyn control 1 Transfer Fcn Scope ubar • Response to a 4.85 degree slope • Anti-windup compensator avoids the velocity overshoot by preventing the error buildup in the integral term November 21, 2010 Fall 2010 16.30/31 23–10 Anti-windup Summary • Anti-wind up schemes guarantee the stability of the compensator when the (original) feedback loop is effectively opened by the satura­ tion. • Prevent divergence of the integral error when the control cannot keep up with the reference. • Maintain the integral errors “small.” November 21, 2010 Fall 2010 16.30/31 23–11 Code: Car Setup 1 2 3 4 5 set(0, 'DefaultAxesFontSize', 12, 'DefaultAxesFontWeight','demi') set(0, 'DefaultTextFontSize', 12, 'DefaultTextFontWeight','demi') set(0,'DefaultAxesFontName','arial') set(0,'DefaultAxesFontSize',12) set(0,'DefaultTextFontName','arial') 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 clear all %close all % Parameters for defining the system dynamics theta=0; gear=5; alpha = [40, 25, 16, 12, 10]; % gear ratios Tm = 190; % engine torque constant, Nm wm = 420; % peak torque rate, rad/sec beta = 0.4; % torque coefficient Cr = 0.01; % coefficient of rolling friction rho = 1.2; % density of air, kg/mˆ3 A = 2.4; % car area, mˆ2 Cd = 0.79/A; % drag coefficient g = 9.8; % gravitational constant m=1500; % mass v=25; 23 24 25 26 27 % Compute the torque produced by the engine, Tm omega = alpha(gear) * v; % engine speed torque = Tm * ( 1 − beta * (omega/wm − 1)ˆ2 ); F = alpha(gear) * torque; 28 29 30 31 32 33 % Compute the external forces on the vehicle Fr = m * g * Cr; % Rolling friction Fa = 0.5 * rho * Cd * A * vˆ2; % Aerodynamic drag Fg = m * g * sin(theta); % Road slope force Fd = Fr + Fa + Fg; % total deceleration 34 35 36 ubar=(Fa+Fr)/(F) vbar=v; 37 38 39 40 dTdv=Tm*(−2*beta*(alpha(gear)*vbar/wm−1)*(alpha(gear)/wm)) Adyn=(alpha(gear)*dTdv*ubar−rho*Cd*A*vbar)/m; Bdyn=F/m; 41 42 43 44 45 46 47 48 49 Hillangle=4; % non−saturating angle sim('cruise control') figure(4); subplot(211); plot(Car1(:,5),Car1(:,[2]),'LineWidth',2);xlabel('Time');ylabel('Speed') title(['Hill response, angle=',num2str(Hillangle)]) subplot(212); plot(Car1(:,5),Car1(:,[4]),'LineWidth',2);xlabel('Time');ylabel('Throttle') 50 51 52 53 54 55 56 57 figure(1); subplot(211); plot(Car1(:,5),Car1(:,[1 2]),'LineWidth',2);xlabel('Time');ylabel('Speed') legend('NL','Lin','Location','SouthEast'); title(['Hill response, angle=',num2str(Hillangle)]) subplot(212); plot(Car1(:,5),Car1(:,[3 4]),'LineWidth',2);xlabel('Time');ylabel('Throttle') 58 59 60 61 62 63 64 65 66 67 68 Antiwindup gain=0; % us esame code with and without AWup Hillangle=4.85; sim('cruise control awup') figure(2); subplot(211); plot(Car2(:,5),Car2(:,[1 2]),'LineWidth',2);xlabel('Time');ylabel('Speed') legend('NL','Lin','Location','SouthEast'); title(['Hill response, angle=',num2str(Hillangle)]) subplot(212); plot(Car2(:,5),Car2(:,[3 4]),'LineWidth',2);xlabel('Time');ylabel('Throttle') 69 70 Car2 no AW=Car2; 71 72 73 74 Antiwindup gain=5; sim('cruise control awup') figure(3); November 21, 2010 Fall 2010 75 76 77 78 79 80 16.30/31 23–12 subplot(211); plot(Car2(:,5),Car2(:,[1 2]),'LineWidth',2);xlabel('Time');ylabel('Speed') legend('NL','Lin','Location','SouthEast'); title(['Hill response, angle=',num2str(Hillangle),' With Anti−windup gain=',num2str(Antiwindup gain)]) subplot(212); plot(Car2(:,5),Car2(:,[3 4]),'LineWidth',2);xlabel('Time');ylabel('Throttle') 81 82 83 84 85 86 87 88 89 90 91 92 figure(5); subplot(211); plot(Car2 no AW(:,5),Car2 no AW(:,[6 7]),Car2 no AW(:,5),10* Car2 no AW(:,[8]),'LineWidth',2);xlabel('Time');yla legend('sat(u)','u','10* e {vel}','Location','NorthEast'); axis([4 20 −0.5 5]) title(['Hill response no AW']) subplot(212); plot(Car2(:,5),Car2(:,[6 7]),Car2(:,5),10*Car2(:,[8]),'LineWidth',2);xlabel('Time');ylabel('Windup error') legend('sat(u)','u','10* e {vel}','Location','NorthEast'); title(['Hill response with AW']) axis([4 20 −0.5 5]) 93 94 95 96 97 98 print print print print print −dpng −dpng −dpng −dpng −dpng −r300 −r300 −r300 −r300 −r300 November 21, 2010 −f1 −f2 −f3 −f4 −f5 AW1.png AW2.png AW3.png AW4.png AW5.png