diff --git a/figs/general_conf_shaping_S_GS.pdf b/figs/general_conf_shaping_S_GS.pdf new file mode 100644 index 0000000..5f5333c Binary files /dev/null and b/figs/general_conf_shaping_S_GS.pdf differ diff --git a/figs/general_conf_shaping_S_GS.png b/figs/general_conf_shaping_S_GS.png new file mode 100644 index 0000000..58ee0cf Binary files /dev/null and b/figs/general_conf_shaping_S_GS.png differ diff --git a/figs/general_conf_shaping_S_GS.svg b/figs/general_conf_shaping_S_GS.svg new file mode 100644 index 0000000..22992e6 --- /dev/null +++ b/figs/general_conf_shaping_S_GS.svg @@ -0,0 +1,228 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/figs/general_conf_shaping_S_T_GS.pdf b/figs/general_conf_shaping_S_T_GS.pdf new file mode 100644 index 0000000..3e97891 Binary files /dev/null and b/figs/general_conf_shaping_S_T_GS.pdf differ diff --git a/figs/general_conf_shaping_S_T_GS.png b/figs/general_conf_shaping_S_T_GS.png new file mode 100644 index 0000000..b79634d Binary files /dev/null and b/figs/general_conf_shaping_S_T_GS.png differ diff --git a/figs/general_conf_shaping_S_T_GS.svg b/figs/general_conf_shaping_S_T_GS.svg new file mode 100644 index 0000000..dc047d4 --- /dev/null +++ b/figs/general_conf_shaping_S_T_GS.svg @@ -0,0 +1,261 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/figs/mixed_sensitivity_control_schematic.pdf b/figs/mixed_sensitivity_control_schematic.pdf new file mode 100644 index 0000000..0f02105 Binary files /dev/null and b/figs/mixed_sensitivity_control_schematic.pdf differ diff --git a/figs/mixed_sensitivity_control_schematic.png b/figs/mixed_sensitivity_control_schematic.png new file mode 100644 index 0000000..5a27bd0 Binary files /dev/null and b/figs/mixed_sensitivity_control_schematic.png differ diff --git a/figs/mixed_sensitivity_control_schematic.svg b/figs/mixed_sensitivity_control_schematic.svg new file mode 100644 index 0000000..9d1950e --- /dev/null +++ b/figs/mixed_sensitivity_control_schematic.svg @@ -0,0 +1,255 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/figs/mixed_sensitivity_ref_tracking.pdf b/figs/mixed_sensitivity_ref_tracking.pdf index ad0fb86..0cfc285 100644 Binary files a/figs/mixed_sensitivity_ref_tracking.pdf and b/figs/mixed_sensitivity_ref_tracking.pdf differ diff --git a/figs/mixed_sensitivity_ref_tracking.png b/figs/mixed_sensitivity_ref_tracking.png index cabe270..efc2a95 100644 Binary files a/figs/mixed_sensitivity_ref_tracking.png and b/figs/mixed_sensitivity_ref_tracking.png differ diff --git a/figs/mixed_sensitivity_ref_tracking.svg b/figs/mixed_sensitivity_ref_tracking.svg index d424b5a..d2e78cd 100644 --- a/figs/mixed_sensitivity_ref_tracking.svg +++ b/figs/mixed_sensitivity_ref_tracking.svg @@ -1,56 +1,50 @@ - + - + - - - - - - - - - - - - - - - - - - - + - + - + - + + + + - + + + + - - + + - + + + + + + + @@ -60,22 +54,22 @@ - + - + - + - + - + @@ -83,6 +77,9 @@ + + + @@ -95,149 +92,126 @@ - - - - - - - - - - + - + - + - + + + + + + + + + + + - - - - + + + + - - - + - - - - - - - - - - - + - - - - + - + - + - + - - - - - + - + - - - - - + - + - + - + - + - - + - - - - - - - - - - - - - + - - - - - - + - + - - + + + + + + + + + + + + + + + + + + + - + - + - - + + - + - - + + - + - + - + - + - + diff --git a/index.html b/index.html index 1d02428..67d77ce 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + A brief and practical introduction to \(\mathcal{H}_\infty\) Control @@ -30,95 +30,111 @@

Table of Contents

-This document is structured as follows: +The purpose of this document is to give a practical introduction to the wonderful world of \(\mathcal{H}_\infty\) Control. +

+ +

+No attend is made to provide an exhaustive treatment of the subject. +\(\mathcal{H}_\infty\) Control is a very broad topic and entire books are written on it. +Therefore, for more advanced discussion, please have a look at the recommended references at the bottom of this document. +

+ +

+When possible, Matlab scripts used for the example/exercises are provided such that you can easily test them on your computer. +

+ + +

+The general structure of this document is as follows:

-
-

1 Introduction to Model Based Control

+
+

1 Introduction to Model Based Control

- +

-
-

1.1 Model Based Control - Methodology

+
+

1.1 Model Based Control - Methodology

- +

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

    @@ -132,7 +148,7 @@ It consists of three steps:
-
+

control-procedure.png

Figure 1: Typical Methodoly for Model Based Control

@@ -143,24 +159,24 @@ In this document, we will mainly focus on steps 2 and 3.

-Step 2 will be discussed in Section 4. +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 6
  • +
  • open loop shaping discussed in Section 2
  • +
  • closed loop shaping discussed in Sections 4 and 6
-
-

1.2 From Classical Control to Robust Control

+
+

1.2 From Classical Control to Robust Control

- +

- +
@@ -293,7 +309,7 @@ There are two main methods for the controller synthesis (step 3):
Table 1: Table summurazing the main differences between classical, modern and robust control
-
+

robustness_performance.png

Figure 2: Comparison of the performance and robustness of classical control methods, modern control methods and robust control methods. The required information on the plant to succesfuly apply each of the control methods are indicated by the colors.

@@ -301,27 +317,27 @@ There are two main methods for the controller synthesis (step 3):
-
-

1.3 Example System

+
+

1.3 Example System

- +

-Let’s consider the model shown in Figure 3. +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. +The notations used are listed in Table 2.

-
+

mech_sys_1dof_inertial_contr.png

Figure 3: Test System consisting of a payload with a mass \(m\) on top of an active system with a stiffness \(k\), damping \(c\) and an actuator. A feedback controller \(K(s)\) is added to position / isolate the payload.

- +
@@ -407,7 +423,7 @@ The notations used are listed in Table 2.
Table 2: Example system variables
-
+

Derive the following open-loop transfer functions:

@@ -438,37 +454,13 @@ You can follow this generic procedure:
-

-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? +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.

-

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

- - -
+

classical_feedback_test_system.png

Figure 4: Block diagram corresponding to the example system

@@ -486,7 +478,7 @@ Let’s define the system parameters on Matlab.

-And now the system dynamics \(G(s)\) and \(G_d(s)\) (their bode plots are shown in Figures 5 and 6). +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
@@ -495,14 +487,14 @@ And now the system dynamics \(G(s)\) and \(G_d(s)\) (their bode plots are shown
 
-
+

bode_plot_example_afm.png

Figure 5: Bode plot of the plant \(G(s)\)

-
+

bode_plot_example_Gd.png

Figure 6: Magnitude of the disturbance transfer function \(G_d(s)\)

@@ -511,44 +503,44 @@ And now the system dynamics \(G(s)\) and \(G_d(s)\) (their bode plots are shown
-
-

2 Classical Open Loop Shaping

+
+

2 Classical Open Loop Shaping

- +

-
-

2.1 Introduction to Loop Shaping

+
+

2.1 Introduction to Loop Shaping

- +

-
+

Loop Shaping refers to a design procedure that involves explicitly shaping the magnitude of the Loop Transfer Function \(L(s)\).

-
+

-The Loop Gain \(L(s)\) usually refers to as the product of the controller and the plant (“Gain around the loop”, see Figure 7): +The Loop Gain \(L(s)\) usually refers to as the product of the controller and the plant (“Gain around the loop”, see Figure 7):

\begin{equation} L(s) = G(s) \cdot K(s) \label{eq:loop_gain} \end{equation} -
+

open_loop_shaping.png

Figure 7: Classical Feedback Architecture

@@ -573,11 +565,11 @@ The Open Loop shape is usually done manually has the loop gain \(L(s)\) depends

-\(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). +\(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).

-
+

open_loop_shaping_shape.png

Figure 8: Typical Wanted Shape for the Loop Gain \(L(s)\)

@@ -585,14 +577,14 @@ The Open Loop shape is usually done manually has the loop gain \(L(s)\) depends
-
-

2.2 Example of Open Loop Shaping

+
+

2.2 Example of Open Loop Shaping

- +

-
+

Let’s take our example system and try to apply the Open-Loop shaping strategy to design a controller that fulfils the following specifications:

@@ -604,7 +596,7 @@ Let’s take our example system and try to apply the Open-Loop shaping strat
-
+

Using SISOTOOL, design a controller that fulfill the specifications.

@@ -621,7 +613,7 @@ In order to have the wanted Roll-off, two integrators are used, a lead is also a

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

K = 14e8 * ... % Gain
@@ -632,7 +624,7 @@ The obtained controller is shown below, and the bode plot of the Loop Gain is sh
 
-
+

loop_gain_manual_afm.png

Figure 9: Bode Plot of the obtained Loop Gain \(L(s) = G(s) K(s)\)

@@ -680,11 +672,11 @@ And we can verify that we have the wanted stability margins:
-
-

2.3 \(\mathcal{H}_\infty\) Loop Shaping Synthesis

+
+

2.3 \(\mathcal{H}_\infty\) Loop Shaping Synthesis

- +

@@ -711,7 +703,7 @@ where:

  • K is the synthesize controller
  • -
    +

    Matlab documentation of loopsyn (link).

    @@ -720,11 +712,11 @@ Matlab documentation of loopsyn ( -

    2.4 Example of the \(\mathcal{H}_\infty\) Loop Shaping Synthesis

    +
    +

    2.4 Example of the \(\mathcal{H}_\infty\) Loop Shaping Synthesis

    - +

    @@ -755,7 +747,7 @@ The \(\mathcal{H}_\infty\) optimal open loop shaping synthesis is performed usin

    -
    +

    It is always important to analyze the controller after the synthesis is performed.

    @@ -767,7 +759,7 @@ In the end, a synthesize controller is just a combination of low pass filters, h

    -Let’s briefly analyze the obtained controller which bode plot is shown in Figure 10: +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
    • @@ -776,28 +768,28 @@ Let’s briefly analyze the obtained controller which bode plot is shown in
    -
    +

    open_loop_shaping_hinf_K.png

    Figure 10: Obtained controller \(K\) using the open-loop \(\mathcal{H}_\infty\) shaping

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

    -
    +

    open_loop_shaping_hinf_L.png

    Figure 11: Obtained Open Loop Gain \(L(s) = G(s) K(s)\) and comparison with the wanted Loop gain \(L_w\)

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

    - +
    @@ -838,29 +830,29 @@ Let’s now compare the obtained stability margins of the \(\mathcal{H}_\inf -
    -

    3 A first Step into the \(\mathcal{H}_\infty\) world

    +
    +

    3 A first Step into the \(\mathcal{H}_\infty\) world

    - +

    -
    -

    3.1 The \(\mathcal{H}_\infty\) Norm

    +
    +

    3.1 The \(\mathcal{H}_\infty\) Norm

    - +

    -
    +

    The \(\mathcal{H}_\infty\) norm is defined as the peak of the maximum singular value of the frequency response

    @@ -877,7 +869,7 @@ For a SISO system \(G(s)\), it is simply the peak value of \(|G(j\omega)|\) as a
    -
    +

    Let’s compute the \(\mathcal{H}_\infty\) norm of our test plant \(G(s)\) using the hinfnorm function:

    @@ -892,11 +884,11 @@ Let’s compute the \(\mathcal{H}_\infty\) norm of our test plant \(G(s)\) u

    -We can see that the \(\mathcal{H}_\infty\) norm of \(G(s)\) does corresponds to the peak value of \(|G(j\omega)|\) as a function of frequency as shown in Figure 12. +We can see that the \(\mathcal{H}_\infty\) norm of \(G(s)\) does corresponds to the peak value of \(|G(j\omega)|\) as a function of frequency as shown in Figure 12.

    -
    +

    hinfinity_norm_siso_bode.png

    Figure 12: Example of the \(\mathcal{H}_\infty\) norm of a SISO system

    @@ -906,14 +898,14 @@ We can see that the \(\mathcal{H}_\infty\) norm of \(G(s)\) does corresponds to
    -
    -

    3.2 \(\mathcal{H}_\infty\) Synthesis

    +
    +

    3.2 \(\mathcal{H}_\infty\) Synthesis

    - +

    -
    +

    \(\mathcal{H}_\infty\) synthesis is a method that uses an algorithm (LMI optimization, Riccati equation) to find a controller that stabilize the system and that minimizes the \(\mathcal{H}_\infty\) norms of defined transfer functions.

    @@ -949,11 +941,11 @@ Note that there are many ways to use the \(\mathcal{H}_\infty\) Synthesis:
    -
    -

    3.3 The Generalized Plant

    +
    +

    3.3 The Generalized Plant

    - +

    @@ -963,9 +955,9 @@ It makes things much easier for the following steps.

    -The generalized plant, usually noted \(P(s)\), is shown in Figure 13. +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. +The meaning of the inputs and outputs are summarized in Table 4.

    @@ -978,14 +970,14 @@ It can indeed represent feedback as well as feedforward control architectures. \end{equation} -

    +

    general_plant.png

    Figure 13: Inputs and Outputs of the generalized Plant

    -
    -
    Table 3: Comparison of the characteristics obtained with the two methods
    +
    +
    @@ -1031,18 +1023,18 @@ It can indeed represent feedback as well as feedforward control architectures. -
    -

    3.4 The \(\mathcal{H}_\infty\) Synthesis applied on the Generalized plant

    +
    +

    3.4 The \(\mathcal{H}_\infty\) Synthesis applied on the Generalized plant

    - +

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

    -
    +
    \(\mathcal{H}_\infty\) Synthesis applied on the generalized plant
    @@ -1057,7 +1049,7 @@ After \(K\) is found, the system is robustified by adjusting the response
    -
    +

    general_control_names.png

    Figure 14: General Control Configuration

    @@ -1089,28 +1081,29 @@ where:
    -
    -

    3.5 From a Classical Feedback Architecture to a Generalized Plant

    +
    +

    3.5 From a Classical Feedback Architecture to a Generalized Plant

    - +

    -The procedure to convert a typical control architecture as the one shown in Figure 15 to a generalized Plant is as follows: +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. -
    3. Remove \(K\) and rearrange the inputs and outputs to match the generalized configuration
    4. +
    5. Remove \(K\) and rearrange the inputs and outputs to match the generalized configuration shown in Figure 13
    -
    -

    -Compute the Generalized plant of corresponding to the tracking control architecture shown in Figure 15 -

    +
    +
      +
    1. Convert the tracking control architecture shown in Figure 15 to a generalized configuration
    2. +
    3. Compute the transfer function matrix using Matlab as a function or \(K\) and \(G\)
    4. +
    -
    +

    classical_feedback_tracking.png

    Figure 15: Classical Feedback Control Architecture (Tracking)

    @@ -1122,29 +1115,27 @@ First, define the signals of the generalized plant:

    • Exogenous inputs: \(w = r\)
    • -
    • Signals to be minimized: \(z_1 = \epsilon\), \(z_2 = u\)
    • -
    • Control signals: \(v = y\)
    • -
    • Control inputs: \(u\)
    • +
    • Signals to be minimized: +Usually, we want to minimize the tracking errors \(\epsilon\) and the control signal \(u\): \(z = [\epsilon,\ u]\)
    • +
    • Controller inputs: this is the signal at the input of the controller: \(v = \epsilon\)
    • +
    • Control inputs: signal generated by the controller: \(u\)

    -Then, Remove \(K\) and rearrange the inputs and outputs. +Then, Remove \(K\) and rearrange the inputs and outputs as in Figure 13.

    Answer

    -The obtained generalized plant shown in Figure 16. +The obtained generalized plant shown in Figure 16.

    -
    +

    mixed_sensitivity_ref_tracking.png

    Figure 16: Generalized plant of the Classical Feedback Control Architecture (Tracking)

    -
    -
    -

    @@ -1157,26 +1148,29 @@ Using Matlab, the generalized plant can be defined as follows: P.InputName = {'w', 'u'}; P.OutputName = {'e', 'u', 'v'}; +

    + +
    -
    -

    4 Modern Interpretation of Control Specifications

    +
    +

    4 Modern Interpretation of Control Specifications

    - +

    -As shown in Section 2, the loop gain \(L(s) = G(s) K(s)\) is a useful and easy tool when manually designing controllers. +As shown in Section 2, the loop gain \(L(s) = G(s) K(s)\) is a useful and easy tool when manually designing controllers. This is mainly due to the fact that \(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)\).

    @@ -1186,11 +1180,11 @@ However, the loop gain \(L(s)\) does not directly give the performances o

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

    -
    Table 4: Notations for the general configuration
    +
    @@ -1232,33 +1226,33 @@ This is summarized in Table 5. - +
    Table 5: Typical Specification and associated closed-loop transfer function
    Robustness (stability margins)Module margin (see Section 4.3)Module margin (see Section 4.3)
    -
    +

    gang_of_four_feedback.png

    Figure 17: Simple Feedback Architecture

    -
    -

    4.1 Closed Loop Transfer Functions

    +
    +

    4.1 Closed Loop Transfer Functions

    - +

    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)\).

    -
    +

    -Write the output signals \([\epsilon, u, y]\) as a function of the systems \(K(s), G(s)\) and of the input signals \([r, d, n]\) as shown in Figure 17. +Write the output signals \([\epsilon, u, y]\) as a function of the systems \(K(s), G(s)\) and of the input signals \([r, d, n]\) as shown in Figure 17.

    Hint @@ -1295,9 +1289,9 @@ The following equations should be obtained:
    -
    +

    -We can see that they are 4 different transfer functions describing the behavior of the system in Figure 17. +We can see that they are 4 different transfer functions describing the behavior of the system in Figure 17. These called the Gang of Four:

    \begin{align} @@ -1309,7 +1303,7 @@ These called the Gang of Four:
    -
    +

    If a feedforward controller is included, a Gang of Six transfer functions can be defined. More on that in this short video. @@ -1333,11 +1327,11 @@ Similarly, to reduce the effect of measurement noise \(n\) on the output \(y\),

    -
    -

    4.2 Sensitivity Function

    +
    +

    4.2 Sensitivity Function

    - +

    @@ -1348,14 +1342,14 @@ In this section, we will see how the shape of the sensitivity function will impa

    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 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. +The comparison of the sensitivity functions shapes and their effect on the step response is summarized in Table 6.

    - +
    @@ -1394,20 +1388,20 @@ The comparison of the sensitivity functions shapes and their effect on the step
    Table 6: Comparison of the sensitivity function shape and the corresponding step response for the three controller variations
    -
    +

    sensitivity_shape_effect.png

    Figure 18: Sensitivity function magnitude \(|S(j\omega)|\) corresponding to the reference controller \(K_r(s)\) and the three modified controllers \(K_i(s)\)

    -
    +

    sensitivity_shape_effect_step.png

    Figure 19: Step response (response from \(r\) to \(y\)) for the different controllers

    -
    +
    Closed-Loop Bandwidth

    The closed-loop bandwidth \(\omega_b\) is the frequency where \(|S(j\omega)|\) first crosses \(1/\sqrt{2} = -3dB\) from below. @@ -1420,9 +1414,9 @@ In general, a large bandwidth corresponds to a faster rise time.

    -
    +

    -From the simple analysis above, we can draw a first estimation of the wanted shape for the sensitivity function in Figure 20. +From the simple analysis above, we can draw a first estimation of the wanted shape for the sensitivity function in Figure 20.

    @@ -1437,7 +1431,7 @@ This will become clear in the next section about the module margin. -

    +

    h-infinity-spec-S.png

    Figure 20: Typical wanted shape of the Sensitivity transfer function

    @@ -1447,18 +1441,18 @@ This will become clear in the next section about the module margin.
    -
    -

    4.3 Robustness: Module Margin

    +
    +

    4.3 Robustness: Module Margin

    - +

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

    -
    +

    Let’s consider the following plant \(G_t(s)\):

    @@ -1472,7 +1466,7 @@ Gt = 1/k*(s

    -Let’s say we have designed a controller \(K_t(s)\) that gives the loop gain shown in Figure 21. +Let’s say we have designed a controller \(K_t(s)\) that gives the loop gain shown in Figure 21.

    @@ -1481,7 +1475,7 @@ Let’s say we have designed a controller \(K_t(s)\) that gives the loop gai

    -The following characteristics can be determined from Figure 21: +The following characteristics can be determined from Figure 21:

    • bandwidth of \(\approx 10\, \text{Hz}\)
    • @@ -1494,7 +1488,7 @@ This might indicate very good robustness properties of the closed-loop system.

      -
      +

      phase_gain_margin_model_plant.png

      Figure 21: Bode plot of the obtained Loop Gain \(L(s)\)

      @@ -1510,7 +1504,7 @@ Now let’s suppose the “real” plant \(G_r(s)\) as a slightly lo

      -The obtained “real” loop gain is shown in Figure 22. +The obtained “real” loop gain is shown in Figure 22. At a frequency little bit above 100Hz, the phase of the loop gain reaches -180 degrees while its magnitude is more than one which indicated instability.

      @@ -1528,7 +1522,7 @@ It is confirmed by checking the stability of the closed loop system: -
      +

      phase_gain_margin_real_plant.png

      Figure 22: Bode plots of \(L(s)\) (loop gain corresponding the nominal plant) and \(L_r(s)\) (loop gain corresponding to the real plant)

      @@ -1548,7 +1542,7 @@ This is due to the fact that the gain and phase margin are robustness indicators Let’s now determine a new robustness indicator based on the Nyquist Stability Criteria.

      -
      +
      Nyquist Stability Criteria (for stable systems)
      If the open-loop transfer function \(L(s)\) is stable, then the closed-loop system is unstable for any encirclement of the point \(−1\) on the Nyquist plot.
      @@ -1557,7 +1551,7 @@ Let’s now determine a new robustness indicator based on the Nyquist Stabil
      -
      +

      For more information about the general Nyquist Stability Criteria, you may want to look at this video.

      @@ -1569,21 +1563,21 @@ From the Nyquist stability criteria, it is clear that we want \(L(j\omega)\) to This minimum distance is called the module margin.

      -
      +
      Module Margin
      The Module Margin \(\Delta M\) is defined as the minimum distance between the point \(-1\) and the loop gain \(L(j\omega)\) in the complex plane.
      -
      +

      -A typical Nyquist plot is shown in Figure 23. +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.

      -
      +

      module_margin_example.png

      Figure 23: Nyquist plot with visual indication of the Gain margin \(\Delta G\), Phase margin \(\Delta \phi\) and Module margin \(\Delta M\)

      @@ -1592,7 +1586,7 @@ The gain, phase and module margins are graphically shown to have an idea of what

      -As expected from Figure 23, there is a close relationship between the module margin and the gain and phase margins. +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:

      \begin{equation} @@ -1611,7 +1605,7 @@ Let’s now try to express the Module margin \(\Delta M\) as an \(\mathcal{H &= \frac{1}{\|S\|_\infty} \end{align*} -
      +

      The \(\mathcal{H}_\infty\) norm of the sensitivity function \(\|S\|_\infty\) is a measure of the Module margin \(\Delta M\) and therefore an indicator of the system robustness.

      @@ -1631,14 +1625,14 @@ Note that this is why large peak value of \(|S(j\omega)|\) usually indicate robu And we know understand why setting an upper bound on the magnitude of \(S\) is generally a good idea.

      -
      +

      Typical, we require \(\|S\|_\infty < 2 (6dB)\) which implies \(\Delta G \ge 2\) and \(\Delta \phi \ge 29^o\)

      -
      +

      To learn more about module/disk margin, you can check out this video.

      @@ -1647,11 +1641,11 @@ To learn more about module/disk margin, you can check out -

      4.4 Other Requirements

      +
      +

      4.4 Other Requirements

      - +

      @@ -1680,7 +1674,7 @@ with: Noise Attenuation: typical wanted shape for \(T\)

      - +
      @@ -1733,11 +1727,11 @@ Noise Attenuation: typical wanted shape for \(T\) -
      -

      5 \(\mathcal{H}_\infty\) Shaping of closed-loop transfer functions

      +
      +

      5 \(\mathcal{H}_\infty\) Shaping of closed-loop transfer functions

      - +

      In the previous sections, we have seen that the performances of the system depends on the shape of the closed-loop transfer function. @@ -1754,22 +1748,22 @@ But don’t worry, the \(\mathcal{H}_\infty\) synthesis will do this job for

      This -Section 5.1 -Section 5.2 -Section 5.3 -Section 5.4 +Section 5.1 +Section 5.2 +Section 5.3 +Section 5.4

      -
      -

      5.1 How to Shape closed-loop transfer function? Using Weighting Functions!

      +
      +

      5.1 How to Shape closed-loop transfer function? Using Weighting Functions!

      - +

      -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. +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\).

      @@ -1779,16 +1773,16 @@ However, as the \(\mathcal{H}_\infty\) norm is the maximum peak value of the tra

      -
      +

      loop_shaping_S_without_W.png

      Figure 24: Generalized Plant

      -
      +

      -The trick is to include a weighting function \(W_S(s)\) in the generalized plant as shown in Figure 25. +The trick is to include a weighting function \(W_S(s)\) in the generalized plant as shown in Figure 25.

      @@ -1806,7 +1800,7 @@ Let’s now show how this is equivalent as shaping the sensitivity fu \Leftrightarrow & \left| S(j\omega) \right| < \frac{1}{\left| W_s(j\omega) \right|} \quad \forall \omega \label{eq:sensitivity_shaping} \end{align} -

      +

      As shown in Equation \eqref{eq:sensitivity_shaping}, the \(\mathcal{H}_\infty\) synthesis applying on the weighted generalized plant allows to shape the magnitude of the sensitivity transfer function.

      @@ -1818,15 +1812,15 @@ Therefore, the choice of the weighting function \(W_s(s)\) is very important: it
      -
      +

      loop_shaping_S_with_W.png

      Figure 25: Weighted Generalized Plant

      -
      +

      -Using matlab, compute the weighted generalized plant shown in Figure 26 as a function of \(G(s)\) and \(W_S(s)\). +Using matlab, compute the weighted generalized plant shown in Figure 26 as a function of \(G(s)\) and \(W_S(s)\).

      Hint @@ -1864,18 +1858,18 @@ The second solution is however more general, and can also be used when weights a
      -
      -

      5.2 Design of Weighting Functions

      +
      +

      5.2 Design of Weighting Functions

      - +

      Weighting function included in the generalized plant 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
      @@ -1900,9 +1894,9 @@ with:
    • hfgain: high frequency gain
    • -
      +

      -The Matlab code below produces a weighting function with the following characteristics (Figure 26): +The Matlab code below produces a weighting function with the following characteristics (Figure 26):

      • Low frequency gain of 100
      • @@ -1916,7 +1910,7 @@ The Matlab code below produces a weighting function with the following character
      -
      +

      first_order_weight.png

      Figure 26: Obtained Magnitude of the Weighting Function

      @@ -1924,7 +1918,7 @@ The Matlab code below produces a weighting function with the following character
      -
      +

      Quite often, higher orders weights are required.

      @@ -1956,7 +1950,7 @@ A Matlab function implementing Equation \eqref{eq:weight_formula_advanced} is sh

      -
      function [W] = generateWeight(args)
      +
      function [W] = generateWeight(args)
           arguments
               args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1
               args.G1 (1,1) double {mustBeNumeric, mustBePositive} = 10
      @@ -1991,11 +1985,11 @@ W3 = generateWeight('G0', 1e2, 27.
      +The obtained shapes are shown in Figure 27.
       

      -
      +

      high_order_weight.png

      Figure 27: Higher order weights using Equation \eqref{eq:weight_formula_advanced}

      @@ -2005,11 +1999,11 @@ The obtained shapes are shown in Figure 27.
      -
      -

      5.3 Shaping the Sensitivity Function

      +
      +

      5.3 Shaping the Sensitivity Function

      - +

      @@ -2022,17 +2016,17 @@ Let’s design a controller using the \(\mathcal{H}_\infty\) synthesis that

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

      -
      +

      Translate the requirements as upper bounds on the Sensitivity function and design the corresponding Weight using Matlab.

      Hint

      -The typical wanted upper bound of the sensitivity function is shown in Figure 28. +The typical wanted upper bound of the sensitivity function is shown in Figure 28.

      @@ -2050,7 +2044,7 @@ Remember that the wanted upper bound of the sensitivity function is defined by t

      -
      +

      h-infinity-spec-S.png

      Figure 28: Typical wanted shape of the Sensitivity transfer function

      @@ -2115,7 +2109,7 @@ And the \(\mathcal{H}_\infty\) synthesis is performed on the weighted gen
      -
      +
       Test bounds:  0.5 <=  gamma  <=  0.51
       
        gamma        X>=0        Y>=0       rho(XY)<1    p/f
      @@ -2136,10 +2130,10 @@ Best performance (actual): 0.503
       \end{aligned}
       
       

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

      -
      +

      Having \(\gamma < 1\) means that the \(\mathcal{H}_\infty\) synthesis found a controller such that the specified closed-loop transfer functions are bellow the specified upper bounds.

      @@ -2152,7 +2146,7 @@ It just means that at some frequency, one of the closed-loop transfer functions
      -
      +

      results_sensitivity_hinf.png

      Figure 29: Weighting function and obtained closed-loop sensitivity

      @@ -2160,16 +2154,16 @@ It just means that at some frequency, one of the closed-loop transfer functions
      -
      -

      5.4 Shaping multiple closed-loop transfer functions

      +
      +

      5.4 Shaping multiple closed-loop transfer functions

      - +

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

      @@ -2177,7 +2171,7 @@ For instance to limit the control input \(u\), \(KS\) should be shaped while to

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

      @@ -2185,14 +2179,14 @@ Depending on the closed-loop transfer function being shaped, different general c

      Shaping of S and KS -
      +

      general_conf_shaping_S_KS.png

      Figure 30: Generalized Plant to shape \(S\) and \(KS\)

      -
      P = [W1 -G*W1
      +
      P = [W1 -G*W1
            0   W2
            1  -G];
       
      @@ -2205,14 +2199,14 @@ Depending on the closed-loop transfer function being shaped, different general c
      Shaping of S and T -
      +

      general_conf_shaping_S_T.png

      Figure 31: Generalized Plant to shape \(S\) and \(T\)

      -
      P = [W1 -G*W1
      +
      P = [W1 -G*W1
            0   G*W2
            1   -G];
       
      @@ -2223,16 +2217,35 @@ Depending on the closed-loop transfer function being shaped, different general c
    • \(W_2\) is used to shape \(T\)
    • -
      Shaping of S, T and KS +
      Shaping of S and GS -
      -

      general_conf_shaping_S_T_KS.png +

      +

      general_conf_shaping_S_GS.png

      -

      Figure 32: Generalized Plant to shape \(S\), \(T\) and \(KS\)

      +

      Figure 32: Generalized Plant to shape \(S\) and \(GS\)

      -
      P = [W1 -G*W1
      +
      P = [W1   -W1
      +     G*W2 -G*W2
      +     G    -G];
      +
      +
      +
        +
      • \(W_1\) is used to shape \(S\)
      • +
      • \(W_2\) is used to shape \(GS\)
      • +
      +
      +
      Shaping of S, T and KS + +
      +

      general_conf_shaping_S_T_KS.png +

      +

      Figure 33: Generalized Plant to shape \(S\), \(T\) and \(KS\)

      +
      + +
      +
      P = [W1 -G*W1
            0   W2
            0   G*W3
            1   -G];
      @@ -2245,16 +2258,38 @@ Depending on the closed-loop transfer function being shaped, different general c
       
    • \(W_3\) is used to shape \(T\)
    • -
      Shaping of S, T, KS and GS +
      Shaping of S, T and GS -
      -

      general_conf_shaping_S_T_KS_GS.png +

      +

      general_conf_shaping_S_T_GS.png

      -

      Figure 33: Generalized Plant to shape \(S\), \(T\), \(KS\) and \(GS\)

      +

      Figure 34: Generalized Plant to shape \(S\), \(T\) and \(GS\)

      -
      P = [ W1  -W1*G*W3 -G*W1
      +
      P = [W1   -W1
      +     G*W2 -G*W2
      +     0     W3
      +     G    -G];
      +
      +
      + +
        +
      • \(W_1\) is used to shape \(S\)
      • +
      • \(W_2\) is used to shape \(GS\)
      • +
      • \(W_3\) is used to shape \(T\)
      • +
      +
      +
      Shaping of S, T, KS and GS + +
      +

      general_conf_shaping_S_T_KS_GS.png +

      +

      Figure 35: Generalized Plant to shape \(S\), \(T\), \(KS\) and \(GS\)

      +
      + +
      +
      P = [ W1  -W1*G*W3 -G*W1
             0    0        W2
             1   -G*W3    -G];
       
      @@ -2278,7 +2313,7 @@ When shaping multiple closed-loop transfer functions, one should be verify caref -
      +

      Mathematical relations are linking the closed-loop transfer functions. For instance, the sensitivity function \(S(s)\) and the complementary sensitivity function \(T(s)\) as link by the following well known relation: @@ -2311,7 +2346,7 @@ The control bandwidth is clearly limited by physical constrains such as sampling Similar relationship can be found for \(T\), \(KS\) and \(GS\).

      -
      +

      Determine the approximate norms of \(T\), \(KS\) and \(GS\) for large loop gains (\(|G(j\omega) K(j\omega)| \gg 1\)) and small loop gains (\(|G(j\omega) K(j\omega)| \ll 1\)).

      @@ -2329,7 +2364,7 @@ You can follows this procedure for \(T\), \(KS\) and \(GS\):
      Answer

      -The obtained constrains are shown in Figure 34. +The obtained constrains are shown in Figure 36.

      @@ -2342,17 +2377,17 @@ However, in some frequency bands, the norms do not depend on the controller and

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

      -
      +

      h-infinity-4-blocs-constrains.png

      -

      Figure 34: Shaping the Gang of Four: Limitations

      +

      Figure 36: Shaping the Gang of Four: Limitations

      -
      +

      The order (resp. number of state) of the controller given by the \(\mathcal{H}_\infty\) synthesis is equal to the order (resp. number of state) of the weighted generalized plant. It is thus equal to the sum of the number of state of the non-weighted generalized plant and the number of state of all the weighting functions. @@ -2371,43 +2406,418 @@ Two approaches can be used to obtain controllers with reasonable order:

      -
      -

      6 Mixed-Sensitivity \(\mathcal{H}_\infty\) Control - Example

      +
      +

      6 Mixed-Sensitivity \(\mathcal{H}_\infty\) Control - Example

      - +

      -
      -

      6.1 Problem

      +
      +

      6.1 Control Problem

      +
      +
        +
      • [ ] 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
      +
      -
      -

      6.2 Typical Procedure

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

      6.3 Step 1 - Shaping of the Sensitivity Function

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

      6.4 Step 2 - Shaping of

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

      7 Conclusion

      +
      +

      6.2 Control Design Procedure

      +
      +
      Table 7: Typical Specifications and corresponding wanted norms of open and closed loop tansfer functions
      + + +++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
         
      Response Time  
      Robustness  
      Limitation of the Command Amplitude  
      Filtering of the measurement noise  
      + + +
      +

      mixed_sensitivity_control_schematic.png +

      +

      Figure 37: Generalized Plant used for the Mixed Sensitivity Synthesis

      +
      + +
        +
      • \(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];
      +
      +
      +
      +
      + +
      +

      6.3 Step 1 - Shaping of \(S\)

      +
      +

      + +

      + +

      +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]');
      +
      +
      +
      +
      + +
      +

      6.4 Step 2 - Shaping of \(KS\)

      +
      +

      + +

      + +
      +
      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]');
      +
      +
      +
      +
      + +
      +

      6.5 Step 3 - Shaping of \(T\)

      +
      +

      + +

      + +
      +
      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]');
      +
      +
      +
      +
      +
      + +
      +

      7 Conclusion

      - +

      -
      -

      Resources

      -
      +
      +

      Resources

      +

      @@ -2420,7 +2830,7 @@ Two approaches can be used to obtain controllers with reasonable order:

      Author: Dehaeze Thomas

      -

      Created: 2020-12-01 mar. 18:49

      +

      Created: 2020-12-02 mer. 11:00

      diff --git a/index.org b/index.org index 9b37339..0867d17 100644 --- a/index.org +++ b/index.org @@ -41,15 +41,24 @@ * Introduction :ignore: -This document is structured as follows: -- As $\mathcal{H}_\infty$ Control is a /model based/ control technique, a short introduction to model based control is given in Section [[sec:model_based_control]] +The purpose of this document is to give a /practical introduction/ to the wonderful world of $\mathcal{H}_\infty$ Control. + +No attend is made to provide an exhaustive treatment of the subject. +$\mathcal{H}_\infty$ Control is a very broad topic and entire books are written on it. +Therefore, for more advanced discussion, please have a look at the recommended references at the bottom of this document. + +When possible, Matlab scripts used for the example/exercises are provided such that you can easily test them on your computer. + + +The general structure of this document is as follows: +- A short introduction to /model based control/ is given in Section [[sec:model_based_control]] - Classical /open/ loop shaping method is presented in Section [[sec:open_loop_shaping]]. - It is also shown that $\mathcal{H}_\infty$ synthesis can be used for /open/ loop shaping. -- $\mathcal{H}_\infty$ - Important concepts such as the $\mathcal{H}_\infty$ norm and the generalized plant are introduced. -- A -- Finally, an complete example of the - is performed in Section [[sec:h_infinity_mixed_sensitivity]]. + It is also shown that $\mathcal{H}_\infty$ synthesis can be used for /open/ loop shaping +- Important concepts indispensable for $\mathcal{H}_\infty$ control such as the $\mathcal{H}_\infty$ norm and the generalized plant are introduced in Section [[sec:h_infinity_introduction]] +- A very important step in $\mathcal{H}_\infty$ control is to express the control specifications (performances, robustness, etc.) as an $\mathcal{H}_\infty$ optimization problem. Such procedure is described in Section [[sec:modern_interpretation_specification]] +- One of the most useful use of the $\mathcal{H}_\infty$ control is the shaping of closed-loop transfer functions. + Such technique is presented in Section [[sec:closed-loop-shaping]] +- Finally, complete examples of the use of $\mathcal{H}_\infty$ Control for practical problems are provided in Section [[sec:h_infinity_mixed_sensitivity]]. * Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) @@ -309,20 +318,6 @@ You can follow this generic procedure: \end{align} #+HTML: #+end_exercice -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
      blocks. - -For now, I am manually adding #+HTML:
      Code and #+HTML:
      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 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 [[fig:mech_sys_1dof_inertial_contr]] into a classical feedback form as shown in Figure [[fig:open_loop_shaping]]. @@ -947,7 +942,7 @@ where: The procedure to convert a typical control architecture as the one shown in Figure [[fig:classical_feedback_tracking]] 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 +2. Remove $K$ and rearrange the inputs and outputs to match the generalized configuration shown in Figure [[fig:general_plant]] #+begin_src latex :file classical_feedback_tracking.pdf \begin{tikzpicture} @@ -994,7 +989,8 @@ The procedure to convert a typical control architecture as the one shown in Figu #+end_src #+begin_exercice - Compute the Generalized plant of corresponding to the tracking control architecture shown in Figure [[fig:classical_feedback_tracking]] + 1. Convert the tracking control architecture shown in Figure [[fig:classical_feedback_tracking]] to a generalized configuration + 2. Compute the transfer function matrix using Matlab as a function or $K$ and $G$ #+name: fig:classical_feedback_tracking #+caption: Classical Feedback Control Architecture (Tracking) @@ -1003,11 +999,12 @@ The procedure to convert a typical control architecture as the one shown in Figu #+HTML:
      Hint First, define the signals of the generalized plant: - Exogenous inputs: $w = r$ - - Signals to be minimized: $z_1 = \epsilon$, $z_2 = u$ - - Control signals: $v = y$ - - Control inputs: $u$ + - Signals to be minimized: + Usually, we want to minimize the tracking errors $\epsilon$ and the control signal $u$: $z = [\epsilon,\ u]$ + - Controller inputs: this is the signal at the input of the controller: $v = \epsilon$ + - Control inputs: signal generated by the controller: $u$ - Then, Remove $K$ and rearrange the inputs and outputs. + Then, Remove $K$ and rearrange the inputs and outputs as in Figure [[fig:general_plant]]. #+HTML:
      #+HTML:
      Answer @@ -1016,18 +1013,18 @@ The procedure to convert a typical control architecture as the one shown in Figu #+name: fig:mixed_sensitivity_ref_tracking #+caption: Generalized plant of the Classical Feedback Control Architecture (Tracking) [[file:figs/mixed_sensitivity_ref_tracking.png]] + + Using Matlab, the generalized plant can be defined as follows: + #+begin_src matlab :tangle no :eval no + P = [1 -G; + 0 1; + 1 -G] + P.InputName = {'w', 'u'}; + P.OutputName = {'e', 'u', 'v'}; + #+end_src #+HTML:
      #+end_exercice -Using Matlab, the generalized plant can be defined as follows: -#+begin_src matlab :tangle no :eval no - P = [1 -G; - 0 1; - 1 -G] - P.InputName = {'w', 'u'}; - P.OutputName = {'e', 'u', 'v'}; -#+end_src - * Modern Interpretation of Control Specifications <> @@ -2087,6 +2084,58 @@ Depending on the closed-loop transfer function being shaped, different general c - $W_2$ is used to shape $T$ #+HTML:
      +*** S GS :ignore: +#+HTML:
      Shaping of S and GS +#+begin_src latex :file general_conf_shaping_S_GS.pdf + \begin{tikzpicture} + % Blocs + \node[block] (G) {$G$}; + + \node[addb={+}{-}{}{}{}, left=0.8 of G] (addw) {}; + \node[block, above right=0.4 and 0.8 of G] (W2) {$W_2$}; + \node[block, above=0.5 of W2] (W1) {$W_1$}; + + \begin{scope}[on background layer] + \node[fit={(addw.west |- G.south) (W1.north east)}, inner sep=10pt, draw, dashed, fill=black!20!white] (P) {}; + \node[above] at (P.north) {Weighted Generalized Plant $P$}; + \end{scope} + + \node[block, below=0.6 of P] (K) {$K$}; + + \coordinate[right=1.0 of W1] (z); + \coordinate[above left=2.4 and 1.0 of addw] (w); + + % Connections + \draw[->] (addw.east) -- (G.west); + \draw[->] ($(addw.east)+(0.3, 0)$)node[branch]{} |- (W1.west); + \draw[->] ($(G.east)+(0.3, 0)$)node[branch]{} |- (W2.west); + + \draw[->] (W1.east) -- (W1-|z) node[above left](z1){$z_1$}; + \draw[->] (W2.east) -- (W2-|z) node[above left](z2){$z_2$}; + + \draw[->] (G.east) -- (G-|z) |- node[near start, right]{$v$} (K.east); + \draw[->] (K.west) -| node[near end, left]{$u$} (addw-|w) -- (addw.west); + + \draw[->] (w) node[above right]{$w$} -| (addw.north); + \end{tikzpicture} +#+end_src + +#+name: fig:general_conf_shaping_S_GS +#+caption: Generalized Plant to shape $S$ and $GS$ +#+RESULTS: +[[file:figs/general_conf_shaping_S_GS.png]] + +#+name: lst:general_plant_S_GS +#+caption: General Plant definition corresponding to Figure [[fig:general_conf_shaping_S_GS]] +#+begin_src matlab :eval no :tangle no + P = [W1 -W1 + G*W2 -G*W2 + G -G]; +#+end_src +- $W_1$ is used to shape $S$ +- $W_2$ is used to shape $GS$ +#+HTML:
      + *** S T KS :ignore: #+HTML:
      Shaping of S, T and KS #+begin_src latex :file general_conf_shaping_S_T_KS.pdf @@ -2147,6 +2196,66 @@ Depending on the closed-loop transfer function being shaped, different general c - $W_3$ is used to shape $T$ #+HTML:
      +*** S T GS :ignore: +#+HTML:
      Shaping of S, T and GS +#+begin_src latex :file general_conf_shaping_S_T_GS.pdf + \begin{tikzpicture} + % Blocs + \node[block] (G) {$G$}; + + \node[addb={+}{-}{}{}{}, left=0.8 of G] (addw) {}; + \node[block, above right=0.4 and 0.8 of G] (W3) {$W_3$}; + \node[block, above=0.2 of W3] (W2) {$W_2$}; + \node[block, above=0.2 of W2] (W1) {$W_1$}; + + \coordinate (addwin) at ($(addw.west)+(-0.5, 0)$); + + \begin{scope}[on background layer] + \node[fit={(addwin|-G.south) (W1.north east)}, inner sep=8pt, draw, dashed, fill=black!20!white] (P) {}; + \node[above] at (P.north) {Weighted Generalized Plant $P$}; + \end{scope} + + \node[block, below=0.6 of P] (K) {$K$}; + + \coordinate[right=0.8 of W1] (z); + \coordinate[above left=2.4 and 0.8 of addwin] (w); + + % Connections + \draw[->] (addw.east) -- (G.west); + \draw[->] ($(addw.east)+( 0.3, 0)$)node[branch]{} |- (W1.west); + \draw[->] ($(G.east)+(0.3, 0)$)node[branch]{} |- (W2.west); + \draw[->] (addwin)node[branch]{} |- (W3.west); + + \draw[->] (W1.east) -- (W1-|z) node[above left](z1){$z_1$}; + \draw[->] (W2.east) -- (W2-|z) node[above left](z2){$z_2$}; + \draw[->] (W3.east) -- (W3-|z) node[above left](z3){$z_3$}; + + \draw[->] (G.east) -- (G-|z) |- node[near start, right]{$v$} (K.east); + \draw[->] (K.west) -| node[near end, left]{$u$} (addw-|w) -- (addw.west); + + \draw[->] (w) node[above right]{$w$} -| (addw.north); + \end{tikzpicture} +#+end_src + +#+name: fig:general_conf_shaping_S_T_GS +#+caption: Generalized Plant to shape $S$, $T$ and $GS$ +#+RESULTS: +[[file:figs/general_conf_shaping_S_T_GS.png]] + +#+name: lst:general_plant_S_T_GS +#+caption: General Plant definition corresponding to Figure [[fig:general_conf_shaping_S_T_GS]] +#+begin_src matlab :eval no :tangle no + P = [W1 -W1 + G*W2 -G*W2 + 0 W3 + G -G]; +#+end_src + +- $W_1$ is used to shape $S$ +- $W_2$ is used to shape $GS$ +- $W_3$ is used to shape $T$ +#+HTML:
      + *** S T KS GS :ignore: #+HTML:
      Shaping of S, T, KS and GS #+begin_src latex :file general_conf_shaping_S_T_KS_GS.pdf @@ -2341,13 +2450,1409 @@ These regions are summarized in Figure [[fig:h-infinity-4-blocs-constrains]]. * Mixed-Sensitivity $\mathcal{H}_\infty$ Control - Example <> -** Problem +** Control Problem -** Typical Procedure +- [ ] Control Diagram -** Step 1 - Shaping of the Sensitivity Function -** Step 2 - Shaping of +- 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 + + +#+begin_src matlab + 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 +#+end_src + +#+begin_src matlab + % 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; +#+end_src + +** Control Design Procedure + +| | | | +|-------------------------------------+---+---| +| Response Time | | | +| Robustness | | | +| Limitation of the Command Amplitude | | | +| Filtering of the measurement noise | | | + +#+begin_src latex :file mixed_sensitivity_control_schematic.pdf + \begin{tikzpicture} + % Blocs + \node[block] (G) {$G$}; + + \node[addb={+}{-}{}{}{}, right=0.8 of G] (addw) {}; + \node[block, above right=0.4 and 0.8 of addw] (W3) {$W_3$}; + \node[block, above=0.2 of W3] (W2) {$W_2$}; + \node[block, above=0.2 of W2] (W1) {$W_1$}; + + \coordinate (Gin) at ($(G.west)+(-0.5, 0)$); + + \begin{scope}[on background layer] + \node[fit={(Gin|-G.south) (W1.north east)}, inner sep=8pt, draw, dashed, fill=black!20!white] (P) {}; + \node[above] at (P.north) {Weighted Generalized Plant $P$}; + \end{scope} + + \node[block, below=0.6 of P] (K) {$K$}; + + \coordinate[right=0.8 of W1] (z); + \coordinate[above left=1.4 and 1.3 of G] (w); + + % Connections + \draw[->] (G.east) -- (addw.west); + \draw[->] ($(addw.east)+(0.3, 0)$)node[branch]{} |- (W1.west); + \draw[->] (Gin)node[branch]{} |- (W2.west); + \draw[->] ($(G.east)+(0.3, 0)$)node[branch]{} |- (W3.west); + + \draw[->] (W1.east) -- (W1-|z) node[above left](z1){$z_1$}; + \draw[->] (W2.east) -- (W2-|z) node[above left](z2){$z_2$}; + \draw[->] (W3.east) -- (W3-|z) node[above left](z3){$z_3$}; + + \draw[->] (addw.east) -- (addw-|z) |- node[near start, right]{$v$} (K.east); + \draw[->] (K.west) -| node[near end, left]{$u$} (G-|w) -- (G.west); + + \draw[->] (w) node[above right]{$w$} -| (addw.north); + \end{tikzpicture} +#+end_src + +#+name: fig:mixed_sensitivity_control_schematic +#+caption: Generalized Plant used for the Mixed Sensitivity Synthesis +#+RESULTS: +[[file:figs/mixed_sensitivity_control_schematic.png]] + +- $W_1$ is used to shape $S$ +- $W_2$ is used to shape $KS$ +- $W_3$ is used to shape $T$ + +#+begin_src matlab + P = [1 -G + 0 1 + 0 G + 1 -G]; +#+end_src + +** Step 1 - Shaping of $S$ +<> + +We start with the shaping of $S$ alone. +To not constrain $KS$ and $T$, we set small values for $W_2$ and $W_3$ + +#+begin_src matlab + W2 = tf(1e-8); + W3 = tf(0.1); +#+end_src + +#+begin_src matlab + W1 = generateWeight('G0', 1e3, ... + 'G1', 1/2, ... + 'Gc', sqrt(2), 'wc', 2*pi*10, ... + 'n', 1); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K1 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +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' +#+end_example + +#+begin_src matlab + Z1 = lft(P, K1); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 4, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('KS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + % 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); +#+end_src + +#+begin_src matlab + z = lsim(Z1sim, [r;n;d], t); + y1 = z(:,1); + u1 = z(:,2); +#+end_src + +#+begin_src matlab + 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]'); +#+end_src + +** Step 2 - Shaping of $KS$ +<> + +#+begin_src matlab + W2 = generateWeight('G0', 5e-7, ... + 'G1', 1e-3, ... + 'Gc', 1e-6, 'wc', 2*pi*100, ... + 'n', 1); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K2 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +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' +#+end_example + +#+begin_src matlab + Z2 = lft(P, K2); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('KS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + Z2sim = lft(Psim, K2); +#+end_src + +#+begin_src matlab + z = lsim(Z2sim, [r;n;d], t); + y2 = z(:,1); + u2 = z(:,2); +#+end_src + +#+begin_src matlab + 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]'); +#+end_src + +** Step 3 - Shaping of $T$ +<> + +#+begin_src matlab + W3 = generateWeight('G0', 1e-1, ... + 'G1', 1e4, ... + 'Gc', 1, 'wc', 2*pi*30, ... + 'n', 3); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K3 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +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' +#+end_example + +#+begin_src matlab + Z3 = lft(P, K3); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('KS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + Z3sim = lft(Psim, K3); +#+end_src + +#+begin_src matlab + z = lsim(Z3sim, [r;n;d], t); + y3 = z(:,1); + u3 = z(:,2); +#+end_src + +#+begin_src matlab + 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]'); +#+end_src + +* Alternative Example - Shaping of $GS$ :noexport: +** Plant +#+begin_src matlab + 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 +#+end_src + +#+begin_src matlab + % Generate Input Signals + t = 0:1e-4:0.5; + + 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-5; + + d = zeros(size(t)); + d(t>0.25) = 5e-4; +#+end_src + +#+begin_src matlab + P = [1 -1 + G -G + 0 1 + G -G]; +#+end_src + +** Step 1 - Shaping of $S$ +#+begin_src matlab + W1 = generateWeight('G0', 1e3, ... + 'G1', 1/2, ... + 'Gc', 1, 'wc', 2*pi*10, ... + 'n', 1); +#+end_src + +#+begin_src matlab + W2 = tf(1e-8); + W3 = tf(1e-8); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K1 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +K1 = hinfsyn(Pw, 1, 1, 'Display', 'on'); + + Test bounds: 18.5 <= gamma <= 18.9 + + gamma X>=0 Y>=0 rho(XY)<1 p/f + 1.87e+01 0.0e+00 0.0e+00 5.225e+01 p + + Best performance (actual): 0.252 +'org_babel_eoe' +ans = + 'org_babel_eoe' +#+end_example + +#+begin_src matlab + Z1 = lft(P, K1); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('GS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + % 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); +#+end_src + +#+begin_src matlab + z = lsim(Z1sim, [r;n;d], t); + y1 = z(:,1); + u1 = z(:,2); +#+end_src + +#+begin_src matlab + 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]'); +#+end_src + +** Step 2 - Shaping of $GS$ +<> + +#+begin_src matlab + W2 = generateWeight('G0', 3e5, ... + 'G1', 1e9, ... + 'Gc', 1e6, 'wc', 2*pi*200, ... + 'n', 1); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K2 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +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' +#+end_example + +#+begin_src matlab + Z2 = lft(P, K2); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('GS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + Z2sim = lft(Psim, K2); +#+end_src + +#+begin_src matlab + z = lsim(Z2sim, [r;n;d], t); + y2 = z(:,1); + u2 = z(:,2); +#+end_src + +#+begin_src matlab + 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]'); +#+end_src + +** Step 3 - Shaping of $T$ +<> + +#+begin_src matlab + W3 = generateWeight('G0', 1e-1, ... + 'G1', 1e4, ... + 'Gc', 1, 'wc', 2*pi*80, ... + 'n', 2); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K3 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +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' +#+end_example + +#+begin_src matlab + Z3 = lft(P, K3); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('GS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + Z3sim = lft(Psim, K3); +#+end_src + +#+begin_src matlab + z = lsim(Z3sim, [r;n;d], t); + y3 = z(:,1); + u3 = z(:,2); +#+end_src + +#+begin_src matlab + 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]'); +#+end_src + +** Evolution of Controller + +#+begin_src matlab :exports none + freqs = logspace(-1, 4, 1000); + + figure; + tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile([2,1]); + hold on; + plot(freqs, abs(squeeze(freqresp(K1, freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(K2, freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(K3, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylabel('Magnitude'); set(gca, 'XTickLabel',[]); + hold off; + + ax2 = nexttile; + hold on; + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(K1, freqs, 'Hz'))))); + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(K2, freqs, 'Hz'))))); + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(K3, freqs, 'Hz'))))); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + yticks(-360:90:360); ylim([-270, 90]); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + linkaxes([ax1,ax2],'x'); + xlim([freqs(1), freqs(end)]); +#+end_src + +Loop Gain +#+begin_src matlab :exports none + freqs = logspace(-1, 4, 1000); + + figure; + tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile([2,1]); + hold on; + plot(freqs, abs(squeeze(freqresp(G*K1, freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(G*K2, freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(G*K3, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylabel('Magnitude'); set(gca, 'XTickLabel',[]); + hold off; + + ax2 = nexttile; + hold on; + plot(freqs, 180/pi*angle(squeeze(freqresp(G*K1, freqs, 'Hz')))); + plot(freqs, 180/pi*angle(squeeze(freqresp(G*K2, freqs, 'Hz')))); + plot(freqs, 180/pi*angle(squeeze(freqresp(G*K3, freqs, 'Hz')))); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + yticks(-360:90:360); ylim([-270, 90]); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + linkaxes([ax1,ax2],'x'); + xlim([freqs(1), freqs(end)]); +#+end_src + +* Alternative Example - Follow Triangle Reference :noexport: +** Plant +#+begin_src matlab + 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 +#+end_src + +#+begin_src matlab + % Generate Input Signals + t = 0:1e-4:0.5; + + r = zeros(size(t)); + r(t>0.1 & t<=0.2) = 1e-2*(t(t>0.1 & t<=0.2)-0.1); + r(t>0.2) = 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.3) = 5e-4; +#+end_src + +#+begin_src matlab + P = [1 -1 + G -G + 0 1 + G -G]; +#+end_src + +** Step 1 - Shaping of $S$ +#+begin_src matlab + W1 = generateWeight('G0', 1e3, ... + 'G1', 1/2, ... + 'Gc', 1, 'wc', 2*pi*10, ... + 'n', 1); +#+end_src + +#+begin_src matlab + W2 = tf(1e-8); + W3 = tf(1e-8); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K1 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +K1 = hinfsyn(Pw, 1, 1, 'Display', 'on'); + + Test bounds: 18.5 <= gamma <= 18.9 + + gamma X>=0 Y>=0 rho(XY)<1 p/f + 1.87e+01 0.0e+00 0.0e+00 5.225e+01 p + + Best performance (actual): 0.252 +'org_babel_eoe' +ans = + 'org_babel_eoe' +#+end_example + +#+begin_src matlab + Z1 = lft(P, K1); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('GS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + % 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); +#+end_src + +#+begin_src matlab + z = lsim(Z1sim, [r;n;d], t); + y1 = z(:,1); + u1 = z(:,2); +#+end_src + +#+begin_src matlab + figure; + tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); + + nexttile; + hold on; + plot(t, y1); + plot(t, r, 'k--'); + 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]'); +#+end_src + +** Step 1 - Shaping of $S$ with 2 integrators +#+begin_src matlab + W1 = generateWeight('G0', 1e3, ... + 'G1', 1/2, ... + 'Gc', 1, 'wc', 2*pi*15, ... + 'n', 2); +#+end_src + +#+begin_src matlab + W2 = tf(1e-8); + W3 = tf(1e-8); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K1b = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +K1 = hinfsyn(Pw, 1, 1, 'Display', 'on'); + + Test bounds: 18.5 <= gamma <= 18.9 + + gamma X>=0 Y>=0 rho(XY)<1 p/f + 1.87e+01 0.0e+00 0.0e+00 5.225e+01 p + + Best performance (actual): 0.252 +'org_babel_eoe' +ans = + 'org_babel_eoe' +#+end_example + +#+begin_src matlab + Z1b = lft(P, K1b); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z1b(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z1b(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('GS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z1b(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), 'k--'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + % 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]; + Z1bsim = lft(Psim, K1b); +#+end_src + +#+begin_src matlab + z = lsim(Z1bsim, [r;n;d], t); + y1b = z(:,1); + u1b = z(:,2); +#+end_src + +#+begin_src matlab + figure; + tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); + + nexttile; + hold on; + plot(t, y1); + plot(t, y1b); + plot(t, r, 'k--'); + hold off; + xlabel('Time [s]'); ylabel('Output $y$ [m]'); + + nexttile; + hold on; + plot(t, u1); + plot(t, u1b); + hold off; + xlabel('Time [s]'); ylabel('Control Input $u$ [N]'); +#+end_src + +** Step 2 - Shaping of $GS$ +<> + +#+begin_src matlab + W2 = tf(4e5); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K2 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +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' +#+end_example + +#+begin_src matlab + Z2 = lft(P, K2); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('GS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + Z2sim = lft(Psim, K2); +#+end_src + +#+begin_src matlab + z = lsim(Z2sim, [r;n;d], t); + y2 = z(:,1); + u2 = z(:,2); +#+end_src + +#+begin_src matlab + figure; + tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); + + nexttile; + hold on; + plot(t, y1); + plot(t, y2); + plot(t, r, 'k--'); + 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]'); +#+end_src + +** Step 3 - Shaping of $T$ +<> + +#+begin_src matlab + W3 = generateWeight('G0', 1e-1, ... + 'G1', 1e4, ... + 'Gc', 1, 'wc', 2*pi*70, ... + 'n', 2); +#+end_src + +#+begin_src matlab + Pw = blkdiag(W1, W2, W3, 1)*P; +#+end_src + +#+begin_src matlab :results output replace + K3 = hinfsyn(Pw, 1, 1, 'Display', 'on'); +#+end_src + +#+RESULTS: +#+begin_example +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' +#+end_example + +#+begin_src matlab + Z3 = lft(P, K3); +#+end_src + +#+begin_src matlab :exports none + freqs = logspace(-1, 3, 1000); + + figure; + tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(1,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(1,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('S'); + + ax2 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(2,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(2,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('GS'); + + ax3 = nexttile; + hold on; + plot(freqs, abs(squeeze(freqresp(Z1(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z2(3,1), freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(Z3(3,1), freqs, 'Hz')))); + plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), 'k-'); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + xlabel('Frequency [Hz]'); ylabel('Magnitude'); + title('T'); +#+end_src + +#+begin_src matlab + Z3sim = lft(Psim, K3); +#+end_src + +#+begin_src matlab + z = lsim(Z3sim, [r;n;d], t); + y3 = z(:,1); + u3 = z(:,2); +#+end_src + +#+begin_src matlab + figure; + tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None'); + + nexttile; + hold on; + plot(t, y1); + plot(t, y2); + plot(t, y3); + plot(t, r, 'k--'); + 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]'); +#+end_src + +** Evolution of Controller + +#+begin_src matlab :exports none + freqs = logspace(-1, 4, 1000); + + figure; + tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile([2,1]); + hold on; + plot(freqs, abs(squeeze(freqresp(K1, freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(K2, freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(K3, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylabel('Magnitude'); set(gca, 'XTickLabel',[]); + hold off; + + ax2 = nexttile; + hold on; + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(K1, freqs, 'Hz'))))); + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(K2, freqs, 'Hz'))))); + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(K3, freqs, 'Hz'))))); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + yticks(-360:90:360); ylim([-270, 90]); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + linkaxes([ax1,ax2],'x'); + xlim([freqs(1), freqs(end)]); +#+end_src + +Loop Gain +#+begin_src matlab :exports none + freqs = logspace(-1, 4, 1000); + + figure; + tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + + ax1 = nexttile([2,1]); + hold on; + plot(freqs, abs(squeeze(freqresp(G*K1, freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(G*K2, freqs, 'Hz')))); + plot(freqs, abs(squeeze(freqresp(G*K3, freqs, 'Hz')))); + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylabel('Magnitude'); set(gca, 'XTickLabel',[]); + hold off; + + ax2 = nexttile; + hold on; + plot(freqs, 180/pi*angle(squeeze(freqresp(G*K1, freqs, 'Hz')))); + plot(freqs, 180/pi*angle(squeeze(freqresp(G*K2, freqs, 'Hz')))); + plot(freqs, 180/pi*angle(squeeze(freqresp(G*K3, freqs, 'Hz')))); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + yticks(-360:90:360); ylim([-270, 90]); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + linkaxes([ax1,ax2],'x'); + xlim([freqs(1), freqs(end)]); +#+end_src * Conclusion <>