Robust Control Toolbox 

This example shows how to use Robust Control Toolbox™ to analyze and quantify the robustness of feedback control systems. It also provides insight into the connection with mu analysis and the mussv function.
On this page… 

Figure 1 shows the block diagram of a closedloop system. The plant model is uncertain and the plant output must be regulated to remain small in the presence of disturbances and measurement noise .
Figure 1: Closedloop system for robustness analysis
Disturbance rejection and noise insensitivity are quantified by the performance objective
where and are weighting functions reflecting the frequency content of and . Here is large at low frequencies and is large at high frequencies.
Wd = makeweight(100,.4,.15); Wn = makeweight(0.5,20,100); bodemag(Wd,'b',Wn,'g') title('Performance Weighting Functions') legend('Input disturbance','Measurement noise')
Creating an Uncertain Plant Model
The uncertain plant model P is a lightlydamped, secondorder system with parametric uncertainty in the denominator coefficients and significant frequencydependent unmodeled dynamics beyond 6 rad/s. The mathematical model looks like:
The parameter k is assumed to be about 40% uncertain, with a nominal value of 16. The frequencydependent uncertainty at the plant input is assumed to be about 30% at low frequency, rising to 100% at 10 rad/s, and larger beyond that. Construct the uncertain plant model P by creating and combining the uncertain elements:
k = ureal('k',16,'Percentage',30); delta = ultidyn('delta',[1 1],'SampleStateDim',4); Wu = makeweight(0.3,10,20); P = tf(16,[1 0.16 k]) * (1+Wu*delta);
We use the controller designed in the example "Improving Stability While Preserving OpenLoop Characteristics". The plant model used there happens to be the nominal value of the uncertain plant model created above. For completeness, we repeat the commands used to generate the controller.
K_PI = pid(1,0.8); K_rolloff = tf(1,[1/20 1]); Kprop = K_PI*K_rolloff; [negK,~,Gamma] = ncfsyn(P.NominalValue,Kprop); K = negK;
Use connect to build an uncertain model of the closedloop system of Figure 1. Name the signals coming in and out of each block and let connect do the wiring:
P.u = 'uP'; P.y = 'yP'; K.u = 'uK'; K.y = 'yK'; S1 = sumblk('uP = yK + D'); S2 = sumblk('uK = yP  N'); Wn.u = 'n'; Wn.y = 'N'; Wd.u = 'd'; Wd.y = 'D'; ClosedLoop = connect(P,K,S1,S2,Wn,Wd,{'d','n'},'yP');
The variable ClosedLoop is an uncertain system with two inputs and one output. It depends on two uncertain elements: a real parameter k and an uncertain linear, timeinvariant dynamic element delta.
ClosedLoop
ClosedLoop = Uncertain continuoustime statespace model with 1 outputs, 2 inputs, 10 states. The model uncertainty consists of the following blocks: delta: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences k: Uncertain real, nominal = 16, variability = [30,30]%, 1 occurrences Type "ClosedLoop.NominalValue" to see the nominal value, "get(ClosedLoop)" to see all properties, and "ClosedLoop.Uncertainty" to interact with the uncertain elements.
The classical margins from allmargin indicate good stability robustness to unstructured gain/phase variations within the loop.
allmargin(P.NominalValue*K)
ans = GainMargin: [6.3299 11.1423] GMFrequency: [1.6110 15.1667] PhaseMargin: [80.0276 99.6641 63.7989] PMFrequency: [0.4472 3.1460 5.2319] DelayMargin: [3.1236 1.4443 0.2128] DMFrequency: [0.4472 3.1460 5.2319] Stable: 1
Does the closedloop system remain stable for all values of k, delta in the ranges specified above? Answering this question requires a more sophisticated analysis using the robuststab function.
[stabmarg,destabunc,report,info] = robuststab(ClosedLoop); stabmarg
stabmarg = LowerBound: 1.4954 UpperBound: 1.4954 DestabilizingFrequency: 5.5678
The variable stabmarg gives upper and lower bounds on the robust stability margin, a measure of how much uncertainty on k, delta the feedback loop can tolerate before becoming unstable. For example, a margin of 0.8 indicates that as little as 80% of the specified uncertainty level can lead to instability. Here the margin is 1.5, which means that the closed loop will remain stable for up to 150% of the specified uncertainty. The report function summarizes these results:
report
report = Uncertain system is robustly stable to modeled uncertainty.  It can tolerate up to 150% of the modeled uncertainty.  A destabilizing combination of 150% of the modeled uncertainty was found.  This combination causes an instability at 5.57 rad/seconds.  Sensitivity with respect to the uncertain elements are: 'delta' is 100%. Increasing 'delta' by 25% leads to a 25% decrease in the margin. 'k' is 26%. Increasing 'k' by 25% leads to a 7% decrease in the margin.
The variable destabunc contains the combination of k and delta closest to their nominal values that causes instability.
destabunc
destabunc = delta: [1x1 ss] k: 23.1778
We can substitute these values into ClosedLoop and verify that these values cause the closedloop system to be unstable.
pole(usubs(ClosedLoop,destabunc))
ans = 1.0e+03 * 0.2094 + 0.0000i 0.0358 + 0.0000i 0.0131 + 0.0095i 0.0131  0.0095i 0.0000 + 0.0056i 0.0000  0.0056i 0.0030 + 0.0017i 0.0030  0.0017i 0.0009 + 0.0000i 0.0004 + 0.0000i 0.0200 + 0.0000i 2.3093 + 0.0000i 0.0000 + 0.0000i
Note that the natural frequency of the unstable closedloop pole is given by stabmarg.DestabilizingFrequency:
stabmarg.DestabilizingFrequency
ans = 5.5678
The structured singular value, or , is the mathematical tool used by robuststab to compute the robust stability margin. If you are comfortable with structured singular value analysis, you can use the mussv function directly to compute mu as a function of frequency and reproduce the results above. The function mussv is the underlying engine for all robustness analysis commands.
To use mussv, we first extract the (M,Delta) decomposition of the uncertain closedloop model ClosedLoop, where Delta is a blockdiagonal matrix of (normalized) uncertain elements. The 3rd output argument of lftdata, BlkStruct, describes the blockdiagonal structure of Delta and can be used directly by mussv
[M,Delta,BlkStruct] = lftdata(ClosedLoop);
For a robust stability analysis, only the channels of M associated with the uncertainty channels are used. Based on the row/column size of Delta, select the proper columns and rows of M. Remember that the rows of Delta correspond to the columns of M, and vice versa. Consequently, the column dimension of Delta is used to specify the rows of M:
szDelta = size(Delta); M11 = M(1:szDelta(2),1:szDelta(1));
Muanalysis is performed on a finite grid of frequencies. For comparison purposes, evaluate the frequency response of M11 over the same frequency grid as used for the robuststab analysis.
omega = info.Frequency; M11_g = frd(M11,omega);
Compute mu(M11) at these frequencies and plot the resulting lower and upper bounds:
mubnds = mussv(M11_g,BlkStruct,'s'); LinMagopt = bodeoptions; LinMagopt.PhaseVisible = 'off'; LinMagopt.XLim = [1e1 1e2]; LinMagopt.MagUnits = 'abs'; bodeplot(mubnds(1,1),mubnds(1,2),LinMagopt); xlabel('Frequency (rad/sec)'); ylabel('Mu upper/lower bounds'); title('Mu plot of robust stability margins (inverted scale)');
Figure 3: Mu plot of robust stability margins (inverted scale)
The robust stability margin is the reciprocal of the structured singular value. Therefore upper bounds from mussv become lower bounds on the stability margin. Make these conversions and find the destabilizing frequency where the mu upper bound peaks (that is, where the stability margin is smallest):
[pkl,wPeakLow] = norm(mubnds(1,2),inf); [pku] = norm(mubnds(1,1),inf); SMfromMU.LowerBound = 1/pku; SMfromMU.UpperBound = 1/pkl; SMfromMU.DestabilizingFrequency = wPeakLow;
Compare SMfromMU to the robust stability margin bounds stabmarg computed with robuststab. They are identical:
stabmarg SMfromMU
stabmarg = LowerBound: 1.4954 UpperBound: 1.4954 DestabilizingFrequency: 5.5678 SMfromMU = LowerBound: 1.4954 UpperBound: 1.4954 DestabilizingFrequency: 5.5678
Finally, note that the same mu bounds mubnd are available in the info output of robuststab:
bodeplot(info.MussvBnds(1,1),info.MussvBnds(1,2),LinMagopt) xlabel('Frequency (rad/sec)'); ylabel('Mu upper/lower bounds'); title('Mu plot of robust stability margins (inverted scale)');
Figure 4: Mu plot of robust stability margins (inverted scale)
The input/output gain of a nominallystable uncertain system model will generally degrade for specific values of its uncertain elements. Moreover, the maximum possible degradation increases as the uncertain elements are allowed to further deviate from their nominal values.
A typical tradeoff curve between allowable deviation of uncertain elements from their nominal values and the worstcase system gain is shown in Figure 5. Robust performance analysis involves determining where the system performance degradation curve crosses the dashedline (the y=1/x hyperbola).
Figure 5: Generic tradeoff between uncertainty level and performance.
The simplest route to analyzing the robust performance margin of the closedloop system is direct use of the robustperf command.
[perfmarg,perfmargunc,report,info] = robustperf(ClosedLoop);
The perfmarg variable has upper and lower bounds on the performance margin.
perfmarg
perfmarg = LowerBound: 0.8502 UpperBound: 0.8505 CriticalFrequency: 6.1898
The report variable summarizes the robust performance analysis.
disp(report)
Uncertain system does not achieve performance robustness to modeled uncertainty.  The tradeoff of model uncertainty and system gain is balanced at a level of 85% of the modeled uncertainty.  A model uncertainty of 85% can lead to input/output gain of 1.18 at 6.19 rad/seconds.  Sensitivity with respect to the uncertain elements are: 'delta' is 49%. Increasing 'delta' by 25% leads to a 12% decrease in the margin. 'k' is 38%. Increasing 'k' by 25% leads to a 10% decrease in the margin.
perfmargunc is a structure of values of uncertain elements associated with the hyperbola crossing. We can substitute the values into ClosedLoop, and verify that this collection of values causes the closed loop system norm to be greater than or equal to the reciprocal of the performance margin upper bound.
norm(usubs(ClosedLoop,perfmargunc),inf)
ans = 1.1848
1/perfmarg.UpperBound
ans = 1.1758
Finally, we plot the bounds from mussv, which is the underlying engine for the robustness analysis. The peak value is the reciprocal of the performance margin, and the frequency at which the peak occurs is the critical frequency.
bodeplot(info.MussvBnds(1,1),info.MussvBnds(1,2),LinMagopt) xlabel('Frequency (rad/sec)'); ylabel('Mu upper/lower bounds'); title('Robust Performance Mu Plot');
If you are comfortable with structured singular value analysis, the method described above may seem too transparent, and leave you wondering what calculations actually took place. Instead, we can choose to extract the (M,Delta) decomposition, and call mussv on the appropriate channels of M.
Robust performance analysis requires a fictitious uncertain element, often referred to as the "performance block." It should be a complex matrixvalued uncertain element (a ucomplexm object), with nominal value of 0, and normbounded by 1). This element is "wrapped" around the input/output channels of the system under consideration, and then a robust stability analysis is performed. Since the rows of uncertain elements correspond to the columns of M, and visaversa, the row dimension of the performance block should be the column (input) dimension of ClosedLoop.
Create PerfBlock, and close the input/output channels of ClosedLoop with PerfBlock, yielding modClosedLoop.
PerfBlock = ucomplexm('PerfBlock',zeros(size(ClosedLoop,2),size(ClosedLoop,1)));
modClosedLoop = lft(PerfBlock,ClosedLoop)
modClosedLoop = Uncertain continuoustime statespace model with 0 outputs, 0 inputs, 10 states. The model uncertainty consists of the following blocks: PerfBlock: Uncertain 2x1 complex matrix, 1 occurrences delta: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences k: Uncertain real, nominal = 16, variability = [30,30]%, 1 occurrences Type "modClosedLoop.NominalValue" to see the nominal value, "get(modClosedLoop)" to see all properties, and "modClosedLoop.Uncertainty" to interact with the uncertain elements.
Note that modClosedLoop has 0 inputs and 0 outputs (this is expected), but has dependence on the original uncertain elements, as well as the new performance block.
We can use the lftdata command to separate the uncertain system, modClosedLoop into a certain system M, in feedback with a normalized, block diagonal uncertain matrix NDelta. The 3rd output argument, BlkStruct concisely describes the blockdiagonal structure of NDelta, and will be used as the block structure argument to mussv.
[M,NDelta,BlkStruct] = lftdata(modClosedLoop);
The mussv methods are frequencybased, so we first compute a frequency response of M. For comparison purposes, use the frequency vector that was used in the robustperf analysis.
omega = info.Frequency; M_g = frd(M,omega);
Generate and plot mussv bounds for M. Note that the plot is identical to the plot obtained from the robustperf analysis.
[mubnds] = mussv(M_g,BlkStruct); bodeplot(mubnds(1,1),mubnds(1,2),LinMagopt) xlabel('Frequency (rad/sec)'); ylabel('Mu upper/lower bounds'); title('Robust Performance Mu Plot');
points completed (of 147) ... 147
The performance margin is the reciprocal of the structured singular value. Therefore upper bounds from mussv become lower bounds on the performance margin. Make these conversions and find the frequency associated with the upper bound of the performance margin.
[pkl,wPeakLow] = norm(mubnds(1,2),inf); [pku] = norm(mubnds(1,1),inf);
To verify that the two calculations are identical, we compare the lower and upper bounds from robustperf and mussv as well as the corresponding critical frequencies:
[ perfmarg.LowerBound 1/pku perfmarg.UpperBound 1/pkl perfmarg.CriticalFrequency wPeakLow ]
ans = 0.8502 0.8502 0.8505 0.8505 6.1898 6.1898
The closedloop transfer function maps [d;n] to e. The worstcase gain of this transfer function indicates where disturbance rejection is worst. You can use wcgain to assess the worst (largest) value of this gain:
[maxgain,maxgainunc,info] = wcgain(ClosedLoop); maxgain
maxgain = LowerBound: 1.5361 UpperBound: 1.5364 CriticalFrequency: 6.1898
The variable maxgainunc contains the values of the uncertain elements associated with maxgain.LowerBound:
maxgainunc
maxgainunc = delta: [1x1 ss] k: 20.8000
We can verify that this perturbation causes the closedloop system to have gain at least equal to the value of maxgain.LowerBound:
maxgain.LowerBound norm(usubs(ClosedLoop,maxgainunc),inf)
ans = 1.5361 ans = 1.5347
Note that there is a difference in the answers, and the gain is actually slightly larger than the lower bound. This is because wcgain uses a finite frequency grid to compute the worstcase gain, whereas norm uses estimates the peak gain more accurately.
Finally, compare the nominal and worstcase closedloop gains:
bodemag(ClosedLoop.Nominal,'b',usubs(ClosedLoop,maxgainunc),'r',{.01,100}) legend('Nominal','Worst case','location','SouthWest')
Figure 6: Bode diagram comparing the nominal and worstcase closedloop gains.
This analysis reveals that the nominal disturbance rejection and noise insensitivity properties of the controller K are not robust to the specified level of uncertainty. The degree to which performance is not robust to the uncertainty level is quantified in the robust performance and worstcase gain calculations.