roban_uncertainty.qmd

---
title: "Uncertainty (in) modelling"
format:
    html:
        html-math-method: katex
        code-fold: true
        code-summary: "Show the code"
execute:
    enabled: true
jupyter: julia-1.10
---

Through this chapter we are stepping into the domain of robust control. We need to define a few keywords first.

::: {#def-uncertainty}
## Uncertainty
Deviation of the mathematical model of the system from the real (mechanical, electrial, chemical, biological, ...) system.
:::

::: {#def-robustness}
## Robustness
Insensitivity of specified properties of the system to the uncertainty.
:::

While these two terms are used in many other fields, here we are tailoring them to the discipline of control systems, in particular their model-based design. 

::: {#def-robust-control}
## Robust control
Not a single type of a control but rather a class of control design methods that aim to ensure robustness of the resulting control system. By convention, a robust controller is a fixed controller, typically designed for a nominal model. This is in contrast with an adaptive controller that adjusts itself in real time to the actual system.
:::

## Origins of uncertainty in models?

- Physical parameters are not known exactly (say, they are known to be within ±10% or ±3σ interval around the nominal value). 
- Even if the physical parameters are initially known with a high accuracy, they can evolve in time, unmeasured. 
- There may be variations among the individual units of the same product.
- If a nonlinear system is planned to be operated around a given operating point, it can be linearized aroud that operating point, which gives a nominal linear model. If the system is then operated in a significantly different operating point, the corresponding linear model is different from the nominal one.
- Our understanding of the underlying physics (or chemistry or biology or ...) is imperfect, hence our model is imperfect too. In fact, our understading can even be incorrect, in which case the model contains some discrepancies too. The imperfections of the model are typically observed at higher frequencies (referring to the frequency-domain modeling such as transfer functions).
- Even if we are able to eventually capture full dynamics of the system in a model, we may opt not to do so. We may want to keep the model simple, even if less accurate, because time invested into modelling is not for free. 
- Even if we can get a high-fidelity model with a reasonable effort, we may still prefer using a simpler (and less accurate) model for a controller design. The reason is that very often the complexity of the model used for model-based control design is reflected by the complexity of the controller – and high-complexity controllers are not particularly appreciated in industry.

## Models of uncertainty

There are several approaches to model the uncertainty (or, in other words, to characterize the uncertainty in the model). They all aim – in one way or another – to express that the controller has to deal no only with the single nominal system, for which it was designed, but a family of systems. Depending on the mathematical frameworks used for characterization of such a family, there are two major classes of approaches. 

- Worst-case models of uncertainty
- Probabilistic models of uncertainty

The former assumes sets of systems with no additional information about the structure of such sets. The latter imposes some probability structure on the set of systems – in other words, although in principle any member of the set possible, some may be more probable than the others. In this course we are focusing on the former, which is also the mainstream in the robust control literature, but note that the latter we already encountered while considering control for systems exposed to random disturbances, namely the LQG control. A possible viewpoint is that as a consequence of the random disturbance, the controller has to deal with a family of systems.

Another classification of models of uncertainty is according to the actual quantity that is uncertain. We distinguish these two

- Parametric uncertainty
- Frequency-dependent (aka dynamical) uncertainty
    - Unstructured uncertainty
    - Structured uncertainty

### Parametric uncertainty

This is obviously straightforward to state: some (real/physical) parameters are uncertain. The conceptually simplest way to characterize such uncertain parameters is by considering intervals instead of just single (nominal) values.

::: {#exm-pendulum-on-cart}
## A pendulum on a cart

$$
\begin{aligned}
{\color{red} m_\mathrm{l}} & \in [m_\mathrm{l}^{-},m_\mathrm{l}^{+}],\\
{\color{red} l} & \in [l^{-}, l^{+}],
\end{aligned}
$$

$$
\dot{\bm x}(t) =
\begin{bmatrix}
0 & 1 & 0 & 0\\
0 & 0 & \frac{\textcolor{red}{m_\mathrm{l}}}{m_\mathrm{c}} g & 0\\
0 & 0 & 0 & 1\\
0 & 0 & -\frac{(\textcolor{red}{m_\mathrm{l}}+m_\mathrm{c})g}{m_\mathrm{c}\textcolor{red}{l}} & 0
\end{bmatrix}
\bm x(t)
+
\begin{bmatrix}
0\\
\frac{1}{m_\mathrm{c}}\\
0\\
-\frac{1}{m_\mathrm{c}\textcolor{red}{l}}
\end{bmatrix}
u(t).
$$
:::

### Unstructured frequency-dependent uncertainty

Not only some enumerated physical parameters but even the order of the system can be uncertain. In other words, there may be some phenomena exhibitted by the system that is not captured by the model at all. Possibly some lightly damped modes, possibly some time delay here and there. The system contains uncertain dynamics. In the linear case, all this can be expressed by regarding the magnitude and phase responses uncertain without mapping these to actual physical parameters.  

![A whole subsystem is uncertain](figures/unstructured_uncertainty.png){#fig-unstructured-uncertainty width=40%}

A popular model for the uncertain subsystem is that of a transfer function $\Delta(s)$, about which we know only that it is stable and that its magnitude is bounded by 1
$$\boxed
{\sup_{\omega}|\Delta(j\omega)|\leq 1,\;\;\Delta \;\text{stable}. }
$$

But typically the uncertainty is higher at higher frequencies. This can be expressed by using some weighting function $w(\omega)$. 

For later theoretical and computational purposes we approximate the real weighting function using a low-order rational stable transfer function $W(s)$. That is, $W(j\omega)\approx w(\omega)$ for $\omega \in \mathbb R$, that is for $s=j\omega$ on the imaginary axis.

The ultimate transfer function model of the uncertainty is then  
$$\boxed{
W(s)\;\Delta(s),\quad \max_{\omega}|\Delta(j\omega)|\leq 1,\;\;\Delta\; \text{stable}. }
$$

#### $\mathcal H_\infty$ norm of an LTI system

##### H-infinity norm of an LTI system interpreted in frequency domain

::: {#def-hinf-norm-siso-system}
## $\mathcal H_\infty$ norm of a SISO LTI system
For a stable LTI system $G$ with a single input and single output, the $\mathcal H_\infty$ norm is defined as 
$$
\|G\|_{\infty} = \sup_{\omega\in\mathbb{R}}|G(j\omega)|.
$$
:::

::: {.callout-note}
## Why supremum and not maximum?
Supremum is uses in the definition because it is not guaranteed that the peak value of the magnitude frequency response is attained at a single frequency. Consider an example of a first-order system $G(s) = \frac{s}{Ts+1}$. The peak gain of $1/T$ is not attained at a single finite frequency.
:::

Having just defined the $\mathcal H_\infty$ norm, the uncertainty model can be expressed compactly as
$$\boxed{
W(s)\;\Delta(s),\quad \|\Delta(j\omega)\|\leq 1. }
$$

::: {.callout-note}
## $\mathcal H_\infty$ as a space of functions
$\mathcal H_\infty$ denotes a normed vector space of functions that are analytic in the closed extended right half plane (of the complex plane). In parlance of control systems, $\mathcal H_\infty$ is the space of proper and stable transfer functions. Poles on the imaginary axis are not allowed. The functions do not need to be rational, but very often we do restrict ourselves to rational functions, in which case we typically write such space as $\mathcal{RH}_\infty$.
:::

We now extend the concept of the $\mathcal H_\infty$ norm to MIMO systems. The extension is perhaps not quite intuitive – certainly it is not computed as the maximum of the norms of individual transfer functions, which may be the first guess.

::: {#def-hinf-norm-mimo-system}
## $\mathcal H_\infty$ norm of a MIMO LTI system
For a stable LTI system $\mathbf G$ with multiple inputs and/or multiple outputs, the $\mathcal H_\infty$ norm is defined as 
$$
\|\mathbf G\|_{\infty} = \sup_{\omega\in\mathbb{R}}\bar{\sigma}(\mathbf{G}(j\omega))
$$
where 
$\bar\sigma$ is the largest singular value.
:::

Here we include a short recap of singular values and singular value decomposition (SVD) of a matrix. Consider a matrix $\mathbf M$, possibly a rectangular one. It can be decomposed as a product of three matrices 
$$
\mathbf M = \mathbf U
\underbrace{
\begin{bmatrix}
\sigma_1 & & & &\\ 
  & \sigma_2 & & &\\
  & &\sigma_3 & &\\
 \\
  & & & & \sigma_n\\
\end{bmatrix}
}_{\boldsymbol\Sigma} 
\mathbf V^{*}.
$$

The two square matrices $\mathbf V$ and $\mathbf U$ are unitary, that is, 
$$
\mathbf V\mathbf V^*=\mathbf I=\mathbf V^*\mathbf V
$$
and
$$
\mathbf U\mathbf U^*=\mathbf I=\mathbf U^*\mathbf U.
$$

The nonnegative diagonal entries $\sigma_i \in \mathbb R_+, \forall i$ of the (possibly rectangular) matrix $\Sigma$ are called *singular values*. Commonly they are ordered in a nonincreasing order, that is
$$
\sigma_1\geq \sigma_2\geq \sigma_3\geq \ldots \geq \sigma_n.
$$

It is also a common notation to denote the largest singular value as $\bar \sigma$, that is, $\bar \sigma \coloneqq \sigma_1$.

##### $\mathcal{H}_{\infty}$ norm of an LTI system interpreted in time domain

We can also view the dynamical system $G$ with inputs and outputs as an operator mapping from some chosen space of functions to another space of functions. A popular model for these spaces are the spaces of square-integrable functions, denoted as $\mathcal{L}_2$, and sometimes interpreted as *bounded-energy* signals
$$
G:\;\mathcal{L}_2\rightarrow \mathcal{L}_2.
$$ 

It is a powerful fact that the $\mathcal{H}_{\infty}$ norm of the system is then defined as the induced norm of the corresponding operator
$$\boxed{
\|G(s)\|_{\infty} = \sup_{u(t)\in\mathcal{L}_2\setminus 0}\frac{\|y(t)\|_2}{\|u(t)\|_2}}.
$$

With the energy interpretation of the input and output variables, this system norm can also be interpreted as the worst-case *energy gain* of the system.

Scaling necessary to get any useful info from MIMO models! See Skogestad's book, section 1.4, pages 5--8. \url{https://folk.ntnu.no/skoge/book/ps/book1-3.pdf}


#### How does the uncertainty enter the model of the system?

##### Additive uncertainty

The transfer function of an uncertain system can be written as a sum of a nominal system and an uncertainty
$$
G(s) = \underbrace{G_0(s)}_{\text{nominal model}}+\underbrace{W(s)\Delta(s)}_{\text{additive uncertainty}}.
$$

The block diagram interpretation is in

![Additive uncertainty](figures/additive_uncertainty.png){#fig-additive-uncertainty width=40%}

The magnitude frequency response of the weighting filter $W(s)$ then serves as an upper bound on the absolute error in the magnitude frequency responses
$$
|G(j\omega)-G_0(j\omega)|<|W(j\omega)|\quad \forall \omega\in\mathbb R.
$$

##### Multiplicative uncertainty

$$
G(s) = (1+W(s)\Delta(s))\,G_0(s).
$$

The block diagram interpretation is in 

![Multiplicative uncertainty](figures/multiplicative_uncertainty.png){#fig-multiplicative-uncertainty width=55%}

::: {.callout-caution}
## For SISO transfer functions no need to bother about the order of terms in the products
Sice we are considering SISO transfer functions, the order of terms in the products is not important. We will have to be more alert to the order of terms when we move to MIMO systems.
:::

The magnitude frequency response of the weighting filter $W(s)$ then serves as an upper bound on the relative error in the magnitude frequency responses
$$\boxed
{\frac{|G(j\omega)-G_0(j\omega)|}{|G_0(j\omega)|}<|W(j\omega)|\quad \forall \omega\in\mathbb R.}
$$

::: {#exm-uncertain-first-order-system}
## Uncertain first-order delayed system
We consider a first-order system with a delay described by
$$
G(s) = \frac{k}{T s+1}e^{-\theta s}, \qquad 2\leq k,\tau,\theta,\leq 3.
$$

We now need to choose the nominal model $G_0(s)$ and then the uncertainty weighting filter $W(s)$. The nominal model corresponds to the nominal values of the parameters, therefore we must choose these. There is no single correct way to do this. Perhaps the most intuitive way is to choose the nominal values as the average of the bounds. But we can also choose the nominal values in a way that makes the nominal system simple. For example, for this system with a delay, we can even choose the nominal value of the delay as zero, which makes the nominal system a first-order system without delay, hence simple enough for application of some basic linear control system design methods. Of course, the price to pay is that the resulting model of an uncertain system, which is actually a set of systems, contains even models of a plant that were not prescribed. 

``` {julia}
using ControlSystems
using Plots

function uncertain_first_order_delayed()
    kmin = 2; kmax = 3; 
    θmin = 2; θmax = 3; 
    Tmin = 2; Tmax = 3;

    k₀ = (kmin+kmax)/2; 
    θ₀ = (θmin+θmax)/2;
    θ₀ = 0 
    T₀ = (Tmin+Tmax)/2;

    G₀ = tf(k₀,[T₀, 1])*delay(θ₀)  

    ω = exp10.(range(-2, 2, length=50))
    G₀ω = freqresp(G₀,ω)
    G₀ω = vec(G₀ω)

    EEω_db = [];
    for k in range(kmin,kmax,length=10)
        for θ in range(θmin,θmax,length=10)
            for T in range(Tmin,Tmax,length=10)
                G = tf(k,[T, 1])*delay(θ)  
                Gω = freqresp(G,ω)
                Gω = vec(Gω)            
                Eω = abs.(Gω-G₀ω)./abs.(G₀ω)
                Eω_db = 20 * log10.(Eω)
                push!(EEω_db,Eω_db)
            end
        end
    end

    plot(ω,EEω_db,xscale=:log10,label="",ylims=(-40,20))
    xlabel!("Frequency [rad/s]")
    ylabel!("Relative error [dB]")
end
uncertain_first_order_delayed()
```

Now we need to find some upper bound on the relative error. Simplicity is a virtue here too, hence we are looking for a rational filter of very low order, say 1 or 2. Speaking of the first-order filter, one useful way to format it is
$$
\boxed{
W(s) = \frac{\tau s+r_0}{(\tau/r_{\infty})s+1}}
$$
where $r_0$ is uncertainty at steady state, $1/\tau$ is the frequency, where the relative uncertainty reaches 100\%, $r_{\infty}$ is relative uncertainty at high frequencies, often $r_{\infty}\geq 2$.

For our example, the parameters of the filter are in the code below and the frequency response follows.
``` {julia}
using ControlSystems
using Plots

function uncertain_first_order_delayed_with_weights()
    kmin = 2; kmax = 3; 
    θmin = 2; θmax = 3; 
    Tmin = 2; Tmax = 3;

    k₀ = (kmin+kmax)/2; 
    θ₀ = (θmin+θmax)/2;
    θ₀ = 0 
    T₀ = (Tmin+Tmax)/2;

    G₀ = tf(k₀,[T₀, 1])*delay(θ₀)  

    ω = exp10.(range(-2, 2, length=50))
    G₀ω = freqresp(G₀,ω)
    G₀ω = vec(G₀ω)

    EEω_db = [];
    for k in range(kmin,kmax,length=10)
        for θ in range(θmin,θmax,length=10)
            for T in range(Tmin,Tmax,length=10)
                G = tf(k,[T, 1])*delay(θ)  
                Gω = freqresp(G,ω)
                Gω = vec(Gω)            
                Eω = abs.(Gω-G₀ω)./abs.(G₀ω)
                Eω_db = 20 * log10.(Eω)
                push!(EEω_db,Eω_db)
            end
        end
    end

    plot(ω,EEω_db,xscale=:log10,label="",ylims=(-40,20))
    xlabel!("Frequency [rad/s]")
    ylabel!("Relative error [dB]")

    τ = 1/0.25
    r₀ = 0.2
    r∞ = 10^(8/20)
    W = tf([τ, r₀],[τ/r∞, 1])
    magW, phaseW = bode(W,ω)
    plot!(ω,20*log10.(vec(magW)),xscale=:log10,lw=3,color=:red,label="W")

    # W2 = W*tf([1, 1.6, 1],[1, 1.4, 1]); 
    # magW2, phaseW2 = bode(W2,ω)
    # plot!(ω,20*log10.(vec(magW2)),xscale=:log10,lw=3,color=:blue,label="W₂")
end
uncertain_first_order_delayed_with_weights()
```

Obviously the filter does not capture the family of systems perfectly. It is now up to the control engineer to decide if this is a problem. If yes, if the control design should be really robust against all uncertainties in the considered set, some more complex (higher-order) filter is needed to described the uncertainty more accurately. The source code shows (in commented lines) one particular candidate, but in general the whole problem boils down to designing a stable filter with a prescribed magnitude frequency response. 
:::




##### Inverse additive uncertainty

...

##### Inverse multiplicative uncertainty

...

##### Linear fractional transformation (LFT)

For a matrix $\mathbf P$ sized $(n_1+n_2)\times(m_1+m_2)$ and divided into blocks like
$$
\mathbf P=
\begin{bmatrix}
\mathbf P_{11} & \mathbf P_{12}\\
\mathbf P_{21} & \mathbf P_{22}
\end{bmatrix},
$$
and a matrix $\mathbf K$ sized $m_2\times n_2$, the *lower LFT* of $\mathbf P$ with respect to $\mathbf K$ is
$$
\boxed{
\mathcal{F}_\mathbf{l}(\mathbf P,\mathbf K) = \mathbf P_{11}+\mathbf P_{12}\mathbf K(\mathbf I-\mathbf P_{22}\mathbf K)^{-1}\mathbf P_{21}}.
$$

It can be viewed as a feedback interconnection of the plant $\mathbf P$ and the controller $\mathbf K$, in which not all plant inputs are used as control inputs and not all plant outputs are measured, as depicted in @fig-lower-LFT

![Lower LFT of $\mathbf P$ with respect to $\mathbf K$](figures/lower_LFT.png){#fig-lower-LFT width=30%}

Similarly, for a matrix $\mathbf N$ sized $(n_1+n_2)\times(m_1+m_2)$ and a matrix $\boldsymbol\Delta$ sized $m_1\times n_1$, the *upper LFT* of $\mathbf N$ with respect to $\mathbf K$ is
$$
\boxed{
\mathcal{F}_\mathbf{u}(\mathbf N,\boldsymbol\Delta) = \mathbf N_{22}+\mathbf N_{21}\boldsymbol\Delta(\mathbf I-\mathbf N_{11}\boldsymbol\Delta)^{-1}\mathbf N_{12}}.
$$

It can be viewed as a feedback interconnection of the nominal plant $\mathbf N$ and the uncertainty block $\boldsymbol\Delta$, as depicted in @fig-upper-LFT

![Upper LFT of $\mathbf N$ with respect to $\boldsymbol \Delta$](figures/upper_LFT.png){#fig-upper-LFT width=30%}

Here we already anticipated MIMO uncertainty blocks. One motivation for them is explained in the very next section on structured uncertainties, another one is appearing once we start formulating robust performance within the same analytical framework as robust stability.


::: {.callout-note}
## Which is lower and which is upper is a matter of convention, but a useful one
Our usage of the lower LFT for a feedback interconnection of a (generalized) plant and a controller and the upper LFT for a feedback interconnection of a nominal system and and uncertainty is completely arbitrary. We could easily use the lower LFT for the uncertainty. But it is a convenient convention to adhere to. The more so that it allows for the combination of both as in the diagram @fig-lower-upper-LFT below, which corresponds to composition of the two LFTs.

![Combination of the lower and upper LFT](figures/lower_upper_LFT.png){#fig-lower-upper-LFT width=30%}
:::

### Structured frequency-domain uncertainty

Not just a single $\Delta(s)$ but several $\Delta_i(s)$, $i=1,\ldots,n$ are considered. Some of them scalar-valued, some of them matrix-valued. 

In the upper LFT, all the individual $\Delta_i$s are collected into a single overall $\boldsymbol \Delta$, which then exhibits some structure. Typically it is block-diagonal as in 
$$
\boldsymbol\Delta = 
\begin{bmatrix}
\Delta_1& 0 & \ldots & 0\\
0 & \Delta_2 & \ldots & 0\\
\vdots\\
0 & 0 & \ldots & \boldsymbol\Delta_n
\end{bmatrix},
$$
with each block (including the MIMO blocks) satisfying the usual condition
$$
\|\Delta_i\|_{\infty}\leq 1, \; i=1,\ldots, n.
$$

#### Structured singular value (SSV, $\mu$, mu)

With this structured uncertainty, how does the small gain theorem look like?