# Uncertainty modeling

We provide two general means of modeling uncertainty, the traditional $M\Delta$ framework ^{[Skogestad]}^{[Doyle91]}, and using parametric uncertainty. Support for parametric uncertainty is almost universal in Julia, not only in ControlSystems.jl, by means of computing with uncertain number types. In this tutorial, we will use a Monte-Carlo approach using uncertain number types from MonteCarloMeasurements.jl.

Both the $M\Delta$ framework and parametric-uncertainty approaches are illustrated below.

- Uncertainty modeling

## Uncertainty API

`δc`

Creates an uncertain complex parameter.`δr`

Creates an uncertain real parameter.`δss`

(Experimental) Creates an uncertain statespace model.`neglected_delay`

Create a multiplicative weight that represents uncertainty from an unmodeled delay.`neglected_lag`

Create a multiplicative weight that represents uncertainty from an unmodeled lag (pole).`gain_and_delay_uncertainty`

Create a multiplicative weight that represents uncertainty from uncertain gains and delay.`makeweight`

Create a custom weighting function.`fit_complex_perturbations`

- See MonteCarloMeasurements.jl to create uncertain parameters that are represented by samples.
`sys_from_particles`

Convert from an uncertain representation using`Particles`

to a "multi-model" representation using multiple`StateSpace`

models.`ss2particles`

Convert a vector of state-space models to a single state-space model with coefficient type`MonteCarloMeasurements.Particles`

.

See example `uncertain.jl`

.

## Parametric uncertainty using MonteCarloMeasurements.jl

The most straightforward way to model uncertainty is to use uncertain parameters, using tools such as IntervalArithmetic (strict, worst case guarantees) or MonteCarloMeasurements (less strict worst-case analysis or probabilistic). In the following, we show an example with MIMO systems with both parametric uncertainty and diagonal, complex uncertainty, adapted from 8.11.3 in Skogestad, "Multivariable Feedback Control: Analysis and Design". This example is also available as a julia script in `uncertain.jl`

.

We will create uncertain parameters using the `δr`

constructor from this package. One may alternatively create uncertain parameters directly using any of the constructors from MonteCarloMeasurements.jl. Most functions from ControlSystemsBase.jl should work with systems containing parameters from MonteCarloMeasurements.jl.

### Basic example

This example shows how to use MonteCarloMeasurements directly to build uncertain systems.

```
using ControlSystemsBase, MonteCarloMeasurements, Plots
ω = 1 ± 0.1 # Create an uncertain Gaussian parameter
```

```
1.0 ± 0.1 MonteCarloMeasurements.Particles{Float64, 2000}
```

`ζ = 0.3..0.4 # Create an uncertain uniform parameter`

```
0.35 ± 0.0289 MonteCarloMeasurements.Particles{Float64, 2000}
```

`G = tf(ω^2, [1, 2ζ*ω, ω^2]) # systems accept uncertain parameters`

```
TransferFunction{Continuous, ControlSystemsBase.SisoRational{MonteCarloMeasurements.Particles{Float64, 2000}}}
1.01 ± 0.2
---------------------------------
1.0s^2 + 0.7 ± 0.09s + 1.01 ± 0.2
Continuous-time transfer function model
```

```
w = exp10.(-2:0.02:2)
bodeplot(G, w)
```

`plot(step(G, 0:0.1:20))`

### Example: Spinning satellite

This example makes use of real-valued uncertain parameters created using `δr`

, it comes from section 3.7.1 of Skogestad's book.

```
using RobustAndOptimalControl, ControlSystemsBase, MonteCarloMeasurements, Plots, LinearAlgebra
unsafe_comparisons(true)
a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
K = ss(1.0I(2))
w = 2π .* exp10.(LinRange(-2, 2, 500))
S, PS, CS, T = gangoffour(P, K)
sigmaplot(S, w, lab="S")
sigmaplot!(T, w, c=2, lab="T", ylims=(0.01, 45))
```

Both sensitivity functions are very large, expect a non-robust system!

Next, we add parametric uncertainty

```
a = 10*(1 + 0.1δr(100)) # Create an uncertain parameter with nominal value 10 and 10% uncertainty, represented by 100 samples
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
Sp, PSp, CSp, Tp = gangoffour(P, K)
sigmaplot(Sp, w, lab="S")
sigmaplot!(Tp, w, c=2, lab="T", ylims=(0.01, 100))
```

Not only are sensitivity functions large, they vary a lot under the considered uncertainty. We can also plot a step response of one of the sensitivity functions to check how the system behaves

`plot(step(c2d(Tp, 0.01), 10))`

This kind of plot is quite useful, it immediately tells you that this transfer function appears stable, and that there is uncertainty in the static gain etc.

Next, we add complex diagonal multiplicative input uncertainty. With input uncertainty of magnitude $ϵ < \dfrac{1}{σ̄(T)}$ we are guaranteed robust stability (even for “full-block complex perturbations")

```
a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
W0 = makeweight(0.2, (20,1), 2)
W = I(2) + W0 .* diagm([δc(100), δc(100)]) # Create a diagonal complex uncertainty weighted in frequency by W0, use 100 samples
Ps = P*W
Ss, PSs, CSs, Ts = gangoffour(Ps, K)
sigmaplot(Ss, w, lab="S")
sigmaplot!(Ts, w, c=2, lab="T", ylims=(0.01, 100))
```

Under this uncertainty, the sensitivity could potentially be sky high., note how some of the 100 realizations peak much higher than the others. This is an indication that the system might be unstable.

With complex entries in the system model, we can't really plot the step response, but we can plot, e.g., the absolute value

```
res = step(c2d(Ts, 0.01), 10)
plot(res.t, [abs.(res.y)[1,:,1] abs.(res.y)[2,:,2]]) # plot only the diagonal response
```

Looks unstable to me. The analysis using $M\Delta$ methodology below will also reach this conclusion.

### Example: Distillation Process

This example comes from section 3.7.2 of Skogestad's book. In this example, we'll explore also complex uncertainties, created using `δc`

.

```
using RobustAndOptimalControl, ControlSystemsBase, MonteCarloMeasurements, Plots, LinearAlgebra
unsafe_comparisons(true)
M = [87.8 -86.4; 108.2 -109.6]
G = Ref(ss(tf(1, [75, 1]))) .* M
RGA = relative_gain_array(G, 0)
sum(abs, RGA) # A good estimate of the true condition number, which is 141.7
```

`138.2752186588932`

large elements in the RGA indicate a process that is difficult to control

We consider the following inverse-based controller, which may also be looked upon as a steady-state decoupler with a PI controller

```
k1 = 0.7
Kinv = Ref(ss(tf(k1*[75, 1], [1, 0]))) .* inv(M)
# reference filter
F = tf(1, [5, 1]) .* I(2)
w = 2π .* exp10.(LinRange(-2, 2, 500))
sigmaplot(input_sensitivity(G, Kinv), w)
sigmaplot!(output_sensitivity(G, Kinv), w, c=2)
```

Sensitivity looks nice, how about step response

`plot(step(feedback(G*Kinv)*F, 20))`

Looks excellent..

We consider again the input gain uncertainty as in the previous example, and we manually select the perturbations to be $ϵ_1 = 0.2$ and $ϵ_2 = 0.2$. We then have

```
G′ = G * diagm([1 + 0.2, 1 - 0.2])
plot!(step(feedback(G′*Kinv)*F, 20), l=:dash)
```

Looks very poor! The system was not robust to simultaneous input uncertainty!

We can also do this with a real, diagonal input uncertainty that grows with frequency

```
W0 = makeweight(0.2, 1, 2.0) # uncertainty goes from 20% at low frequencies to 200% at high frequencies
W = I(2) + W0 .* diagm([δr(100), δr(100)])
Gs = G*W
plot(step(feedback(G*Kinv)*F, 20))
plot!(step(feedback(G′*Kinv)*F, 20), l=:dash)
res = step(c2d(feedback(Gs*Kinv)*F, 0.01), 20)
mcplot!(res.t, abs.(res.y[:, :, 1]'), alpha=0.3)
mcplot!(res.t, abs.(res.y[:, :, 2]'), alpha=0.3)
```

The system is very sensitive to real input uncertainty!

With a complex, diagonal uncertainty, modeling both gain and phase variations, it looks slightly worse, but not much worse than with real uncertainty.

```
W = I(2) + W0 .* diagm([δc(100), δc(100)]) # note δc instead of δr above
Gs = G*W
res = step(c2d(feedback(Gs*Kinv)*F, 0.01), 20)
mcplot!(res.t, abs.(res.y[:, :, 1]'), alpha=0.3)
mcplot!(res.t, abs.(res.y[:, :, 2]'), alpha=0.3)
```

How about the sensitivity functions?

```
Si = input_sensitivity(Gs, Kinv)
sigmaplot(Si, w, c=1, lab="Si")
So = output_sensitivity(Gs, Kinv)
sigmaplot!(So, w, c=2, lab="So")
```

The sensitivity at the plant output is enormous. A low sensitivity with the nominal system does not guarantee robustness!

## Model-order reduction for uncertain models

The $\nu$-gap metric is a measure of distance between models when they are used in a feedback loop. This metric has the nice property that a controller designed for a process $P$ that achieves a normalized coprime factor margin (`ncfmargin`

) of $m$, will stabilize all models that are within a $\nu$-gap distance of $m$ from $P$. This can be used to reduce the number of uncertain realizations for a model represented with `Particles`

like above in a smart way. Say that we have a plant model $P$

```
using RobustAndOptimalControl, ControlSystemsBase, MonteCarloMeasurements
ω = with_nominal(0.9 .. 1.1, 1)
ζ = with_nominal(0.5 ± 0.01, 0.5)
P = tf([ω^2], [1, 2*ζ*ω, ω^2]) |> ss
```

```
StateSpace{Continuous, MonteCarloMeasurements.Particles{Float64, 2000}}
A =
0.0 1.0
-1.0 ± 0.12 -1.0 ± 0.061
B =
0.0
1.0
C =
1.0 ± 0.12 0.0
D =
0.0
Continuous-time state-space model
```

represented by 2000 samples (indicated by the displayed type `Particles{Float64, 2000}`

). If we plot $P$, it looks something like this:

```
w = exp10.(LinRange(-1, 0.5, 150))
# nyquistplot(P, w, lab="Original P", xlims=(-1.1,1.1), ylims=(-1.5,0.7), points=true, format=:png, dpi=80)
bodeplot(P, w, lab="Original P", plotphase = false, format=:png, dpi=80, ri=false, c=1, legend=true)
```

We can compute the $\nu$-gap metric between each realization in $P$ and the nominal value (encoded using `with_nominal`

above):

`gap = nugap(P)`

```
0.0606433 ± 0.0335 MonteCarloMeasurements.Particles{Float64, 2000}
```

The worst-case gap is:

`pmaximum(gap) # p for "particle" maximum`

`0.12986707291999236`

That means that if we design a controller for the nominal $P$ without any uncertainty, and make sure that it achieves an `ncfmargin`

of at least this value, it will stabilize all realizations in $P$.

We can also reduce the number of realizations in $P$ by discarding those that are close in the $\nu$-gap sense to the nominal value:

`Pr = nu_reduction(P, 0.1)`

```
StateSpace{Continuous, MonteCarloMeasurements.Particles{Float64, 316}}
A =
0.0 1.0
-0.972 ± 0.18 -0.98 ± 0.088
B =
0.0
1.0
C =
0.972 ± 0.18 0.0
D =
0.0
Continuous-time state-space model
```

here, all realizations that were within a $\nu$-gap distance of 0.1 from the nominal value were discarded. `nu_reduction`

usually reduces the number of realizations substantially. The plot of $P_r$ looks like

```
# nyquistplot(Pr, lab="Reduced P", xlims=(-1.1,1.1), ylims=(-1.5,0.7), points=true, format=:png, dpi=80)
bodeplot!(Pr, w, lab="Reduced P", plotphase = false, format=:png, dpi=80, ri=false, c=2, l=2)
```

we see that the reduction kept the realizations that were furthest away from the nominal value.

We can reduce the number of realizations even further using `nu_reduction_recursive`

:

`Prr = nu_reduction_recursive(P, 0.1)`

```
StateSpace{Continuous, MonteCarloMeasurements.Particles{Float64, 3}}
A =
0.0 1.0
-1.0 ± 0.2 -0.979 ± 0.055
B =
0.0
1.0
C =
1.0 ± 0.2 0.0
D =
0.0
Continuous-time state-space model
```

```
# nyquistplot(Prr, lab="Recursively reduced P", xlims=(-1.1,1.1), ylims=(-1.5,0.7), points=true, format=:png, dpi=80)
bodeplot!(Prr, w, lab="Recursively reduced P", plotphase = false, format=:png, dpi=80, ri=false, c=3, l=3)
```

We now have only three realizations left, the nominal one and the two extreme cases (in the $\nu$-gap sense).

The algorithm used in `nu_reduction_recursive`

has a worst-case complexity of $O(N^2)$ where $N$ is the number of realizations (particles) in $P$, but this complexity is only problematic for small gaps and large number of realizations, say, $\nu < 0.02$ and $N > 50$.

With the reduced model, we can more easily perform loop-shaping or other control design tasks, as long as we keep track of `ncfmargin`

staying above our $\nu$-gap.

If `P`

has a stochastic interpretation, i.e., the coefficients come from some distribution, this interpretation will be lost after reduction, mean values and standard deviations will not be preserved. The reduced system should instead be interpreted as preserving worst-case uncertainty.

For more background on the $\nu$-gap metric, see When are two systems similar? and the book by Skogestad and Postlethwaite or by Åström and Murray.

## Using the $M\Delta$ framework

The examples above never bothered with things like the "structured singular value", $\mu$ or linear-fractional transforms. We do, however, provide some elementary support for this modeling framework.

In robust control, we often find ourselves having to consider the feedback interconnections below.

```
┌─────────┐
zΔ◄───┤ │◄────wΔ
│ │
z◄───┤ P │◄────w
│ │
y◄───┤ │◄────u
└─────────┘
```

```
┌─────────┐
zΔ◄───┤ │◄────wΔ
│ │
z◄───┤ P │◄────w
│ │
y┌───┤ │◄───┐u
│ └─────────┘ │
│ ┌───┐ │
└─────►│ K ├───────┘
└───┘
```

```
┌───┐
┌─────►│ Δ ├───────┐
│ └───┘ │
│ ┌─────────┐ │
zΔ└───┤ │◄───┘wΔ
│ │
z◄───┤ P │◄────w
│ │
y┌───┤ │◄───┐u
│ └─────────┘ │
│ ┌───┐ │
└─────►│ K ├───────┘
└───┘
```

The first block diagram denotes an open-loop system $P$ with an uncertainty mapping $w_\Delta = \Delta z_\Delta$, a *performance mapping* from $w$ to $z$ and a input-output mapping between $u$ and $y$. Such a system $P$ can be partitioned as

\[P = \begin{bmatrix} P_{11} & P_{12} & P_{13}\\ P_{21} & P_{22} & P_{23}\\ P_{31} & P_{32} & P_{33}\\ \end{bmatrix}\]

where each $P(s)_{ij}$ is a transfer matrix. The type `UncertainSS`

with constructor `uss`

represents the block

\[P = \begin{bmatrix} P_{11} & P_{12}\\ P_{21} & P_{22}\\ \end{bmatrix}\]

while an `ExtendedStateSpace`

object represents the block

\[P = \begin{bmatrix} P_{22} & P_{23}\\ P_{32} & P_{33}\\ \end{bmatrix}\]

there is thus no type that represents the full system $P$ above. However, we provide the function `partition`

which allows you to convert from a regular statespace system to an extended statespace object, and it is thus possible to represent $P$ by placing the whole block

\[P = \begin{bmatrix} P_{22} & P_{23}\\ P_{32} & P_{33}\\ \end{bmatrix}\]

into $P_{22}$ for the purposes of uncertainty analysis (use `ss`

to convert it to a standard statespace object), and later use `partition`

to recover the internal block structure.

Given an `UncertainSS`

$P$, we can close the loop around $\Delta$ by calling `lft(P, Δ, :u)`

, and given an `ExtendedStateSpace`

, we can close the loop around `K`

by calling `starprod(P, K)`

or `lft(P, K)`

(using positive feedback). This works even if `P`

is a regular statespace object, in which case the convention is that the inputs and outputs are ordered as in the block diagrams above. The number of signals that will be connected by `lft`

is determined by the input-output arity of $K$ and $\Delta$ respectively.

We have the following methods for `lft`

(in addition to the standard ones in ControlSystemsBase.jl)

`lft(G::UncertainSS, K::LTISystem)`

forms the lower LFT closing the loop around $K$.`lft(G::UncertainSS, Δ::AbstractArray=G.Δ)`

forms the upper LFT closing the loop around $\Delta$.`lft(G::ExtendedStateSpace, K)`

forms the lower LFT closing the loop around $K$.

### Robust stability and performance

To check robust stability of the system in the last block diagram (with or without $z$ and $w$), we can use the functions `structured_singular_value`

, `robstab`

and `diskmargin`

.

Currently, `structured_singular_value`

is rather limited and supports diagonal complex blocks only. If $\Delta$ is a single full complex block, `opnorm(P.M) < 1`

is the condition for stability.

Robust performance can be verified by introducing an additional fictitious "performance perturbation" $\Delta_p$ which is a full complex block, around which we close the loop from $z$ to $w$ and check the `structured_singular_value`

with the augmented perturbation block

\[\Delta_a = \begin{bmatrix} \Delta & 0\\ 0 & \Delta_p\\ \end{bmatrix}\]

### Examples

We repeat the first example here, but using $M\Delta$ formalism rather than direct Monte-Carlo modeling.

When we call `δc`

without any arguments, we get a symbolic (or structured) representation of the uncertainty rather than the sampled representation we got from calling `δc(100)`

.

```
a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
W0 = makeweight(0.2, (1,1), 2) |> ss
W = I(2) + W0 .* uss([δc(), δc()]) # Create a diagonal complex uncertainty weighted in frequency by W0
Ps = P*W
```

```
UncertainSS{Continuous}
A =
0.0 10.0 -1.590990257669732 0.0
-10.0 0.0 0.0 -1.590990257669732
0.0 0.0 -1.7677669529663689 0.0
0.0 0.0 0.0 -1.7677669529663689
B =
2.0 0.0 1.0 0.0
0.0 2.0 0.0 1.0
2.0 0.0 0.0 0.0
0.0 2.0 0.0 0.0
C =
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
1.0 10.0 0.0 0.0
-10.0 1.0 0.0 0.0
D =
0.0 0.0 1.0 0.0
0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
Continuous-time state-space model
```

`Ps`

is now represented as a upper linear fractional transform (upper LFT).

We can draw samples from this uncertainty representation (sampling of $\Delta$ and closing the loop `starprod(Δ, Ps)`

) like so

```
Psamples = rand(Ps, 100)
sigmaplot(Psamples, w)
```

We can extract the nominal model using

`system_mapping(Ps)`

```
StateSpace{Continuous, Float64}
A =
0.0 10.0
-10.0 0.0
B =
1.0 0.0
0.0 1.0
C =
1.0 10.0
-10.0 1.0
D =
0.0 0.0
0.0 0.0
Continuous-time state-space model
```

And obtain $M$ and $\Delta$ when the loop is closed with $K$ like this:

`lft(Ps, K).M`

```
StateSpace{Continuous, Float64}
A =
1.0 20.0 -1.590990257669732 0.0
-20.0 1.0 0.0 -1.590990257669732
0.0 0.0 -1.7677669529663689 0.0
0.0 0.0 0.0 -1.7677669529663689
B =
2.0 0.0
0.0 2.0
2.0 0.0
0.0 2.0
C =
1.0 10.0 0.0 0.0
-10.0 1.0 0.0 0.0
D =
0.0 0.0
0.0 0.0
Continuous-time state-space model
```

`Ps.Δ # Ps.delta also works`

```
2-element Vector{Any}:
δ{ComplexF64, Float64}(0.0 + 0.0im, 1.0, Symbol("##δ#561"))
δ{ComplexF64, Float64}(0.0 + 0.0im, 1.0, Symbol("##δ#562"))
```

We can evaluate the frequency response of $M$ and calculate the structured singular value $\mu$

```
M = freqresp(lft(Ps, -K).M, w) # -K to get negative feedback
μ = structured_singular_value(M)
plot(w, μ, xscale=:log10, title="Structured singular value μ", xlabel="Frequency [rad/s]", ylabel="μ")
```

$\mu$ is very high, whenever $\mu > 1$, the system is not stable with respect to the modeled uncertainty. The tolerated uncertainty is only about $\dfrac{1}{||\mu||_\infty}$

`1/norm(μ, Inf)`

`0.13647901995727751`

of the modeled uncertainty. Another way of calculating this value is

`robstab(lft(Ps, -K))`

`0.136477575525228`

### Internals of the $M\Delta$ framework

TODO

## Uncertain time delays

Modeling uncertain time delays can be done in several ways, one approach is to make use of a multiplicative uncertainty weight created using `neglected_delay`

multiplied by an uncertain element created using `δc`

, example:

```
using RobustAndOptimalControl, ControlSystemsBase, MonteCarloMeasurements, Plots, LinearAlgebra
a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) # Plant
W0 = neglected_delay(0.005) |> ss # Weight
W = I(2) + W0 .* uss([δc(), δc()]) # Create a diagonal real uncertainty weighted in frequency by W0
Ps = P*W # Uncertain plant
Psamples = rand(Ps, 1000) # Sample the uncertain plant for plotting
w = exp10.(LinRange(-1, 3, 300)) # Frequency vector
bodeplot(Psamples, w, legend=false, N=0, quantile=0)
bodeplot!(P*[delay(0.005) tf(0); tf(0) delay(0.005)], w) # Compare to the plant with a model of the delay
```

We see that the uncertain model set includes the model with the delay. Note how this approximation approach imparts some uncertainty also in the gain.

More details on this approach can be found in Skogestad sec. 7.4.

The other alternative is to use use sampled uncertain delays. The next example shows how we can create a system with an uncertain delay, where we know that the delay is an integer number of milliseconds between 1ms and 4ms.

```
using RobustAndOptimalControl, ControlSystemsBase, MonteCarloMeasurements, Plots, LinearAlgebra
unsafe_comparisons(true)
L = Particles(collect((1:4) ./ 1000)) # Uncertain time delay, an integer number of milliseconds between 1ms and 4ms
P = delay(L)*tf(1, [0.01, 1])
C = pid(2, 1)
w = exp10.(-1:0.01:4)
plot(
bodeplot(P, exp10.(-1:0.001:3), legend=false),
# plot(step(feedback(P, C), 0:0.0001:0.05), lab="L = " .* string.(P.Tau[].particles'), title="Disturbance response"), # This simulation requires using ControlSystems
nyquistplot(P*C, w[1:10:end], points=true, xlims=(-3.5, 2.5), ylims=(-5, 1.5), Ms_circles=[1.5, 2], alpha=1) # Note, the nyquistplot with uncertain coefficients requires manual selection of plot limits
)
```

Notice how the gain is completely certain, while the phase starts becoming very uncertain for high frequencies.

## Models of uncertain dynamics

This section goes through a number of uncertainty descriptions in block-diagram form and shows the equivalent transfer function appearing in feedback with the uncertain element. A common approach is to model an uncertain element as $W(s)\Delta$ where $||\Delta|| \leq 1$ and $W(s)$ is a frequency-dependent weighting function that is large for frequencies where the uncertainty is large.

#### Additive uncertainty

```
┌────┐ ┌────┐
┌►│ WΔ ├─┐ ┌─►│ WΔ ├──┐
│ └────┘ │ │ └────┘ │
│ │ │ │
┌───┐ │ ┌───┐ ▼ │ ┌──────┐ │
┌►│ C ├─┴─►│ P ├─+ │ │ C │ │
│ └───┘ └───┘ │ └─┤ ──── │◄┘
│ │ │ I+PC │
└────────────────┘ └──────┘
```

The system is made robust with respect to this uncertainty by making

\[W\dfrac{C}{I+PC} < 1\]

for all frequencies.

#### Multiplicative uncertainty

At the process output

```
┌────┐ ┌────┐
┌►│ WΔ ├┐ ┌─►│ WΔ ├──┐
│ └────┘│ │ └────┘ │
┌───┐ ┌───┐ │ ▼ │ │
┌►│ C ├──►│ P ├─┴───────+─► │ ┌──────┐ │
│ └───┘ └───┘ │ │ │ PC │ │
│ │ └─┤ ──── │◄┘
└───────────────────────┘ │ I+PC │
└──────┘
```

The system is made robust with respect to this uncertainty by making the complimentary sensitivity function $T$ satisfy

\[W\dfrac{PC}{I+PC} = WT < 1\]

for all frequencies.

This means that we must make the transfer function $T$ small for frequencies where the relative uncertainty is large. The relative uncertainty is always > 1 for sufficiently large frequencies, and this gives rise to the common adage of "apply lowpass filtering to avoid exciting higher-order dynamics at high frequencies".

This uncertainty representation was used in the examples above where we spoke about *multiplicative uncertainty*. For MIMO systems, uncertainty appearing on the plant input may behave different than if it appears on the plant output. In general, the loop-transfer function $L_o = PC$ denotes the *output loop-transfer function* (the loop is broken at the output of the plant) and $L_i = CP$ denotes the *input loop-transfer function* (the loop is broken at the input of the plant). For multiplicative uncertainty at the plant input, the corresponding transfer function to be constrained is

\[\dfrac{CP}{I+CP}W_i = TW_i < 1\]

with corresponding diagram

```
┌───┐ ┌────┐
┌►│ ΔW├─┐ ┌─►│ ΔW ├───┐
│ └───┘ │ │ └────┘ │
│ │ │ │
┌───┐ │ ▼ ┌───┐ │ ┌──────┐ │
┌►│ C ├─┴─────────►│ P ├─┐ │ │ CP │ │
│ └───┘ └───┘ │ └─┤ ──── │◄─┘
│ │ │ I+CP │
└────────────────────────┘ └──────┘
```

The input version represents uncertainty in the actuator, and is a particularly attractive object for analysis of SIMO systems, where this transfer function is SISO.

Skogestad and Postlethwaite Sec. 7.5.1

##### Example

See the example with Uncertain time delays above.

#### Additive feedback uncertainty

```
┌────┐ ┌────┐
┌─┤ WΔ │◄┐ ┌─►│ WΔ ├──┐
│ └────┘ │ │ └────┘ │
│ │ │ │
┌───┐ ▼ ┌───┐ │ │ ┌──────┐ │
┌►│ C ├─+─►│ P ├─┤ │ │ P │ │
│ └───┘ └───┘ │ └─┤ ──── │◄┘
│ │ │ I+PC │
└────────────────┘ └──────┘
```

The system is made robust with respect to this uncertainty by making

\[W\dfrac{P}{I+PC} < 1\]

for all frequencies.

This kind of uncertainty can represent uncertainty regarding presence of feedback loops, and uncertainty regarding the implementation of the controller (this uncertainty is equivalent to additive uncertainty in the controller).

#### Multiplicative feedback uncertainty

At the **process output**

```
┌────┐ ┌────┐
┌┤ WΔ │◄┐ ┌─►│ WΔ ├──┐
│└────┘ │ │ └────┘ │
┌───┐ ┌───┐ ▼ │ │ │
┌►│ C ├──►│ P ├─+───────┼─► │ ┌──────┐ │
│ └───┘ └───┘ │ │ │ I │ │
│ │ └─┤ ──── │◄┘
└───────────────────────┘ │ I+PC │
└──────┘
```

The system is made robust with respect to this uncertainty by making the (output) sensitivity function $S$ satisfy

\[W\dfrac{I}{I+PC} = WS < 1\]

for all frequencies.

This kind of uncertainty can represent uncertainty regarding which half plane poles are located. For frequencies where $W$ is larger than 1, poles can move from the left to the right half plane, and we thus need to make $S$ small (use lots of feedback) for those frequencies.

At the **process input**:

```
┌────┐ ┌───┐
┌─┤ WΔ │◄┐ ┌─►│ WΔ├───┐
│ └────┘ │ │ └───┘ │
│ │ │ │
┌───┐ ▼ │ ┌───┐ │ ┌──────┐ │
┌►│ C ├─+────────┴─►│ P ├┐ │ │ I │ │
│ └───┘ └───┘│ └─┤ ──── │◄┘
│ │ │ I+CP │
└────────────────────────┘ └──────┘
```

The system is made robust with respect to this uncertainty by making the (input) sensitivity function $S_i$ satisfy

\[\dfrac{I}{I+CP} W = S_i W < 1\]

for all frequencies.

Skogestad and Postlethwaite Sec. 7.5.3

#### Uncertainty through disturbances

Uncertainty can of course also be modeled as disturbances acting on the system. Similar to above, we may model disturbances as a signal that has $L_2$ norm less than 1, scaled by a weight $W(s)$. Additive, norm-bounded disturbances can never make a stable linear system unstable, the uncertainty does not appear *in the loop. The analysis of such disturbances can thus be focused on making the transfer function from the disturbance to the performance output small.

Some convenient facts when working with disturbances are that the $H_\infty$ norm of a transfer function $e = G_{ed}d$ is equal to the worst-case gain in $L_2$ norm of signals

\[\|G_{ed}\|_\infty = \sup_{d} \dfrac{||e||_2}{||d||_2}\]

and that the $H_2$ norm is equal to the gain in variance

\[\sigma_e^2 = \|G_{ed}\|_2 \sigma_d^2\]

#### Visualizing uncertainty

For any of the uncertainty descriptions above, we may plot the total loop gain excluding the uncertain element $\Delta$, that is, the weight $W(s)$ multiplied by the equivalent transfer function of the nominal part of the loop. For example, for multiplicative uncertainty at the plant output, we would plot $WT$ in a `sigmaplot`

and verify that all singular values are smaller than 1 for all frequencies. Alternatively, for SISO systems, we may plot $T$ and $W^{-1}$ in a Bode plot and verify that $T < W^{-1}$ for all frequencies. This latter visualization usually provides better intuition.

##### Example: Bode plot

Below, we perform this procedure for an multiplicative (relative) uncertainty model at the plant output. The uncertainty weight $W(s)$ is chosen to give 10% uncertainty at low frequencies and 10x uncertainty at high frequencies, indicating that we are absolutely oblivious to the behavior of the plant at high frequencies. This is often the case, either because identification experiments did not contain excitation for high frequencies, or because the plant had nonlinear behavior at higher frequencies.

```
using ControlSystemsBase, RobustAndOptimalControl, Plots
P = tf(1, [1, 2, 1]) # Plant model
C = pid(19.5, 0) # Controller
W = makeweight(0.1, 10, 10) # Low uncertainty (0.1) at low frequencies, large (10) at high frequencies.
bodeplot([comp_sensitivity(P, C), inv(W)], lab=["\$S\$" "\$W^{-1}\$"], linestyle=[:solid :dash], plotphase=false)
```

As long as the complimentary sensitivity function $T(s)$ stays below the inverse weight $W^{-1}(s)$, the closed-loop system is robust with respect to the modeled uncertainty.

##### Example: Nyquist plot

Continuing from the Bode-plot example above, we can translate the multiplicative weight $W(s)$ to a set of circles we could plot in the Nyquist diagram, one for each frequency, that covers the true open-loop system. For sampled representations of uncertainty, this is done using `fit_complex_perturbations`

, but here, we do it manually. For a given frequency $\omega$, the radius of the circle for an additive uncertainty in the loop gain is given by $|W(i\omega)|$, and for a multiplicative (relative) uncertainty, it is scaled by the loop gain $|W(i\omega) P(i\omega) C(i\omega)|$.^{[circ]} The center of the circle is simply given by the nominal value of the loop-gain $P(i\omega)C(i\omega)$.

```
w = exp10.(LinRange(-2, 2, 200))
centers = freqrespv(P*C, w)
radii = abs.(freqrespv(W*P*C, w))
nyquistplot(P*C, w, xlims=(-4,0.1), ylims=(-4,0.1))
nyquistcircles!(w, centers, radii)
```

If the plots above are created using the `plotly()`

backend, each circle is associated with hover information that is accessible by hovering the mouse over the plot. This indicates that the circle that touches the critical point is the one at $\omega \approx 4.5$, which coincides exactly with the point at thich the Bode plot above touches the inverse weight$W^{-1}$.

#### Converting between uncertainty descriptions

Any of the representations above, if modeled using uncertainty elements (`UncertainSS`

, `δc`

, `δr`

), may be converted to a sampled uncertainty representation using `rand(P_uncertain, 100)`

. The sampled representation can be further converted using `fit_complex_perturbations`

which results in a set of circles, additive or multiplicative (relative), one for each frequency considered, that covers the true system. These can be plotted in a Nyquist diagram using `nyquistcircles`

(see Example: Nyquist plot above).

A sampled representation can also be converted to a nominal value and a maximum $\nu$-gap, see Model-order reduction for uncertain models for an example of this

#### More references

Skogestad Sec. 8.5.3 contains result for moving uncertainty descriptions between input and output for MIMO systems as well as some additional forms of uncertainty descriptions, with robust stability conditions listed in Sec. 8.6.1.

- SkogestadSkogestad, "Multivariable Feedback Control: Analysis and Design"
- Doyle91Doyle, Packard, Zhou, "Review of LFTs, LMIs and μ".
`https://www.researchgate.net/publication/257200344_Review_of_LFTs_LMIs_and_mu`

- circA relative uncertainty does not apply to $P$ only, it appears after $P$ and thus models the relative uncertainty in the entire loop gain $PC$.