▲
│uw      │dw              │do
┌─┴─┐    ┌─▼─┐            ┌─▼─┐
│   │    │   │            │   │
│ Wu│    │ Wi│            │ Wd│
│   │    │   │            │   │
┌────┐       ┌─────┐ └─▲─┘    └───┘   ┌─────┐  └─┬─┘    ┌───┐
r  │    │       │     │   │      di▼     │     │   d│  y   │   │ e
──►│ Wr ├──►+──►│  K  ├───┴────────+────►│  P  ├────▼───┬─►│ We├──►
│    │   ▲   │     │   u              │     │ Py     │  │   │
└────┘   │-  └─────┘                  └─────┘        │  └───┘
│                                           │
│                                           │  ┌───┐
│                                           ▼  │   │
└───────────────────────────────────────────+◄─┤ Wn│◄─
n │   │ ni
└───┘

In this example we will show how to manually construct the following transfer function, the $H_\infty$ norm of which is commonly a performance metric for $H_\infty$ synthesis.

$$$\begin{bmatrix} W_e S_o W_d \\ W_u K S_o W_d \end{bmatrix}$$$

(Note, the function hinfpartition exists to make special instances of this kind of partitioning smoother, this example shows how general transfer matrices can be constructed for more complicated $H_\infty$ designs)

We note that the first transfer function is the transfer function from $d_o \rightarrow e$ and the second transfer function is from $d_o \rightarrow u_w$. It's thus enough for this example to define the blocks $P,W_e,W_d,W_u$. To form the feedback interconnection, we also need to define a summation node after $P$. Note, $K$ is the output of the synthesis procedure, we're thus going to break the loop at $K$ and place the inputs and outputs to $K$ ($u, y$) in the system_mapping, while the external inputs and outputs ($w, z$) go in the performance_mapping.

Using the connect function, we can then form the closed-loop transfer function automatically, without forming the intermediate output sensitivity function $S_o$.

The code below constructs all transfer functions using the same names as in the block diagram above. Before we start, it really helps to have a block diagram in front of you to not get lost. We redraw the diagram focusing on the relevant signals and blocks only

               ▲
│uw              │do
┌─┴─┐            ┌─▼─┐
│   │            │   │
│ Wu│            │ Wd│
│   │            │   │
┌─────┐ └─▲─┘   ┌─────┐  └─┬─┘    ┌───┐
│     │   │     │     │   d│  y   │   │ e
┌───►│  K  ├───┴─────┤  P  ├────▼───┬─►│ We├──►
│    │     │   u     │     │        │  │   │
│    └─────┘         └─────┘        │  └───┘
│                                   │
│                                   │
│                                   │
└───────────────────────────────────┘
using RobustAndOptimalControl, ControlSystemsBase, Plots
P = named_ss(ss(tf(1, [1, 0.02, 1])), :P)
We = named_ss(10P.sys,  :We, u=:y, y=:e)  # Weigh sensitivity by magnitude of plant
Wu = named_ss(makeweight(1e-3, 10, 10), :Wu, u=:Wu, y=:uw) # Above ω=100, we want to limit the control effort
Wd = named_ss(makeweight(10, 1, 1e-3),  :Wd, u=:do, y=:d) # d is low frequency
Wd: NamedStateSpace{Continuous, Float64}
A =
-0.01
B =
0.25
C =
0.39996000000000004
D =
0.001

Continuous-time state-space model
With state  names: Wdx
input  names: do
output names: d


We also create the summation node after $P$

sumP = sumblock("y = Py + d")
sumblock: NamedStateSpace{Continuous, Float64}
D =
1.0  1.0

Continuous-time state-space model
With state  names:
input  names: Py d
output names: y


We also need a splitter to split the external input :u (remember, the loop was broken at $K$) into two signals with unique names.

split_u = splitter(:u, 2)
splitter: NamedStateSpace{Continuous, Float64}
D =
1.0
1.0

Continuous-time state-space model
With state  names:
input  names: u
output names: u1 u2


The next step is to list all the connections

connections = [
:u1 => :Wu # splitter to input of Wu
:u2 => :Pu # splitter to input of P
:Py => :Py # P output to first input of sumblock
:d => :d   # output of Wd to second input of sumblock
:y => :y   # output of sumblock to input of We
]

When we specify the external inputs and outputs to connect, we include $y$ and $u$ since they are external from the view of connect:

w1 = [ # External inputs
:do, :u
]
z1 = [ # External outputs
:e, :uw, :y
]

We are now ready to form the system we want to minimize the norm of

G = connect([P,We,Wu,Wd,sumP,split_u], connections; z1, w1)
NamedStateSpace{Continuous, Float64}
A =
0.0   1.0    0.0   0.0     0.0                0.0
-1.0  -0.02   0.0   0.0     0.0                0.0
0.0   0.0    0.0   1.0     0.0                0.0
1.0   0.0   -1.0  -0.02    0.0                0.39996000000000004
0.0   0.0    0.0   0.0   -99.49879346007116   0.0
0.0   0.0    0.0   0.0     0.0               -0.01
B =
0.0     0.0
0.0     1.0
0.0     0.0
0.001   0.0
0.0    32.0
0.25    0.0
C =
0.0  0.0  10.0  0.0    0.0               0.0
0.0  0.0   0.0  0.0  -31.09026361897661  0.0
1.0  0.0   0.0  0.0    0.0               0.39996000000000004
D =
0.0     0.0
0.0    10.0
0.001   0.0

Continuous-time state-space model
With state  names: Px1##feedback#477 Px2##feedback#477 Wex1##feedback#477 Wex2##feedback#477 Wux##feedback#477 Wdx##feedback#477
input  names: do u
output names: e uw y


At this stage, it's good practice to check that the poles, inputs and outputs of $G$ looks correct, it's easy to forget some signal..

Before we perform the $H_\infty$ synthesis, we need to partition this system to indicate which outputs go to and from the controller and which are part of the performance specification

Gsyn = partition(G, u = [:u], y = [:y]) # You can provide either u or w, and either y or z
ExtendedStateSpace{Continuous, Float64}
A =
0.0   1.0    0.0   0.0     0.0                0.0
-1.0  -0.02   0.0   0.0     0.0                0.0
0.0   0.0    0.0   1.0     0.0                0.0
1.0   0.0   -1.0  -0.02    0.0                0.39996000000000004
0.0   0.0    0.0   0.0   -99.49879346007116   0.0
0.0   0.0    0.0   0.0     0.0               -0.01
B1 =
0.0
0.0
0.0
0.001
0.0
0.25
B2 =
0.0
1.0
0.0
0.0
32.0
0.0
C1 =
0.0  0.0  10.0  0.0    0.0               0.0
0.0  0.0   0.0  0.0  -31.09026361897661  0.0
C2 =
1.0  0.0  0.0  0.0  0.0  0.39996000000000004
D11 =
0.0
0.0
D12 =
0.0
10.0
D21 =
0.001
D22 =
0.0
Continuous-time extended state-space model

Let's perform some checks to see that we have a sound problem specification

hinfassumptions(Gsyn)
true

we also check the equivalences we know should hold

system_mapping(Gsyn) == P.sys
true
G[:uw, :u].sys == Wu.sys
false

These will not be identical, the realization might differ, but they should represent the same system

hinfnorm2(G[:e, :do].sys - (We * Wd).sys)[1] < 1e-8
true
hinfnorm2(G[:e, :u].sys - (We * P).sys)[1] < 1e-8
false

Now, let's synthesize!

K, γ = hinfsynthesize(Gsyn, γrel=1.1, transform=false)
γ, K
(0.1972551401220788, StateSpace{Continuous, Float64}
A =
0.0                  1.0                    0.0                   0.0                  0.0                   0.0
-5.029564643626379   -1.4683706824648812    -3.344193594283276    -6.880021065211896    2.283886347486903    -1.7120108840893369
0.0                  0.0                    0.0                   1.0                  0.0                   0.0
0.0                  0.0                   -1.0                  -0.02                 0.0                   0.0
-128.9460685960441    -46.3478618388762     -107.01419501706484   -220.16067408678066   -26.414430340490256   -54.78434829085878
-250.00000000000006     0.0                    0.0                   0.0                  0.0                -100.00000000000004
B =
0.0
0.0
0.0
1.0
0.0
250.00000000000006
C =
-4.029564643626379  -1.4483706824648812  -3.344193594283276  -6.880021065211896  2.283886347486903  -1.7120108840893369
D =
0.0

Continuous-time state-space model)

We may now form the closed-loop system and check the calculated $H_\infty$ norm

Gcl = lft(Gsyn, K)
hinfnorm2(Gcl) # Should be equal to γ
(0.1972551401220788, 2.557529082231186)

The calculated $\gamma$ is quite small, indicating that our performance specifications were easy to satisfy. Let's plot the resulting performance mapping with the lower loop closed around $K$ as well as some interesting sensitivity functions

w = exp10.(LinRange(-3, 3, 300))
bodeplot(Gcl, w, plotphase=false, title=["do -> e" "do -> uw"])
S, PS, CS, T = gangoffour(P, K)
specificationplot([S, CS, T], ylims=(1e-3, 1e2))