MATLAB Examples

# MIMO Robustness Analysis

You can create and analyze uncertain state-space models made up of uncertain state-space matrices. In this example, create a MIMO system with parametric uncertainty and analyze it for robust stability and worst-case performance.

## Contents

Consider a two-input, two-output, two-state system whose model has parametric uncertainty in the state-space matrices. First create an uncertain parameter p. Using the parameter, make uncertain A and C matrices. The B matrix happens to be not-uncertain, although you will add frequency-domain input uncertainty to the model later.

```p = ureal('p',10,'Percentage',10); A = [0 p;-p 0]; B = eye(2); C = [1 p;-p 1]; H = ss(A,B,C,[0 0;0 0]) ```
```H = Uncertain continuous-time state-space model with 2 outputs, 2 inputs, 2 states. The model uncertainty consists of the following blocks: p: Uncertain real, nominal = 10, variability = [-10,10]%, 2 occurrences Type "H.NominalValue" to see the nominal value, "get(H)" to see all properties, and "H.Uncertainty" to interact with the uncertain elements. ```

You can view the properties of the uncertain system H using the get command.

```get(H) ```
``` NominalValue: [2x2 ss] Uncertainty: [1x1 struct] A: [2x2 umat] B: [2x2 double] C: [2x2 umat] D: [2x2 double] E: [] StateName: {2x1 cell} StateUnit: {2x1 cell} InternalDelay: [0x1 double] InputDelay: [2x1 double] OutputDelay: [2x1 double] Ts: 0 TimeUnit: 'seconds' InputName: {2x1 cell} InputUnit: {2x1 cell} InputGroup: [1x1 struct] OutputName: {2x1 cell} OutputUnit: {2x1 cell} OutputGroup: [1x1 struct] Notes: [0x1 string] UserData: [] Name: '' SamplingGrid: [1x1 struct] ```

Most properties behave in the same way as the corresponding properties of ss objects. The property NominalValue is itself an ss object.

## Adding Independent Input Uncertainty to Each Channel

The model for H does not include actuator dynamics. Said differently, the actuator models are unity-gain for all frequencies.

Nevertheless, the behavior of the actuator for channel 1 is modestly uncertain (say 10%) at low frequencies, and the high-frequency behavior beyond 20 rad/s is not accurately modeled. Similar statements hold for the actuator in channel 2, with larger modest uncertainty at low frequency (say 20%) but accuracy out to 45 rad/s.

Use ultidyn objects Delta1 and Delta2 along with shaping filters W1 and W2 to add this form of frequency domain uncertainty to the model.

```W1 = makeweight(.1,20,50); W2 = makeweight(.2,45,50); Delta1 = ultidyn('Delta1',[1 1]); Delta2 = ultidyn('Delta2',[1 1]); G = H*blkdiag(1+W1*Delta1,1+W2*Delta2) ```
```G = Uncertain continuous-time state-space model with 2 outputs, 2 inputs, 4 states. The model uncertainty consists of the following blocks: Delta1: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences Delta2: Uncertain 1x1 LTI, peak gain = 1, 1 occurrences p: Uncertain real, nominal = 10, variability = [-10,10]%, 2 occurrences Type "G.NominalValue" to see the nominal value, "get(G)" to see all properties, and "G.Uncertainty" to interact with the uncertain elements. ```

Note that G is a two-input, two-output uncertain system, with dependence on three uncertain elements, Delta1, Delta2, and p. It has four states, two from H and one each from the shaping filters W1 and W2, which are embedded in G.

You can plot a 2-second step response of several samples of G The 10% uncertainty in the natural frequency is obvious.

```stepplot(G,2) ```

You can create a Bode plot of samples of G. The high-frequency uncertainty in the model is also obvious. For clarity, start the Bode plot beyond the resonance.

```bodeplot(G,{13 100}) ```

## Closed-Loop Robustness Analysis

Load the controller and verify that it is two-input and two-output.

```load(fullfile(matlabroot,'examples','robust','mimoKexample.mat')) size(K) ```
```State-space model with 2 outputs, 2 inputs, and 9 states. ```

You can use the command loopsens to form all the standard plant/controller feedback configurations, including sensitivity and complementary sensitivity at both the input and output. Because G is uncertain, all the closed-loop systems are uncertain as well.

```F = loopsens(G,K) ```
```F = struct with fields: Si: [2x2 uss] Ti: [2x2 uss] Li: [2x2 uss] So: [2x2 uss] To: [2x2 uss] Lo: [2x2 uss] PSi: [2x2 uss] CSo: [2x2 uss] Poles: [13x1 double] Stable: 1 ```

F is a structure with many fields. The poles of the nominal closed-loop system are in F.Poles, and F.Stable is 1 if the nominal closed-loop system is stable. In the remaining 10 fields, S stands for sensitivity, T or complementary sensitivity, and L for open-loop gain. The suffixes i and o refer to the input and output of the plant. Finally, P and C refer to the plant and controller.

Hence, Ti is mathematically the same as:

Lo is G*K, and CSo is mathematically the same as

Examine the transmission of disturbances at the plant input to the plant output by plotting responses of F.PSi. Graph some samples along with the nominal.

```bodemag(F.PSi.NominalValue,'r+',F.PSi,'b-',{1e-1 100}) ```

## Nominal Stability Margins

You can use loopmargin to investigate loop-at-a-time gain and phase margins, loop-at-a-time disk margins, and simultaneous multivariable margins. They are computed for the nominal system and do not reflect the uncertainty models within G.

Explore the simultaneous margins individually at the plant input, individually at the plant output, and simultaneously at both input and output.

```[I,DI,SimI,O,DO,SimO,Sim] = loopmargin(G,K); ```

The third output argument is the simultaneous gain and phase variations allowed in all input channels to the plant.

```SimI ```
```SimI = struct with fields: GainMargin: [0.1187 8.4215] PhaseMargin: [-76.4565 76.4565] Frequency: 18.5851 ```

This information implies that the gain at the plant input can vary in both channels independently by factors between (approximately) 1/8 and 8, as well as phase variations up to 76 degrees.

The sixth output argument is the simultaneous gain and phase variations allowed in all output channels to the plant.

```SimO ```
```SimO = struct with fields: GainMargin: [0.1202 8.3178] PhaseMargin: [-76.2892 76.2892] Frequency: 17.8500 ```

Note that the simultaneous margins at the plant output are similar to those at the input. This is not always the case in multiloop feedback systems.

The last output argument is the simultaneous gain and phase variations allowed in all input and output channels to the plant. As expected, when you consider all such variations simultaneously, the margins are somewhat smaller than those at the input or output alone.

```Sim ```
```Sim = struct with fields: GainMargin: [0.5677 1.7616] PhaseMargin: [-30.8362 30.8362] Frequency: 9.1399 ```

Nevertheless, these numbers indicate a generally robust closed-loop system, able to tolerate significant gain (more than +/-50% in each channel) and 30 degree phase variations simultaneously in all input and output channels of the plant.

## Robust Stability Margin

With loopmargin, you determined various stability margins of the nominal multiloop system. These margins are computed only for the nominal system and do not reflect the uncertainty explicitly modeled by the ureal and ultidyn objects. When you work with a detailed uncertainty model, the conventional stability margins computed by loopmargin may not accurately reflect how close the system is from being unstable. You can then use robstab to compute the robust stability margin for the specified uncertainty.

In this example, use robstab to compute the robust stability margin for the uncertain feedback loop comprised of G and K. You can use any of the closed-loop transfer functions in F = loopsens(G,K). All of them, F.Si, F.To, etc., have the same internal dynamics, and hence their stability properties are the same.

```opt = robOptions('Display','on'); stabmarg = robstab(F.So,opt) ```
```Computing peak... Percent completed: 100/100 System is robustly stable for the modeled uncertainty. -- It can tolerate up to 221% of the modeled uncertainty. -- There is a destabilizing perturbation amounting to 222% of the modeled uncertainty. -- This perturbation causes an instability at the frequency 13.7 rad/seconds. stabmarg = struct with fields: LowerBound: 2.2129 UpperBound: 2.2173 CriticalFrequency: 13.7286 ```

This analysis confirms what the loopmargin analysis suggested. The closed-loop system is quite robust, in terms of stability, to the variations modeled by the uncertain parameters Delta1, Delta2, and p. In fact, the system can tolerate more than twice the modeled uncertainty without losing closed-loop stability.

## Worst-Case Gain Analysis

You can plot the Bode magnitude of the nominal output sensitivity function. It clearly shows decent disturbance rejection in all channels at low frequency.

```bodemag(F.So.NominalValue,{1e-1 100}) ```

You can compute the peak value of the maximum singular value of the frequency response matrix using norm.

```[PeakNom,freq] = getPeakGain(F.So.NominalValue) ```
```PeakNom = 1.1288 freq = 6.7969 ```

The peak is about 1.13. What is the maximum output sensitivity gain that is achieved when the uncertain elements Delta1, Delta2, and p vary over their ranges? You can use wcgain to answer this.

```[maxgain,wcu] = wcgain(F.So); maxgain ```
```maxgain = struct with fields: LowerBound: 2.1599 UpperBound: 2.1644 CriticalFrequency: 8.3356 ```

The analysis indicates that the worst-case gain is somewhere between 2.1 and 2.2. The frequency where the peak is achieved is about 8.5.

Use usubs to replace the values for Delta1, Delta2, and p that achieve the gain of 2.1. Make the substitution in the output complementary sensitivity, and do a step response.

```step(F.To.NominalValue,usubs(F.To,wcu),5) ```

The perturbed response, which is the worst combination of uncertain values in terms of output sensitivity amplification, does not show significant degradation of the command response. The settling time is increased by about 50%, from 2 to 4, and the off-diagonal coupling is increased by about a factor of about 2, but is still quite small.

You can also examine the worst-case frequency response alongside the nominal and sampled systems using wcsigma.

```wcsigma(F.To,{1e-1,100}) ```