Quantcast

Robust Control Toolbox

Robust Stability, Robust Performance and Mu Analysis

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.

System Description

Figure 1 shows the block diagram of a closed-loop system. The plant model $P$ is uncertain and the plant output $y$ must be regulated to remain small in the presence of disturbances $d$ and measurement noise $n$ .

Figure 1: Closed-loop system for robustness analysis

Disturbance rejection and noise insensitivity are quantified by the performance objective

$$\left\| \left( P(1+KP)^{-1}W_d  , (1+PK)^{-1}W_n \right) \right\|_{\infty}$$

where $W_d$ and $W_n$ are weighting functions reflecting the frequency content of $d$ and $n$ . Here $W_d$ is large at low frequencies and $W_n$ 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:

$$P(s) = \frac{16}{s^2+0.16 s+k} (1+W_u(s) \delta(s))$$

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);

Designing a Controller

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;

Closing the Loop

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.

Robust Stability Analysis

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

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');

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)

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 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');

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)
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

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.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.