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

The typical methodology when applying Model Based Control to a plant is schematically shown in Figure 1. It consists of three steps:

1. Identification or modeling: \(\Longrightarrow\) mathematical model
2. Translate the specifications into mathematical criteria:
   • Specifications: Response Time, Noise Rejection, Maximum input amplitude, Robustness, …
   • Mathematical Criteria: Cost Function, Shape of TF
3. Synthesis: research of \(K\) that satisfies the specifications for the model of the system

In this document, we will mainly focus on steps 2 and 3.

Step 2 will be discussed in Section 4. There are two main methods for the controller synthesis (step 3):

• open loop shaping discussed in Section 2
• closed loop shaping discussed in Sections 4 and 5

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.

Hi Musa, Thank you very much for sharing this awesome package. For a long time, I am dreaming of being abble to export source blocks to HTML tha are surounded by <details> blocks.

For now, I am manually adding #+HTML: <details><summary>Code</summary> and #+HTML: </details> around the source blocks I want to hide… This is a very simple solution, but not so elegent nor practical.

Do you have any idea if it would be easy to extend to org-mode export of source blocks to add such functionallity?

Similarly, I would love to be able to export a <span> block with the name of the file corresponding to the source block. For instance, if a particular source block is tangled to script.sh, it would be so nice to display the filename when exporting!

Thanks in advance

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);

The Bode plot of the obtained controller is shown in Figure 10.

Let’s briefly analyze this controller:

• 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.

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

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

Using Matlab, the generalized plant can be defined as follows:

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

As shown in Section 2, the loop gain \(L(s) = G(s) K(s)\) is a useful and easy tool for the manual design of controllers.

\(L(s)\) is very easy to shape as it depends linearly on \(K(s)\). Moreover, important quantities such as the stability margins and the control bandwidth can be estimated from the shape/phase of \(L(s)\).

However, the loop gain \(L(s)\) does not directly give the performances of the closed-loop system. The closed loop system behavior is indeed determined by the closed-loop transfer functions.

If we consider the feedback system shown in Figure 17, we can link to the following specifications to closed-loop transfer functions. This is summarized in Table 5.

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 want 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 want to act on the complementary sensitivity function \(T(s)\).

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;
Gr = 1/k*(s/w0/4 + 1)/(s^2/w0^2 + 2*xi*s/w0 + 1);

isstable(feedback(Gp,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 away from the \(-1\) point (called the unstable point) in the complex plane. From this, we define the module margin.

A typical Nyquist plot is shown in Figure 23. The gain, phase and module margins are graphically shown to have an idea of what they represent.

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.

• [ ] Maybe put this section in Previous chapter

Let’s say we want to shape the sensitivity transfer function corresponding to the transfer function from \(r\) to \(\epsilon\) of the control architecture shown in Figure 24.

If the \(\mathcal{H}_\infty\) synthesis is directly applied on the generalized plant \(P(s)\) shown in Figure 24, if will minimize the \(\mathcal{H}_\infty\) norm of transfer function from \(r\) to \(\epsilon\) (the sensitivity transfer function).

However, as the \(\mathcal{H}_\infty\) norm is the maximum peak value of the transfer function’s magnitude, it does not allow to shape the norm over all frequencies.

A trick is to include a weighting function in the generalized plant as shown in Figure 25. Applying the \(\mathcal{H}_\infty\) synthesis to the weighted generalized plant \(\tilde{P}(s)\) (Figure 25) will generate a controller \(K(s)\) that minimizes the \(\mathcal{H}_\infty\) norm between \(r\) and \(\tilde{\epsilon}\):

Once the weighting function is designed, it should be added to the generalized plant as shown in Figure 25.

The weighted generalized plant can be defined in Matlab by either re-defining all the inputs or by pre-multiplying the (non-weighted) generalized plant by a block-diagonal MIMO transfer function containing the weights for the outputs \(z\) and 1 for the outputs \(v\).

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

% Alternative
Pw = blkdiag(Ws, 1)*P;

Weighting function used must be proper, stable and minimum phase transfer functions.

proper

more poles than zeros, this implies \(\lim_{\omega \to \infty} |W(j\omega)| < \infty\)

stable

no poles in the right half plane

minimum phase

no zeros in the right half plane

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

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

with:

• dcgain
• freq
• mag
• hfgain

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);

• Robustness: Module margin > 2 (\(\Rightarrow \text{GM} > 2 \text{ and } \text{PM} > 29^o\))
• Bandwidth:
• Slope of -2

First, the weighting functions is generated.

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

It is then added to the generalized plant.

Pw = blkdiag(Ws, 1)*P;

And the \(\mathcal{H}_\infty\) synthesis is performed.

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

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     4.497e-28     p
  Limiting gains...
  5.05e-01     0.0e+00     0.0e+00     0.000e+00     p
  5.05e-01    -1.8e+01 #  -2.9e-15     1.514e-15     f

  Best performance (actual): 0.504

The obtained \(\gamma \approx 0.5\) means that it found a controller \(K(s)\) that stabilize the closed-loop system, and such that:

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