A brief and practical introduction to \(\mathcal{H}_\infty\) Control

The typical methodology for Model Based Control techniques is schematically shown in Figure 1.

It consists of three steps:

1. Identification or modeling: a mathematical model \(G(s)\) representing the plant dynamics is obtained
2. Translate the specifications into mathematical criteria:
   • Specifications: Response Time, Noise Rejection, Maximum input amplitude, Robustness, …
   • Mathematical Criteria: Cost Function, Shape of transfer function, Phase/Gain margin, Roll-Off, …
3. Synthesis: research of a controller \(K(s)\) that satisfies the specifications for the model of the system

In this document, we will suppose a model of the plant is available (step 1 already performed), and we will focus on steps 2 and 3.

In Section 2, steps 2 and 3 will be described for a control techniques called classical (open-)loop shaping.

Then, steps 2 and 3 for the \(\mathcal{H}_\infty\) Loop Shaping of closed-loop transfer functions will be discussed in Sections 4, 5 and 6.

Many different model based control techniques have been developed since the birth of classical control theory in the ’30s.

Classical control methods were developed starting from 1930 based on tools such as the Laplace and Fourier transforms. It was then natural to study systems in the frequency domain using tools such as the Bode and Nyquist plots. Controllers were manually tuned to optimize criteria such as control bandwidth, gain and phase margins.

The ’60s saw the development of control techniques based on a state-space. Linear algebra and matrices were used instead of the frequency domain tool of the class control theory. This allows multi-inputs multi-outputs systems to be easily treated. Kalman introduced the well known Kalman estimator as well the notion of optimality by minimizing quadratic cost functions. This set of developments is loosely termed Modern Control theory.

By the 1980’s, modern control theory was shown to have some robustness issues and to lack the intuitive tools that the classical control methods were offering. This lead to a new control theory called Robust control that blends the best features of classical and modern techniques. This robust control theory is the subject of this document.

The three presented control methods are compared in Table 1.

Note that in parallel, there have been numerous other developments, including non-linear control, adaptive control, machine-learning control just to name a few.

Let’s consider the model shown in Figure 3. It could represent a suspension system with a payload to position or isolate using an force actuator and an inertial sensor. The notations used are listed in Table 2.

Having obtained \(G(s)\) and \(G_d(s)\), we can transform the system shown in Figure 3 into a classical feedback form as shown in Figure 7.

Let’s define the system parameters on Matlab.

1: k = 1e6; % Stiffness [N/m]
2: c = 4e2; % Damping [N/(m/s)]
3: m = 10; % Mass [kg]

And now the system dynamics \(G(s)\) and \(G_d(s)\) (their bode plots are shown in Figures 5 and 6).

4: G = 1/(m*s^2 + c*s + k); % Plant
5: Gd = (c*s + k)/(m*s^2 + c*s + k); % Disturbance

This synthesis method is widely used as many characteristics of the closed-loop system depend on the shape of the open loop gain \(L(s)\) such as:

• Performance: \(L\) large
• Good disturbance rejection: \(L\) large
• Limitation of measurement noise on plant output: \(L\) small
• Small magnitude of input signal: \(K\) and \(L\) small
• Nominal stability: \(L\) small (RHP zeros and time delays)
• Robust stability: \(L\) small (neglected dynamics)

The Open Loop shape is usually done manually has the loop gain \(L(s)\) depends linearly on \(K(s)\) \eqref{eq:loop_gain}.

\(K(s)\) then consists of a combination of leads, lags, notches, etc. such that \(L(s)\) has the wanted shape (an example is shown in Figure 8).

sisotool(G)

In order to have the wanted Roll-off, two integrators are used, a lead is also added to have sufficient phase margin.

The obtained controller is shown below, and the bode plot of the Loop Gain is shown in Figure 9.

K = 14e8 * ... % Gain
    1/(s^2) * ... % Double Integrator
    1/(1 + s/2/pi/40) * ... % Low Pass Filter
    (1 + s/(2*pi*10/sqrt(8)))/(1 + s/(2*pi*10*sqrt(8))); % Lead

And we can verify that we have the wanted stability margins:

[Gm, Pm, ~, Wc] = margin(G*K)

The Open Loop Shaping synthesis can be performed using the \(\mathcal{H}_\infty\) Synthesis.

Even though we will not go into details, we will provide one example.

Using Matlab, the \(\mathcal{H}_\infty\) Loop Shaping Synthesis can be performed using the loopsyn command:

K = loopsyn(G, Gd);

where:

• G is the (LTI) plant
• Gd is the wanted loop shape
• K is the synthesize controller

Let’s reuse the previous plant.

Translate the specification into the wanted shape of the open loop gain.

• Performance: Bandwidth of approximately 10Hz: \(|L_w(j2 \pi 10)| = 1\)
• Noise Attenuation: Roll-off of -40dB/decade past 30Hz
• Robustness: Gain margin > 3dB and Phase margin > 30 deg

Lw = 2.3e3 * ...
     1/(s^2) * ... % Double Integrator
     (1 + s/(2*pi*10/sqrt(3)))/(1 + s/(2*pi*10*sqrt(3))); % Lead

The \(\mathcal{H}_\infty\) optimal open loop shaping synthesis is performed using the loopsyn command:

[K, ~, GAM] = loopsyn(G, Lw);

Let’s briefly analyze the obtained controller which bode plot is shown in Figure 10:

• two integrators are used at low frequency to have the wanted low frequency high gain
• a lead is added centered with the crossover frequency to increase the phase margin
• a notch is added at the resonance of the plant to increase the gain margin (this is very typical of \(\mathcal{H}_\infty\) controllers, and can be an issue, more info on that latter)

The obtained Loop Gain is shown in Figure 11 and matches the specified one by a factor \(\gamma\).

Let’s now compare the obtained stability margins of the \(\mathcal{H}_\infty\) controller and of the manually developed controller in Table 3.

Why optimizing the \(\mathcal{H}_\infty\) norm of transfer functions is a pertinent choice will become clear when we will translate the typical control specifications into the \(\mathcal{H}_\infty\) norm of transfer functions.

Then applying the \(\mathcal{H}_\infty\) synthesis to a plant, the engineer work usually consists of the following steps

1. Write the problem as standard \(\mathcal{H}_\infty\) problem
2. Translate the specifications as \(\mathcal{H}_\infty\) norms of transfer functions
3. Make the synthesis and analyze the obtain controller
4. Reduce the order of the controller for implementation

Note that there are many ways to use the \(\mathcal{H}_\infty\) Synthesis:

• Traditional \(\mathcal{H}_\infty\) Synthesis (hinfsyn doc)
• Open Loop Shaping \(\mathcal{H}_\infty\) Synthesis (loopsyn doc)
• Mixed Sensitivity Loop Shaping (mixsyn doc)
• Fixed-Structure \(\mathcal{H}_\infty\) Synthesis (hinfstruct doc)
• Signal Based \(\mathcal{H}_\infty\) Synthesis, and many more…

The first step when applying the \(\mathcal{H}_\infty\) synthesis is usually to write the problem as a standard \(\mathcal{H}_\infty\) problem. This consist of deriving the Generalized Plant for the current problem. It makes things much easier for the following steps.

The generalized plant, usually noted \(P(s)\), is shown in Figure 13. It has two inputs and two outputs (both could contains many signals). The meaning of the inputs and outputs are summarized in Table 4.

Note that this generalized plant is as its name implies, quite general. It can indeed represent feedback as well as feedforward control architectures.

Once the generalized plant is obtained, the \(\mathcal{H}_\infty\) synthesis problem can be stated as follows:

Note that the closed-loop transfer function from \(w\) to \(z\) is:

Using Matlab, the \(\mathcal{H}_\infty\) Synthesis applied on a Generalized plant can be applied using the hinfsyn command (documentation):

K = hinfsyn(P, nmeas, ncont);

where:

• P is the generalized plant transfer function matrix
• nmeas is the number of sensed output (size of \(v\))
• ncont is the number of control signals (size of \(u\))
• K obtained controller that minimized the \(\mathcal{H}_\infty\) norm from \(w\) to \(z\)

The procedure to convert a typical control architecture as the one shown in Figure 15 to a generalized Plant is as follows:

1. Define signals (\(w\), \(z\), \(u\) and \(v\)) of the generalized plant
2. Remove \(K\) and rearrange the inputs and outputs to match the generalized configuration shown in Figure 13

P = [1 -G;
     0  1;
     1 -G]
P.InputName = {'w', 'u'};
P.OutputName = {'e', 'u', 'v'};

As the performances of a controlled system depend on the closed loop transfer functions, it is very important to derive these closed-loop transfer functions as a function of the plant \(G(s)\) and controller \(K(s)\).

And we have:

Thus, for reference tracking, we have to shape the closed-loop transfer function from \(r\) to \(\epsilon\), that is the sensitivity function \(S(s)\). Similarly, to reduce the effect of measurement noise \(n\) on the output \(y\), we have to act on the complementary sensitivity function \(T(s)\).

The sensitivity function is indisputably the most important closed-loop transfer function of a feedback system. In this section, we will see how the shape of the sensitivity function will impact the performances of the closed-loop system.

Suppose we have developed a “reference” controller \(K_r(s)\) and made three small changes to obtained three controllers \(K_1(s)\), \(K_2(s)\) and \(K_3(s)\). The obtained sensitivity functions are shown in Figure 18 and the corresponding step responses are shown in Figure 19.

The comparison of the sensitivity functions shapes and their effect on the step response is summarized in Table 6.

Let’s start by an example demonstrating why the phase and gain margins might not be good indicators of robustness.

w0 = 2*pi*100;
xi = 0.1;
k = 1e7;

Gt = 1/k*(s/w0/4 + 1)/(s^2/w0^2 + 2*xi*s/w0 + 1);

Kt = 1.2e6*(s + w0)/s;

xi = 0.03;

isstable(feedback(Gr,K))

0

Let’s now determine a new robustness indicator based on the Nyquist Stability Criteria.

From the Nyquist stability criteria, it is clear that we want \(L(j\omega)\) to be as far as possible from the \(-1\) point (called the unstable point) in the complex plane. This minimum distance is called the module margin.

As expected from Figure 23, there is a close relationship between the module margin and the gain and phase margins. We can indeed show that for a given value of the module margin \(\Delta M\), we have:

Let’s now try to express the Module margin \(\Delta M\) as an \(\mathcal{H}_\infty\) norm of a closed-loop transfer function:

Note that this is why large peak value of \(|S(j\omega)|\) usually indicate robustness problems. And we know understand why setting an upper bound on the magnitude of \(S\) is generally a good idea.

Interpretation of the \(\mathcal{H}_\infty\) norm of systems:

• frequency by frequency attenuation / amplification

Let’s note \(G_t(s)\) the closed-loop transfer function from \(w\) to \(z\).

Consider an input sinusoidal signal \(w(t) = \sin\left( \omega_0 t \right)\), then the output signal \(z(t)\) will be equal to: \[ z(t) = A \sin\left( \omega_0 t + \phi \right) \] with:

• \(A = |G_t(j\omega_0)|\) is the magnitude of \(G_t(s)\) at \(\omega_0\)
• \(\phi = \angle G_t(j\omega_0)\) is the phase of \(G_t(s)\) at \(\omega_0\)

Noise Attenuation: typical wanted shape for \(T\)

If the \(\mathcal{H}_\infty\) synthesis is applied on the generalized plant \(P(s)\) shown in Figure 24, it will generate a controller \(K(s)\) such that the \(\mathcal{H}_\infty\) norm of closed-loop transfer function from \(r\) to \(\epsilon\) is minimized. This closed-loop transfer function actually correspond to the sensitivity function. Therefore, it will minimize the the \(\mathcal{H}_\infty\) norm of the sensitivity function: \(\|S\|_\infty\).

However, as the \(\mathcal{H}_\infty\) norm is the maximum peak value of the transfer function’s magnitude, this synthesis is quite useless and clearly does not allow to shape the norm of \(S(j\omega)\) over all frequencies.

Let’s now show how this is equivalent as shaping the sensitivity function:

Pw = [Ws -Ws*G;
      1  -G];

Pw = blkdiag(Ws, 1)*P;

Weighting function included in the generalized plant must be proper, stable and minimum phase transfer functions.

Matlab is providing the makeweight function that allows to design first-order weights by specifying the low frequency gain, high frequency gain, and the gain at a specific frequency:

W = makeweight(dcgain,[freq,mag],hfgain)

with:

• dcgain: low frequency gain
• [freq,mag]: frequency freq at which the gain is mag
• hfgain: high frequency gain

Ws = makeweight(1e2, [2*pi*10, 1], 1/2);

function [W] = generateWeight(args)
    arguments
        args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1
        args.G1 (1,1) double {mustBeNumeric, mustBePositive} = 10
        args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1
        args.wc (1,1) double {mustBeNumeric, mustBePositive} = 2*pi
        args.n  (1,1) double {mustBeInteger, mustBePositive} = 1
    end

    if (args.Gc <= args.G0 && args.Gc <= args.G1) || (args.Gc >= args.G0 && args.Gc >= args.G1)
        eid = 'value:range';
        msg = 'Gc must be between G0 and G1';
        throwAsCaller(MException(eid,msg))
    end

    s = zpk('s');

    W = (((1/args.wc)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.G1)^(2/args.n)))*s + (args.G0/args.Gc)^(1/args.n))/((1/args.G1)^(1/args.n)*(1/args.wc)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.G1)^(2/args.n)))*s + (1/args.Gc)^(1/args.n)))^args.n;

end

W1 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 1);
W2 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 2);
W3 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 3);

Let’s design a controller using the \(\mathcal{H}_\infty\) synthesis that fulfils the following requirements:

1. Bandwidth of at least 10Hz
2. Small static errors for step responses
3. Robustness: Large module margin \(\Delta M > 0.5\) (\(\Rightarrow \Delta G > 2\) and \(\Delta \phi > 29^o\))

As usual, the plant used is the one presented in Section 1.3.

Ws = makeweight(1e3, [2*pi*10, sqrt(2)], 1/2);

Ws = generateWeight('G0', 1e3, ...
                    'G1', 1/2, ...
                    'Gc', sqrt(2), 'wc', 2*pi*10, ...
                    'n', 2);

Let’s say we came up with the following weighting function:

Ws = generateWeight('G0', 1e3, ...
                    'G1', 1/2, ...
                    'Gc', sqrt(2), 'wc', 2*pi*10, ...
                    'n', 2);

The weighting function is then added to the generalized plant.

P = [1 -G;
     1 -G];
Pw = blkdiag(Ws, 1)*P;

And the \(\mathcal{H}_\infty\) synthesis is performed on the weighted generalized plant.

K = hinfsyn(Pw, 1, 1, 'Display', 'on');

Test bounds:  0.5 <=  gamma  <=  0.51

 gamma        X>=0        Y>=0       rho(XY)<1    p/f
5.05e-01     0.0e+00     0.0e+00     3.000e-16     p
Limiting gains...
5.05e-01     0.0e+00     0.0e+00     3.461e-16     p
5.05e-01    -3.5e+01 #  -4.9e-14     1.732e-26     f

Best performance (actual): 0.503

\(\gamma \approx 0.5\) means that the \(\mathcal{H}_\infty\) synthesis generated a controller \(K(s)\) that stabilizes the closed-loop system, and such that:

This is indeed what we can see by comparing \(|S|\) and \(|W_S|\) in Figure 29.

As was shown in Section 4, depending on the specifications, up to four closed-loop transfer function may be shaped (the Gang of four). This was summarized in Table 7.

For instance to limit the control input \(u\), \(KS\) should be shaped while to filter measurement noise, \(T\) should be shaped.

When multiple closed-loop transfer function are shaped at the same time, it is refereed to as “Mixed-Sensitivity \(\mathcal{H}_\infty\) Control” and is the subject of Section 6.

Depending on the closed-loop transfer function being shaped, different general control configuration are used and are described below.

P = [W1 -G*W1
     0   W2
     1  -G];

P = [W1 -G*W1
     0   G*W2
     1   -G];

P = [W1   -W1
     G*W2 -G*W2
     G    -G];

P = [W1 -G*W1
     0   W2
     0   G*W3
     1   -G];

P = [W1   -W1
     G*W2 -G*W2
     0     W3
     G    -G];

P = [ W1  -W1*G*W3 -G*W1
      0    0        W2
      1   -G*W3    -G];

When shaping multiple closed-loop transfer functions, one should be verify careful about the three following points that are further discussed:

• The shaped closed-loop transfer functions are linked by mathematical relations and cannot be shaped
• Closed-loop transfer function can only be shaped in certain frequency range.
• The size of the obtained controller may be very large and not implementable in practice

The control bandwidth is clearly limited by physical constrains such as sampling frequency, electronics bandwidth,

Similar relationship can be found for \(T\), \(KS\) and \(GS\).

Depending on the frequency band, the norms of the closed-loop transfer functions depend on the controller \(K\) and therefore can be shaped. However, in some frequency bands, the norms do not depend on the controller and therefore cannot be shaped.

Therefore the weighting functions should only focus on certainty frequency range depending on the transfer function being shaped. These regions are summarized in Figure 36.

• [ ] Control Diagram

• Inputs Signals
  • Reference steps of 1mm
  • Measurement noise is considered to be a white noise with a power spectral density of …
• Specifications
  • Follow reference steps with a response time of 0.1s and static error less than \(1 \mu m\)
  • Maximum control signal of 10N
  • Robustness
  • Reduce the effect of measurement noise on the position

k = 1e6; % Stiffness [N/m]
c = 4e2; % Damping [N/(m/s)]
m = 10; % Mass [kg]

G = 1/(m*s^2 + c*s + k); % Plant
Gd = (c*s + k)/(m*s^2 + c*s + k); % Disturbance

% Generate Input Signals
t = 0:1e-3:1;

r = zeros(size(t));
r(t>0.1) = 1e-3;


Fs = 1e3; % Sampling Frequency [Hz]
Ts = 1/Fs; % Sampling Time [s]
n = sqrt(Fs/2)*randn(1, length(t)); % Signal with an ASD equal to one
n = n*1e-6;

d = zeros(size(t));
d(t>0.5) = 5e-4;

• \(W_1\) is used to shape \(S\)
• \(W_2\) is used to shape \(KS\)
• \(W_3\) is used to shape \(T\)

P = [1 -G
     0  1
     0  G
     1 -G];

We start with the shaping of \(S\) alone. To not constrain \(KS\) and \(T\), we set small values for \(W_2\) and \(W_3\)

W2 = tf(1e-8);
W3 = tf(0.1);

W1 = generateWeight('G0', 1e3, ...
                    'G1', 1/2, ...
                    'Gc', sqrt(2), 'wc', 2*pi*10, ...
                    'n', 1);

Pw = blkdiag(W1, W2, W3, 1)*P;

K1 = hinfsyn(Pw, 1, 1, 'Display', 'on');

K1 = hinfsyn(Pw, 1, 1, 'Display', 'on');

  Test bounds:  0.5 <=  gamma  <=  0.552

   gamma        X>=0        Y>=0       rho(XY)<1    p/f
  5.25e-01     0.0e+00     0.0e+00     6.061e-16     p
  5.13e-01     0.0e+00     0.0e+00     0.000e+00     p
  5.06e-01    -2.5e+00 #  -6.4e-14     1.440e-15     f
  5.09e-01    -5.5e+00 #  -4.1e-14     9.510e-16     f
  Limiting gains...
  5.14e-01     0.0e+00     0.0e+00     1.039e-25     p
  5.14e-01     0.0e+00     0.0e+00     1.040e-25     p

  Best performance (actual): 0.514
'org_babel_eoe'
ans =
    'org_babel_eoe'

Z1 = lft(P, K1);

% Create model with r and n as inputs and y and u as outputs
Psim = [0  0  Gd  G
        0  0  0   1
        1 -1 -Gd -G];
Z1sim = lft(Psim, K1);

z = lsim(Z1sim, [r;n;d], t);
y1 = z(:,1);
u1 = z(:,2);

figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');

nexttile;
hold on;
plot(t, y1);
hold off;
xlabel('Time [s]'); ylabel('Output $y$ [m]');

nexttile;
hold on;
plot(t, u1);
hold off;
xlabel('Time [s]'); ylabel('Control Input $u$ [N]');

W2 = generateWeight('G0', 5e-7, ...
                    'G1', 1e-3, ...
                    'Gc', 1e-6, 'wc', 2*pi*100, ...
                    'n', 1);

Pw = blkdiag(W1, W2, W3, 1)*P;

K2 = hinfsyn(Pw, 1, 1, 'Display', 'on');

K1 = hinfsyn(Pw, 1, 1, 'Display', 'on');

  Test bounds:  0.51 <=  gamma  <=  1.2

   gamma        X>=0        Y>=0       rho(XY)<1    p/f
  7.82e-01     0.0e+00     0.0e+00     0.000e+00     p
  6.32e-01    -2.2e+00 #  -3.0e-14     -1.136e-32     f
  7.03e-01    -5.1e+00 #  -2.9e-14     5.128e-17     f
  7.41e-01    -1.1e+01 #  -1.0e-14     1.431e-22     f
  7.61e-01    -2.6e+01 #  -9.1e-15     1.215e-21     f
  7.72e-01    -6.9e+01 #  -1.9e-14     2.828e-16     f
  7.77e-01    -3.6e+02 #  -1.2e-14     1.213e-15     f

  Best performance (actual): 0.782
'org_babel_eoe'
ans =
    'org_babel_eoe'

Z2 = lft(P, K2);

Z2sim = lft(Psim, K2);

z = lsim(Z2sim, [r;n;d], t);
y2 = z(:,1);
u2 = z(:,2);

figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');

nexttile;
hold on;
plot(t, y1);
plot(t, y2);
hold off;
xlabel('Time [s]'); ylabel('Output $y$ [m]');

nexttile;
hold on;
plot(t, u1);
plot(t, u2);
hold off;
xlabel('Time [s]'); ylabel('Control Input $u$ [N]');

W3 = generateWeight('G0', 1e-1, ...
                    'G1', 1e4, ...
                    'Gc', 1, 'wc', 2*pi*30, ...
                    'n', 3);

Pw = blkdiag(W1, W2, W3, 1)*P;

K3 = hinfsyn(Pw, 1, 1, 'Display', 'on');

K3 = hinfsyn(Pw, 1, 1, 'Display', 'on');

  Test bounds:  0.578 <=  gamma  <=  1.66

   gamma        X>=0        Y>=0       rho(XY)<1    p/f
  9.78e-01    -1.3e+01 #  -6.2e-15     1.141e-18     f
  1.27e+00     0.0e+00     0.0e+00     1.524e-15     p
  1.12e+00     0.0e+00     0.0e+00     4.481e-15     p
  1.04e+00    -4.0e+01 #  -1.0e-13     1.102e-42     f
  1.08e+00    -6.6e+02 #  -4.4e-15     3.641e-18     f
  1.10e+00     0.0e+00     0.0e+00     6.052e-21     p
  1.09e+00     0.0e+00     0.0e+00     5.005e-15     p

  Best performance (actual): 1.09
'org_babel_eoe'
ans =
    'org_babel_eoe'

Z3 = lft(P, K3);

Z3sim = lft(Psim, K3);

z = lsim(Z3sim, [r;n;d], t);
y3 = z(:,1);
u3 = z(:,2);

figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');

nexttile;
hold on;
plot(t, y1);
plot(t, y2);
plot(t, y3);
hold off;
xlabel('Time [s]'); ylabel('Output $y$ [m]');

nexttile;
hold on;
plot(t, u1);
plot(t, u2);
plot(t, u3);
hold off;
xlabel('Time [s]'); ylabel('Control Input $u$ [N]');