▲
│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 frequencyWd: 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:
external_inputs = [
:do, :u
]
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; external_outputs, external_inputs)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 zExtendedStateSpace{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 modelLet's perform some checks to see that we have a sound problem specification
hinfassumptions(Gsyn)truewe also check the equivalences we know should hold
system_mapping(Gsyn) == P.systruesminreal(G[:uw, :u].sys) == Wu.systrueThese 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-8truehinfnorm2(G[:e, :u].sys - (We * P).sys)[1] < 1e-8trueNow, let's synthesize!
K, γ = hinfsynthesize(Gsyn, γrel=1.1, transform=false)
γ, K(0.19725514012208165, 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.19725514012208165, 2.557529082231149)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))