Robust MIMO Controller for Two-Loop Autopilot
This example shows how to design a robust controller for a two-loop autopilot that controls the pitch rate and vertical acceleration of an airframe. The controller is robust against gain and phase variations in the multichannel feedback loop.
Linearized Airframe Model
The airframe dynamics and the autopilot are modeled in Simulink®. See Tuning of a Two-Loop Autopilot for more information about this model and for additional design and tuning options.
open_system('rct_airframe2')
As in the example Tuning of a Two-Loop Autopilot, trim the airframe for and . The trim condition corresponds to zero normal acceleration and pitching moment ( and steady). Use findop
to compute the corresponding operating condition. Then, linearize the airframe model at this trim condition.
opspec = operspec('rct_airframe2'); % Specify trim condition % Xe,Ze: known, not steady opspec.States(1).Known = [1;1]; opspec.States(1).SteadyState = [0;0]; % u,w: known, w steady opspec.States(3).Known = [1 1]; opspec.States(3).SteadyState = [0 1]; % theta: known, not steady opspec.States(2).Known = 1; opspec.States(2).SteadyState = 0; % q: unknown, steady opspec.States(4).Known = 0; opspec.States(4).SteadyState = 1; % controller states unknown, not steady opspec.States(5).SteadyState = [0;0]; op = findop('rct_airframe2',opspec); G = linearize('rct_airframe2','rct_airframe2/Airframe Model',op); G.InputName = 'delta'; G.OutputName = {'az','q'};
Operating point search report: --------------------------------- opreport = Operating point search report for the Model rct_airframe2. (Time-Varying Components Evaluated at time t=0) Operating point specifications were successfully met. States: ---------- Min x Max dxMin dx dxMax ___________ ___________ ___________ ___________ ___________ ___________ (1.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/Position 0 0 0 -Inf 984 Inf -3047.9999 -3047.9999 -3047.9999 -Inf 0 Inf (2.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/Theta 0 0 0 -Inf -0.0097235 Inf (3.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/U,w 984 984 984 -Inf 22.6897 Inf 0 0 0 0 2.4588e-11 0 (4.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/q -Inf -0.0097235 Inf 0 -4.0169e-16 0 (5.) rct_airframe2/MIMO Controller -Inf 0.00065361 Inf -Inf -0.0089973 Inf -Inf 2.6167e-19 Inf -Inf 0.030259 Inf Inputs: ---------- Min u Max __________ __________ __________ (1.) rct_airframe2/delta trim -Inf 0.00043574 Inf Outputs: None ----------
Nominal Controller Design
Design an H-infinity controller that responds to a step change in vertical acceleration under 1 second. Use a mixed-sensitivity formulation with a lowpass weight wS
that penalizes low-frequency tracking error and a highpass weight wT
that enforces adequate roll-off. Build an augmented plant with azref,delta
as inputs and the filtered wS*e,wT*az,e,q
as outputs. (For information about the mixed-sensitivity formulation, see Mixed-Sensitivity Loop Shaping.)
sb = sumblk('e = azref-az'); wS = makeweight(1e2,5,0.1); wS.u = 'e'; wS.y = 'we'; wT = makeweight(0.1,50,1e2); wT.u = 'az'; wT.y = 'waz'; Paug = connect(G,wS,wT,sb,{'azref','delta'},{'we','waz','e','q'});
Compute the optimal H-infinity controller using hinfsyn
. In this configuration there are two measurements e,q
and one control delta
.
[Knom,~,gam] = hinfsyn(Paug,2,1); gam
gam = 0.5928
Verify the open-loop gain against wS,wT
.
f = figure(); sigma(Knom*G,wS,'r--',1/wT,'g--'), grid legend('open-loop gain','> wS at low freq','< 1/wT at high freq')
ans = Legend (open-loop gain, > wS at low freq, < 1/wT at high freq) with properties: String: {'open-loop gain' '> wS at low freq' '< 1/wT at high freq'} Location: 'northeast' Orientation: 'vertical' FontSize: 9 Position: [0.5863 0.7614 0.2996 0.1144] Units: 'normalized' Use GET to show all properties
Plot the closed-loop response.
CL = feedback(G*Knom,diag([1 -1])); step(CL(:,1)), grid
Robustness Analysis
Compute the disk margins at the plant input and outputs. That the az
loop uses negative feedback while the q
loop uses positive feedback, so change the sign of the loop gain of the q
loop before using diskmargin
.
loopsgn = diag([1 -1]);
Examine the margins at the plant input.
DM = diskmargin(Knom*loopsgn*G)
DM = struct with fields: GainMargin: [0.3918 2.5524] PhaseMargin: [-47.2105 47.2105] DiskMargin: 0.8740 LowerBound: 0.8740 UpperBound: 0.8740 Frequency: 28.8842 WorstPerturbation: [1x1 ss]
For the margins at plant outputs, use the multiloop margin to account for simultaneous, independent variations in both output channels.
[~,MM] = diskmargin(loopsgn*G*Knom)
MM = struct with fields: GainMargin: [0.4994 2.0025] PhaseMargin: [-36.9261 36.9261] DiskMargin: 0.6678 LowerBound: 0.6678 UpperBound: 0.6691 Frequency: 23.6855 WorstPerturbation: [2x2 ss]
Finally, compute the margins against simultaneous variations at the plant input and outputs.
MMIO = diskmargin(loopsgn*G,Knom)
MMIO = struct with fields: GainMargin: [0.6866 1.4565] PhaseMargin: [-21.0565 21.0565] DiskMargin: 0.3717 LowerBound: 0.3717 UpperBound: 0.3725 Frequency: 24.8030 WorstPerturbation: [1x1 struct]
The disk-based gain and phase margins exceed 2 (6dB) and 35 degrees at the plant input and the plant outputs. Moreover, MMIO
indicates that this feedback loop can withstand gain variations by a factor 1.45 or phase variations of 21 degrees affecting all plant inputs and outputs at once.
Next, investigate the impact of gain and phase uncertainty on performance. Use the umargin
control design block to model a gain change factor of 1.4 (3dB) and phase change of 20 degrees in each feedback channel. Use getDGM
to fit an uncertainty disk to these amounts of gain and phase change.
GM = 1.4; PM = 20; DGM = getDGM(GM,PM,'balanced'); ue = umargin('e',DGM); uq = umargin('q',DGM); Gunc = blkdiag(ue,uq)*G; Gunc.OutputName = {'az','q'};
Rebuild the closed-loop model and sample the gain and phase uncertainty to gauge the impact on the step response.
CLunc = feedback(Gunc*Knom,loopsgn); step(CLunc(:,1),3) grid
The plot shows significant variability in performance, with large overshoot and steady-state error for some combinations of gain and phase variations.
Robust Controller Synthesis
Redo the controller synthesis to account for gain and phase variations using musyn
. This synthesis optimizes performance uniformly for the specified range of gain and phase uncertainty.
Punc = connect(Gunc,wS,wT,sb,{'azref','delta'},{'we','waz','e','q'}); [Kr,gam] = musyn(Punc,2,1); gam
D-K ITERATION SUMMARY: ----------------------------------------------------------------- Robust performance Fit order ----------------------------------------------------------------- Iter K Step Peak MU D Fit D 1 51.06 26.33 26.62 4 2 7.028 5.24 5.303 8 3 1.681 1.156 1.157 4 4 0.9705 0.9702 0.9822 10 5 0.962 0.9619 0.9622 10 6 0.9625 0.9623 0.9663 10 Best achieved robust performance: 0.962 gam = 0.9619
Compare the sampled step responses for the nominal and robust controllers. The robust design reduces both overshoot and steady-state errors and gives more consistent performance across the modeled range of gain and phase uncertainty.
CLr = feedback(Gunc*Kr,loopsgn);
rng(0) % for reproducibility
step(CLunc(:,1),3)
hold
rng(0)
step(CLr(:,1),3)
grid
Current plot held
The robust controller achieves this performance by increasing the (nominal) disk margins at the plant output and, to a lesser extent, the I/O disk margin. For instance, compare the disk-based margins at the plant outputs for the nominal and robust designs. Use diskmarginplot
to see the variations of the margins with frequency.
Lnom = loopsgn*G*Knom; Lrob = loopsgn*G*Kr; clf diskmarginplot(Lnom,Lrob) title('Disk margins at plant outputs') grid legend('Nominal','Robust')
ans = Legend (Nominal, Robust) with properties: String: {'Nominal' 'Robust'} Location: 'northeast' Orientation: 'vertical' FontSize: 8.1000 Position: [0.7756 0.8583 0.2018 0.0884] Units: 'normalized' Use GET to show all properties
Examine the margins against variations simultaneous variations at the plant inputs and outputs with the robust controller.
MMIO = diskmargin(loopsgn*G,Kr)
MMIO = struct with fields: GainMargin: [0.6492 1.5404] PhaseMargin: [-24.0167 24.0167] DiskMargin: 0.4254 LowerBound: 0.4254 UpperBound: 0.4263 Frequency: 19.6043 WorstPerturbation: [1x1 struct]
Recall that with the nominal controller, the feedback loop could withstand gain variations of a factor of 1.45 or phase variations of 21 degrees affecting all plant inputs and outputs at once. With the robust controller, those margins increase to about 1.54 and 24 degrees.
View the range of simultaneous gain and phase variations that the robust design can sustain at all plant inputs and outputs.
diskmarginplot(MMIO.GainMargin)
Nonlinear Simulation of Worst-Case Scenario
diskmargin
returns the smallest change in gain and phase that destabilizes the feedback loop. It can be insightful to inject this perturbation in the nonlinear simulation to further analyze the implications for the real system. For example, compute the destabilizing perturbation at the plant outputs for the nominal controller.
[~,MM] = diskmargin(Lnom);
WP = MM.WorstPerturbation;
bode(WP)
title('Smallest destabilizing perturbation')
The worst perturbation is a diagonal, dynamic perturbation that multiplies the open-loop response at the plant outputs. With this perturbation, the closed-loop system becomes unstable with an undamped pole at the frequency w0 = MM.Frequency
.
damp(feedback(WP*G*Knom,loopsgn))
Pole Damping Frequency Time Constant (rad/seconds) (seconds) -1.88e-03 1.00e+00 1.88e-03 5.33e+02 -2.55e-02 1.00e+00 2.55e-02 3.92e+01 -3.20e-02 1.00e+00 3.20e-02 3.12e+01 -5.16e-02 1.00e+00 5.16e-02 1.94e+01 -1.25e-01 1.00e+00 1.25e-01 7.98e+00 -3.85e+00 + 5.04e+00i 6.07e-01 6.34e+00 2.60e-01 -3.85e+00 - 5.04e+00i 6.07e-01 6.34e+00 2.60e-01 -8.38e+00 + 1.17e+01i 5.81e-01 1.44e+01 1.19e-01 -8.38e+00 - 1.17e+01i 5.81e-01 1.44e+01 1.19e-01 4.88e-14 + 2.37e+01i -2.06e-15 2.37e+01 -2.05e+13 4.88e-14 - 2.37e+01i -2.06e-15 2.37e+01 -2.05e+13 -2.95e+01 1.00e+00 2.95e+01 3.39e-02 -3.33e+01 1.00e+00 3.33e+01 3.00e-02 -6.04e+01 + 2.28e+01i 9.36e-01 6.46e+01 1.66e-02 -6.04e+01 - 2.28e+01i 9.36e-01 6.46e+01 1.66e-02 -3.86e+01 + 5.36e+01i 5.85e-01 6.61e+01 2.59e-02 -3.86e+01 - 5.36e+01i 5.85e-01 6.61e+01 2.59e-02 -1.10e+03 1.00e+00 1.10e+03 9.05e-04
w0 = MM.Frequency
w0 = 23.6855
Note that the gain and phase variations induced by this perturbation lie on the boundary of the range of safe gain/phase variations computed by diskmargin
. To confirm this result, compute the frequency response of the worst perturbation at the frequency w0
, convert it to a gain and phase variation, and visualize it along with the range of safe gain and phase variations represented by the multiloop disk margin.
DELTA = freqresp(WP,w0); clf diskmarginplot(MM.GainMargin) title('Range of stable gain and phase variations') hold on plot(20*log10(abs(DELTA(1,1))),abs( angle(DELTA(1,1))*180/pi),'ro') plot(20*log10(abs(DELTA(2,2))),abs( angle(DELTA(2,2))*180/pi),'ro')
To simulate the effect of this worst perturbation on the full nonlinear model in Simulink, insert it as a block before the controller block, as done in the modified model rct_airframeWP
.
close(f)
open_system('rct_airframeWP')
Here the MIMO Controller block is set to the nominal controller Knom
. To simulate the nonlinear response with this controller, compute the trim deflection and q
initial value from the operating condition op
.
delta_trim = op.Inputs.u + [1.5 0]*op.States(5).x; q_ini = op.States(4).x;
By commenting the WorstPerturbation block in and out, you can simulate the response with or without this perturbation. The results are shown below and confirm the destabilizing effect of the gain and phase variation computed by diskmargin
.
Figure 1: Nominal response.
Figure 2: Response with destabilizing gain/phase perturbation.
See Also
umargin
| diskmargin
| diskmarginplot
| musyn