HH_\infty control design

In this tutorial, we will design a controller for a DC-servo (electrical motor).

The servo takes a torque command as input and we are interested in controlling the angle of the output shaft. An approximate model for a DC servo, valid for low frequencies, is

τ=k(Jϕ¨+cϕ˙)\tau = k ( J \ddot{\phi} + c\dot{\phi} )

where τ\tau is the torque and ϕ\phi is the angle. The parameter JJ denotes the inertia of the motor and load, cc is a viscous friction parameter and kk is a torque constant (gain).

This is a simple SISO example with a pole in the origin. Poles on the stability boundary are problematic for several numerical routines, and this tutorial demonstrates a common trick applied in these situations.

We start by defining the process model. We will use the parameter values J=1,c=0.12,k=11.2J=1, c=0.12, k=11.2 which yields the transfer function

G(s)=ks(Js+c)=11.2s2+0.12s+0G(s) = \dfrac{k}{s(Js + c)} = \dfrac{11.2}{s^2 + 0.12s + 0}

using DyadControlSystems, Plots
Gtrue = tf([11.2], [1, 0.12, 0])

When designing a controller using HH_\infty synthesis, we formally specify an optimization problem where the cost function is the HH_\infty norm of a suitably chosen transfer function, and the optimization variable is the controller. In this example we will optimize

minimizeKWSSWUKS\operatorname{minimize}_K \begin{Vmatrix} W_S S \\ W_U KS \end{Vmatrix}_\infty

where KK is the controller and SS is the sensitivity function (1+GK)1(1+GK)^{-1}. The transfer functions WSW_S and WUW_U are weight functions that emphasize different frequency ranges, for example, if WS(iω)W_S(i\omega) is large for a particular frequency ω\omega, then SS is forced to be small at ω\omega in order for the HH_\infty norm to be small.

# Sensitivity weight function
WS = makeweight(1e5, 0.05, 0.5) |> tf

# Output sensitivity weight function. Increase this value to penalize controller effort more
WU = ss(1.0)

# Complementary sensitivity weight function
WT = [] # We do not put any weight on T in this example

To solve this problem, we partition the system PP in a two-input, two output configuration such that lftl(P,K)\operatorname{lft}_{l}(P,K) forms the system we want to minimize the HH_\infty of. To help with this partitioning, we have the function hinfpartition:

P = hinfpartition(Gtrue, WS, WU, WT)

Before solving, we may check if the synthesis problem is feasible

hinfassumptions(P, verbose=true)
false

The problem is not feasible, in this case due to the integrator pole on the stability margin. We thus modify the plant description by moving the integrator pole in the origin slightly ($\vareplsilon$) into the stable region

ε = 1e-4
G = tf([11.2], [1, 0.12]) * tf([1], [1, ε])

If we now perform the assumption check again, it passes

P = hinfpartition(G, WS, WU, WT)
hinfassumptions(P)
true

Synthesize the H-infinity optimal controller

With the problem properly defined, we may call hinfsyn_lmi to solve it. The result is a controller K and a performance index γ\gamma:

K, γ = hinfsyn_lmi(P, γrel=1.05, ftype=BigFloat)
γ
3.146036497316298

The achieved performance level is indicated by γ\gamma, this number is the norm we are optimizing and it should be as low as possible.

For verification purposes, we may extract some transfer functions defining common sensitivity functions:

Pcl, S, KS, T = hinfsignals(P, G, K)

We may also verify that the closed-loop system has HH_\infty norm γ\gamma

Pcl == lft(P, K)
isapprox(hinfnorm2(Pcl)[1], γ, rtol=1e-5)
true

Plot the specifications

The resulting sensitivity functions can be plotted together with the inverse weighting functions multiplied by γ\gamma to get a feeling for where the constraints are active.

specificationplot([S, KS, T], [WS, WU, WT], γ)
Example block output

In this case, the noise amplification constraint is active between about 100101.510^0 - 10^{1.5} rad/s and the sensitivity function is pushed down for lower frequencies. In this case, the complimentary sensitivity function T=(1+GK)1GKT = (1+GK)^{-1}GK has desireable properties without us specifying and penalizing this function in the optimization problem. If we would like to push this function down at certain frequencies, we may specify a non-empty weight WTW_T as well.