## 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 closed-loop 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:** Closed-loop 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 lightly-damped, second-order system with parametric uncertainty in the denominator coefficients and significant frequency-dependent 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 frequency-dependent 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 Open-Loop 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 closed-loop 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, time-invariant dynamic element `delta`.

ClosedLoop

ClosedLoop = Uncertain continuous-time state-space 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 closed-loop 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 closed-loop 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 closed-loop 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 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'); LinMagopt = bodeoptions; LinMagopt.PhaseVisible = 'off'; LinMagopt.XLim = [1e-1 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 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 elements 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 = 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 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)
```

modClosedLoop = Uncertain continuous-time state-space 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 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); 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 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.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 closed-loop 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 worst-case gain, whereas `norm` uses 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','location','SouthWest')

**Figure 6:** Bode diagram comparing the nominal and worst-case closed-loop 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 worst-case gain calculations.