Robust Stability and Mu Analysis
This demo shows how to use Robust Control Toolbox™ to analyze and quantify the robustness of feedback control systems, including modeling errors and parameter variations. We’ll look at how to test for robust stability with the robuststab function and gain insight into the connection with mu analysis and the mussv function.
Contents
System Description
Figure 1 shows the block diagram of a closed-loop system. The plant model P is uncertain, the signals d and n are disturbances, and the output e must be kept small for adequate disturbance rejection.
Figure 1: Closed-loop system for robustness analysis
Creating An Uncertain Plant Model
An uncertain plant model is a lightly-damped, second-order system with parametric uncertainty in the denominator coefficients and significant frequency-dependent unmodeled dynamics beyond 6 radians/second. The mathematical model looks like this:
The parameter k is assumed to be about 40% uncertain, with a nominal value of 16. The frequency-dependent uncertainty at the plant input is assumed to be about 20% at low frequency, rising to 100% at 6 radians/second. In this M-Code example, we create and combine uncertain elements to construct the uncertain plant model, P.
k = ureal('k',16,'Percentage',40); delta = ultidyn('delta',[1 1]); Wu = makeweight(0.2,6,10); P = tf(16,[1 0.16 k])*(1+Wu*delta);
Creating a Controller
For this example we use the following third-order controller:
num = [ -12.56 17.32 67.28]; den = [ 1 20.37 136.74 179.46]; C = tf(num,den);
Creating a Closed-Loop System with SYSIC
To build an uncertain model of the closed-loop system, we'll write the interconnection equations directly from the block diagram in Figure 1 and use the sysic function to compute the interconnection:
systemnames = 'P C'; inputvar = '[d; n]'; input_to_P = '[d + C]'; input_to_C = '[P + 0.2*n]'; outputvar = '[0.125*P]'; cleanupsysic = 'yes'; ClosedLoop = sysic;
Alternatively, we can use the icsignal and iconnect functions to build the closed-loop model:
IC = iconnect; d = icsignal(1); n = icsignal(1); y = icsignal(1); IC.Input = [d;n]; IC.Output = [0.125*y]; IC.Equation{1} = equate(y,P*(d+C*(0.2*n+y))); ClosedLoop = IC.System; ClosedLoop.InputName = {'d' 'n'}; ClosedLoop.OutputName = {'e'};
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, time-invariant dynamic element delta.
ClosedLoop
USS: 6 States, 1 Output, 2 Inputs, Continuous System delta: 1x1 LTI, max. gain = 1, 1 occurrence k: real, nominal = 16, variability = [-40 40]%, 1 occurrence
The nominal value of ClosedLoop is stable. The peak input/output gain of the nominal model is about 0.22:
pole(ClosedLoop.NominalValue) norm(ClosedLoop.NominalValue,inf)
ans = -11.0706 -1.8083 + 4.4576i -1.8083 - 4.4576i -4.1576 -1.6852 -60.9303 ans = 0.2168
Pick 10 random samples of the uncertain closed-loop systems
clf
step(usample(ClosedLoop,10),5)
title('Variability of closed-loop response due to model uncertainty')

Figure 2: Variability of closed-loop response due to model uncertainty
Robust Stability Analysis
Does the closed-loop system remain stable for all values of k, delta in the ranges specified above? We use robuststab to answer this basic robustness question:
[stabmarg,destabunc,report,info] = robuststab(ClosedLoop); stabmarg
stabmarg = UpperBound: 1.1437 LowerBound: 1.1436 DestabilizingFrequency: 5.8179
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.14, which means that the closed loop will remain stable for up to 114% 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 114% of the modeled uncertainty. -- A destabilizing combination of 114% of the modeled uncertainty exists, causing an instability at 5.82 rad/s. -- Sensitivity with respect to uncertain element ... 'delta' is 100%. Increasing 'delta' by 25% leads to a 25% decrease in the margin. 'k' is 84%. Increasing 'k' by 25% leads to a 21% decrease in the margin.
The variable destabunc contains the smallest combination of (k,|delta|) perturbations that causes instability.
destabunc
destabunc = delta: [1x1 ss] k: 23.3197
We can substitute these values into ClosedLoop, and verify that these values cause the closed-loop system to be unstable:
pole(usubs(ClosedLoop,destabunc))
ans = -59.7975 -48.5398 -10.4942 + 8.6663i -10.4942 - 8.6663i 0.0000 + 5.8179i 0.0000 - 5.8179i -2.4292 + 1.9138i -2.4292 - 1.9138i -1.7017
Note that the natural frequency of the unstable closed-loop pole is given by stabmarg.DestabilizingFrequency:
stabmarg.DestabilizingFrequency
ans = 5.8179
Robust Stability Analysis: Connection with Mu Analysis
The structured singular value, or mu, 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 closed-loop model ClosedLoop, where Delta is a block-diagonal matrix of (normalized) uncertain elements. The 3rd output argument of lftdata, BlkStruct, describes the block-diagonal 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));
Mu-analysis 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'); semilogx(mubnds) 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,pklidx] = max(mubnds(1,2).ResponseData(:)); [pku,pkuidx] = max(mubnds(1,1).ResponseData(:)); SM.UpperBound = 1/pkl; SM.LowerBound = 1/pku; SM.DestabilizingFrequency = omega(pklidx);
Compare SM to the robust stability margin bounds stabmarg computed with robuststab; they are identical:
stabmarg SM
stabmarg = UpperBound: 1.1437 LowerBound: 1.1436 DestabilizingFrequency: 5.8179 SM = UpperBound: 1.1437 LowerBound: 1.1436 DestabilizingFrequency: 5.8179
Finally, note that the mu bounds mubnd are available in the info output of robuststab:
semilogx(info.MussvBnds) 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)
Robust Performance Analysis
The input/output gain of a nominally-stable 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 elments from their nominal values and the worst-case system gain is shown in Figure 5. Robust performance analysis involves determining where the system performance degradation curve crosses the dashed-line (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 closed-loop 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 = UpperBound: 0.9891 LowerBound: 0.9886 CriticalFrequency: 5.8179
The report variable summarizes the robust performance analysis.
disp(report)
Uncertain System achieves a robust performance margin of 0.9891. -- A model uncertainty exists of size 98.9% resulting in a performance margin of 1.01 at 5.82 rad/sec. -- Sensitivity with respect to uncertain element ... 'delta' is 84%. Increasing 'delta' by 25% leads to a 21% decrease in the margin. 'k' is 62%. Increasing 'k' by 25% leads to a 16% 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.0133
1/perfmarg.UpperBound
ans = 1.0111
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.
semilogx(info.MussvBnds) xlabel('Frequency (rad/sec)'); ylabel('Mu upper/lower bounds'); title('Robust Performance Mu Plot');

Robust Performance: Connection with Mu Analysis
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 matrix-valued uncertain element (a ucomplexm object), with nominal value of 0, and norm-bounded 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 visa-versa, 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)
USS: 6 States, 0 Outputs, 0 Inputs, Continuous System PerfBlock: 2x1 complex matrix, 1 occurrence delta: 1x1 LTI, max. gain = 1, 1 occurrence k: real, nominal = 16, variability = [-40 40]%, 1 occurrence
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 block-diagonal structure of NDelta, and will be used as the block structure argument to mussv.
[M,NDelta,BlkStruct] = lftdata(modClosedLoop);
The mussv methods are frequency-based, 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,'sc'); semilogx(mubnds) xlabel('Frequency (rad/sec)'); ylabel('Mu upper/lower bounds'); title('Robust Performance Mu Plot');

The performance margin is the reciprical 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 upperbound of the performance margin.
[pkl,pklidx] = max(mubnds(1,2).ResponseData(:)); [pku,pkuidx] = max(mubnds(1,1).ResponseData(:));
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 omega(pklidx) ]
ans = 0.9886 0.9886 0.9891 0.9891 5.8179 5.8179
Worst-Case Gain Analysis
The closed-loop transfer function maps [d;n] to e. The worst-case 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.0886 UpperBound: 1.1273 CriticalFrequency: 5.8179
The variable maxgainunc contains the values of the uncertain elements associated with maxgain.LowerBound:
maxgainunc
maxgainunc = delta: [1x1 ss] k: 22.4000
We can verify that this perturbation causes the closed-loop system to have gain = maxgain.LowerBound:
maxgain.LowerBound norm(usubs(ClosedLoop,maxgainunc),inf)
ans = 1.0886 ans = 1.0913
Note that there is a slight difference in the answers. This is because wcgain uses a finite frequency grid to compute the worst-case gain, whereas norm uses an iterative refinement algorithm that estimates the peak gain more accurately.
Finally, compare the nominal and worst-case closed-loop gains:
bodemag(ClosedLoop.Nominal,'b-',usubs(ClosedLoop,maxgainunc),'r--',{.01,100}) legend('Nominal','Worst case')

Figure 6: Bode diagram comparing the nominal and worst-case closed-loop gains.
This analysis reveals that the nominal disturbance rejection properties of the controller C are not robust to the specified level of uncertainty.