lecture-h-infinity/index.org

#+TITLE: A brief and practical introduction to $\mathcal{H}_\infty$ Control
:DRAWER:
#+STARTUP: overview

#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas

#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP:   ../index.html

#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="https://research.tdehaeze.xyz/css/style.css"/>
#+HTML_HEAD: <script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.js"></script>

#+HTML_MATHJAX: align: center tagside: right font: TeX

#+CSL_STYLE: ieee.csl

#+PROPERTY: header-args:matlab  :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :output-dir figs
#+PROPERTY: header-args:matlab+ :tangle no
#+PROPERTY: header-args:matlab+ :mkdirp yes

#+PROPERTY: header-args:shell  :eval no-export

#+PROPERTY: header-args:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/tikz/org/}{config.tex}")
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
#+PROPERTY: header-args:latex+ :results file raw replace
#+PROPERTY: header-args:latex+ :buffer no
#+PROPERTY: header-args:latex+ :tangle no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports results
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
#+PROPERTY: header-args:latex+ :post pdf2svg(file=*this*, ext="png")
:END:

* Introduction                                                        :ignore:

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
- 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)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

#+begin_src matlab :tangle no
  addpath('matlab')
#+end_src

* Introduction to Model Based Control
<<sec:model_based_control>>

** Introduction                                                      :ignore:

# - Section [[sec:model_based_control_methodology]]
# - Section [[sec:comp_classical_modern_robust_control]]
# - Section [[sec:example_system]]

** Model Based Control - Methodology
<<sec:model_based_control_methodology>>

The typical methodology for *Model Based Control* techniques is schematically shown in Figure [[fig:control-procedure]].

It consists of three steps:
1. *Identification or modeling*: a mathematical model $G(s)$ representing the plant dynamics is obtained
2. *Translate the specifications into mathematical criteria*:
   - _Specifications_: Response Time, Noise Rejection, Maximum input amplitude, Robustness, ...
   - _Mathematical Criteria_: Cost Function, Shape of transfer function, Phase/Gain margin, Roll-Off, ...
3. *Synthesis*: research of a controller $K(s)$ that satisfies the specifications for the model of the system

#+begin_src latex :file control-procedure.pdf
  \begin{tikzpicture}
    \node[addb={+}{}{}{}{-}] (addsub) at (0, 0){};

    \node[block, right=1.5 of addsub] (controller) {Controller};
    \node[block, right=1.5 of controller] (plant) {Plant};

    \node[block, above=1 of controller] (controller_design) {Synthesis};
    \node[block, above=1 of plant] (model_plant) {Model};

    \draw[<-] (addsub.west) -- ++(-1, 0) node[above right]{$r$};

    \draw[->] (addsub) -- (controller.west) node[above left]{$\epsilon$};
    \draw[->] (controller) -- (plant.west) node[above left]{$u$};
    \draw[->] (plant.east) -- ++(1, 0) node[above left]{$y$};
    \draw[] ($(plant.east) + (0.5, 0)$) -- ++(0, -1);
    \draw[->] ($(plant.east) + (0.5, -1)$) -| (addsub.south);

    \draw[->, dashed] (plant) -- node[midway, right, labelc, solid]{1} (model_plant);
    \draw[->, dashed] (controller_design) --node[midway, right, labelc, solid]{3} (controller);
    \draw[->, dashed] (model_plant) -- (controller_design);
    \draw[<-, dashed] (controller_design.west) -- node[midway, above, labelc, solid]{2} ++(-1, 0) node[left, style={align=center}]{Specifications};
  \end{tikzpicture}
#+end_src

#+name: fig:control-procedure
#+caption: Typical Methodoly for Model Based Control
#+RESULTS:
[[file:figs/control-procedure.png]]

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


In Section [[sec:open_loop_shaping]], steps 2 and 3 will be described for a control techniques called *classical (open-)loop shaping*.

Then, steps 2 and 3 for the *$\mathcal{H}_\infty$ Loop Shaping* of closed-loop transfer functions will be discussed in Sections [[sec:modern_interpretation_specification]], [[sec:closed-loop-shaping]] and [[sec:h_infinity_mixed_sensitivity]].

** From Classical Control to Robust Control
<<sec:comp_classical_modern_robust_control>>

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


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


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


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


The three presented control methods are compared in Table [[tab:comparison_control_methods]].

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

#+name: tab:comparison_control_methods
#+caption: Table summurazing the main differences between classical, modern and robust control
| <l>                     |                <c>                 |                 <c>                  |                                   <c>                                   |
|                         |        *Classical Control*         |           *Modern Control*           |                            *Robust Control*                             |
|-------------------------+------------------------------------+--------------------------------------+-------------------------------------------------------------------------|
| *Date*                  |               1930-                |                1960-                 |                                  1980-                                  |
|-------------------------+------------------------------------+--------------------------------------+-------------------------------------------------------------------------|
| *Tools*                 |         Transfer Functions         |       State Space formulation        | Systems and Signals Norms ($\mathcal{H}_\infty$, $\mathcal{H}_2$ Norms) |
|                         |           Nyquist Plots            |          Riccati Equations           |                     Closed Loop Transfer Functions                      |
|                         |             Bode Plots             |                                      |                        Open/Closed Loop Shaping                         |
|                         |       Phase and Gain margins       |                                      |                           Weighting Functions                           |
|                         |                                    |                                      |                               Disk margin                               |
|                         |                                    |                                      |                      Singular Value Decomposition                       |
|-------------------------+------------------------------------+--------------------------------------+-------------------------------------------------------------------------|
| *Control Architectures* | Proportional, Integral, Derivative |       Full State Feedback, LQR       |                      General Control Configuration                      |
|                         |            Leads, Lags             |         Kalman Filters, LQG          |                            Generalized Plant                            |
|-------------------------+------------------------------------+--------------------------------------+-------------------------------------------------------------------------|
| *Advantages*            |          Study Stability           |         Automatic Synthesis          |                           Automatic Synthesis                           |
|                         |               Simple               |                 MIMO                 |                                  MIMO                                   |
|                         |              Natural               |         Optimization Problem         |                          Optimization Problem                           |
|                         |                                    |                                      |                          Guaranteed Robustness                          |
|                         |                                    |                                      |                   Easy specification of performances                    |
|-------------------------+------------------------------------+--------------------------------------+-------------------------------------------------------------------------|
| *Disadvantages*         |           Manual Method            |       No Guaranteed Robustness       |                  Required knowledge of specific tools                   |
|                         |             Only SISO              | Difficult Rejection of Perturbations |               Need a reasonably good model of the system                |
|                         | No clear way to limit input usage  |                                      |                                                                         |

#+begin_src latex :file robustness_performance.pdf
  \begin{tikzpicture}
    % Scale
    \def\yscale{0.8}
    \def\xscale{1.0}

    % Colors
    \def\colorstart{blue}
    \def\colorend{red}

    % Axis
    \draw [->] (-0.5,0) -- (10*\xscale,0) node[below left]{Robustness};
    \draw [->] (0,-0.5) -- (0,10*\yscale) node[below left, rotate=90, anchor=south east]{Performance};

    % Color Bar
    \shade[draw, bottom color=\colorstart, top color=\colorend, fill opacity=0.5] (10*\xscale, 1*\yscale) rectangle (11*\xscale, 9*\yscale);
    \node[rotate=90, above] at (10*\xscale, 5*\yscale) {Required information on plant};
    \node[above] at (10.5*\xscale, 1*\yscale) {little};
    \node[below] at (10.5*\xscale, 9*\yscale) {large};

    % ===================================
    % Classical Control
    % ===================================
    % Control Types
    \node[align=center] (pid)  at (7.0*\xscale, 1.2*\yscale) {PID\\Lead-Lag};

    \begin{scope}[on background layer]
      % Control Families
      \node[ellipse, draw, dashed, minimum width=3.0*\xscale cm, minimum height=2.0*\yscale cm,
          fill=\colorstart!90!\colorend, fill opacity=0.5, text opacity=1]
          (classicalcontrol) at (pid) {};
    \end{scope}
    \node[above, align=center] at (classicalcontrol.north) {\textbf{Classical control} (1930)\\{\small SISO, Manual Method}};
    % ===================================


    % ===================================
    % Modern Control
    % ===================================
    % Control Types
    \node[align=center] (lqg)  at (2.0*\xscale, 7.5*\yscale) {LQR\\LQG};

    \begin{scope}[on background layer]
      \node[ellipse, draw, dashed, minimum width=2.0*\xscale cm, minimum height=2.0*\yscale cm,
          fill=\colorstart!20!\colorend, fill opacity=0.5, text opacity=1]
          (moderncontrol) at (lqg) {};
    \end{scope}
    \node[above, align=center] at (moderncontrol.north) {\textbf{Modern control} (1960)\\{\small MIMO, Optimal}};
    % ===================================


    % ===================================
    % Robust Control
    % ===================================
    % Control Types
    \node[align=center] (hinf) at (4.5*\xscale, 4.8*\yscale) {$H_\infty$\\$H_2$};
    \node[] (mu) at (5.5*\xscale, 4.8*\yscale) {$\mu$};


    \begin{scope}[on background layer]
      \node[ellipse, draw, dashed, minimum width=3.0*\xscale cm, minimum height=2.5*\yscale cm,
          shade, left color=\colorstart!50!\colorend, right color=\colorstart!10!\colorend, fill opacity=0.5, text opacity=1]
          (robustcontrol) at ($0.5*(hinf)+0.5*(mu)$) {};
    \end{scope}
    \node[above, align=center] at (robustcontrol.north) {\textbf{Robust control} (1990)\\{\small MIMO, Robust}};
    % ===================================
  \end{tikzpicture}
#+end_src

# #+name: fig:robustness_performance
# #+caption: 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.
# #+RESULTS:
# [[file:figs/robustness_performance.png]]

** Example System
<<sec:example_system>>

Throughout this document, multiple examples and practical application of presented control strategies will be provided.
Most of them will be applied on a physical system presented in this section.

This system is shown in Figure [[fig:mech_sys_1dof_inertial_contr]].
It could represent an active suspension stage supporting a payload.
The /inertial/ motion of the payload is measured using an inertial sensor and this is feedback to a force actuator.
Such system could be used to actively isolate the payload (disturbance rejection problem) or to make it follow a trajectory (tracking problem).

The notations used on Figure [[fig:mech_sys_1dof_inertial_contr]] are listed and described in Table [[tab:example_notations]].

#+begin_src latex :file mech_sys_1dof_inertial_contr.pdf
  \begin{tikzpicture}
    % Parameters
    \def\massw{3}
    \def\massh{1}
    \def\spaceh{1.8}

    % Ground
    \draw[] (-0.5*\massw, 0) -- (0.5*\massw, 0);
    % Mass
    \draw[fill=white] (-0.5*\massw, \spaceh) rectangle (0.5*\massw, \spaceh+\massh) node[pos=0.5](m){$m$};

    % Spring, Damper, and Actuator
    \draw[spring]   (-0.3*\massw, 0) -- (-0.3*\massw, \spaceh) node[midway, left=0.1]{$k$};
    \draw[damper]   ( 0, 0) -- ( 0, \spaceh) node[midway, left=0.3]{$c$};
    \draw[actuator] ( 0.3*\massw, 0) -- (0.3*\massw, \spaceh) node[midway](F){};

    % Displacements
    \draw[dashed] (0.5*\massw, 0) -- ++(0.5, 0);
    \draw[->] (0.6*\massw, 0) -- ++(0, 0.5) node[below right]{$d$};

    % Inertial Sensor
    \node[inertialsensor] (inertials) at (0.5*\massw, \spaceh+\massh){};
    \node[addb={+}{-}{}{}{}, right=0.8 of inertials] (subf) {};

    \node[block, below=0.4 of subf] (K){$K(s)$};

    \draw[->] (inertials.east) node[above right]{$y$} -- (subf.west);
    \draw[->] (subf.south) -- (K.north) node[above right]{$\epsilon$};
    \draw[<-] (subf.north) -- ++(0, 0.6) node[below right]{$r$};
    \draw[->] (K.south) |- (F.east) node[above right]{$u$};
  \end{tikzpicture}
#+end_src

#+name: fig:mech_sys_1dof_inertial_contr
#+caption: 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.
#+RESULTS:
[[file:figs/mech_sys_1dof_inertial_contr.png]]

#+name: tab:example_notations
#+caption: Example system variables
| *Notation*         | *Description*                                                  | *Value*        | *Unit*    |
|--------------------+----------------------------------------------------------------+----------------+-----------|
| $m$                | Payload's mass to position / isolate                           | $10$           | [kg]      |
| $k$                | Stiffness of the suspension system                             | $10^6$         | [N/m]     |
| $c$                | Damping coefficient of the suspension system                   | $400$          | [N/(m/s)] |
| $y$                | Payload absolute displacement (measured by an inertial sensor) |                | [m]       |
| $d$                | Ground displacement, it acts as a disturbance                  |                | [m]       |
| $u$                | Actuator force                                                 |                | [N]       |
| $r$                | Wanted position of the mass (the reference)                    |                | [m]       |
| $\epsilon = r - y$ | Position error                                                 |                | [m]       |
| $K$                | Feedback controller                                            | to be designed | [N/m]     |

#+begin_exercice
Derive the following open-loop transfer functions:
  \begin{align}
    G(s)   &= \frac{y}{u} \\
    G_d(s) &= \frac{y}{d}
  \end{align}

#+HTML: <details><summary>Hint</summary>
You can follow this generic procedure:
1. List all applied forces ot the mass: Actuator force, Stiffness force (Hooke's law), ...
2. Apply the Newton's Second Law on the payload
   \[ m \ddot{y} = \Sigma F \]
3. Transform the differential equations into the Laplace domain:
   \[ \frac{d\ \cdot}{dt} \Leftrightarrow \cdot \times s \]
4. Write $y(s)$ as a function of $u(s)$ and $w(s)$
#+HTML: </details>

#+HTML: <details><summary>Results</summary>
\begin{align}
  G(s)   &= \frac{1}{m s^2 + cs + k} \\
  G_d(s) &= \frac{cs + k}{m s^2 + cs + k}
\end{align}
#+HTML: </details>
#+end_exercice

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 architecture as shown in Figure [[fig:open_loop_shaping]].

#+begin_src latex :file classical_feedback_test_system.pdf
  \begin{tikzpicture}
    \node[addb={+}{}{}{}{-}] (addfb) at (0, 0){};
    \node[block, right=0.8 of addfb] (K){$K(s)$};
    \node[block, right=0.8 of K] (G){$G(s)$};
    \node[addb={+}{}{}{}{}, right=0.8 of G] (addd){};
    \node[block, above=0.5 of addd] (Gd){$G_d(s)$};

    \draw[<-] (addfb.west) -- ++(-0.8, 0) node[above right]{$r$};
    \draw[->] (addfb.east) -- (K.west) node[above left]{$\epsilon$};
    \draw[->] (K.east) -- (G.west) node[above left]{$u$};
    \draw[->] (G.east) -- (addd.west);
    \draw[<-] (Gd.north) -- ++(0, 0.8) node[below right]{$d$};
    \draw[->] (Gd.south) -- (addd.north);
    \draw[->] (addd.east) -- ++(1.2, 0);
    \draw[->] ($(addd.east) + (0.6, 0)$) node[branch]{} node[above]{$y$} -- ++(0, -1.0) -| (addfb.south);
  \end{tikzpicture}
#+end_src

#+name: fig:classical_feedback_test_system
#+caption: Block diagram corresponding to the example system of Figure [[fig:mech_sys_1dof_inertial_contr]]
#+RESULTS:
[[file:figs/classical_feedback_test_system.png]]


Let's define the system parameters on Matlab.
#+begin_src matlab +n
  k = 1e6; % Stiffness [N/m]
  c = 4e2; % Damping [N/(m/s)]
  m = 10; % Mass [kg]
#+end_src

And now the system dynamics $G(s)$ and $G_d(s)$.
#+begin_src matlab +n -r
  G = 1/(m*s^2 + c*s + k); % Plant
  Gd = (c*s + k)/(m*s^2 + c*s + k); % Disturbance
#+end_src

The Bode plots of $G(s)$ and $G_d(s)$ are shown in Figures [[fig:bode_plot_example_afm]] and [[fig:bode_plot_example_Gd]].

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

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

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G, 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(G, 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

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/bode_plot_example_afm.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:bode_plot_example_afm
#+caption: Bode plot of the plant $G(s)$
#+RESULTS:
[[file:figs/bode_plot_example_afm.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  tiledlayout(1, 1, 'TileSpacing', 'None', 'Padding', 'None');
  nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(Gd, freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Magnitude');
  hold off;
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/bode_plot_example_Gd.pdf', 'width', 'wide', 'height', 'small');
#+end_src

#+name: fig:bode_plot_example_Gd
#+caption: Magnitude of the disturbance transfer function $G_d(s)$
#+RESULTS:
[[file:figs/bode_plot_example_Gd.png]]

* Classical Open Loop Shaping
<<sec:open_loop_shaping>>

** Introduction                                                      :ignore:

After an introduction to classical Loop Shaping in Section [[sec:open_loop_shaping_introduction]], a practical example is given in Section [[sec:loop_shaping_example]].
Such Loop Shaping is usually performed manually with tools coming from the classical control theory.

However, the $\mathcal{H}_\infty$ synthesis can be used to automate the Loop Shaping process.
This is presented in Section [[sec:h_infinity_open_loop_shaping]] and applied on the same example in Section [[sec:h_infinity_open_loop_shaping_example]].

** Introduction to Loop Shaping
<<sec:open_loop_shaping_introduction>>

#+begin_definition
*Loop Shaping* refers to a control design procedure that involves explicitly shaping the magnitude of the *Loop Transfer Function* $L(s)$.
#+end_definition

#+begin_definition
The *Loop Gain* (or Loop transfer function) $L(s)$ usually refers to as the product of the controller and the plant (see Figure [[fig:open_loop_shaping]]):
\begin{equation}
  L(s) = G(s) \cdot K(s) \label{eq:loop_gain}
\end{equation}

Its name comes from the fact that this is actually the "gain around the loop".

#+begin_src latex :file open_loop_shaping.pdf
  \begin{tikzpicture}
    \node[addb={+}{}{}{}{-}] (addsub) at (0, 0){};

    \node[block, right=0.8 of addsub] (K) {$K(s)$};
    \node[below] at (K.south) {Controller};
    \node[block, right=0.8 of K] (G) {$G(s)$};
    \node[below] at (G.south) {Plant};

    \draw[<-] (addsub.west) -- ++(-0.8, 0) node[above right]{$r$};

    \draw[->] (addsub) -- (K.west) node[above left]{$\epsilon$};
    \draw[->] (K.east) -- (G.west) node[above left]{$u$};
    \draw[->] (G.east) -- ++(0.8, 0) node[above left]{$y$};
    \draw[] ($(G.east) + (0.5, 0)$) -- ++(0, -1.4);
    \draw[->] ($(G.east) + (0.5, -1.4)$) -| (addsub.south);

    \draw [decoration={brace, raise=5pt}, decorate] (K.north west) -- node[above=6pt]{$L(s)$} (G.north east);
  \end{tikzpicture}
#+end_src

#+name: fig:open_loop_shaping
#+caption: Classical Feedback Architecture
[[file:figs/open_loop_shaping.png]]
#+end_definition

This synthesis method is one of main way controllers are design in the classical control theory.
It is widely used and generally successful as many characteristics of the closed-loop system depend on the shape of the open loop gain $L(s)$ such as:
- *Good Tracking*: $L$ large
- *Good disturbance rejection*: $L$ large
- *Attenuation of measurement noise on plant output*: $L$ small
- *Small magnitude of input signal*: $L$ small
- *Nominal stability*: $L$ small (RHP zeros and time delays)
- *Robust stability*: $L$ small (neglected dynamics)

The shaping of the Loop Gain is done manually by combining several leads, lags, notches...
This process is very much simplified by the fact that the loop gain $L(s)$ depends *linearly* on $K(s)$ eqref:eq:loop_gain.
A typical wanted Loop Shape is shown in Figure [[fig:open_loop_shaping_shape]].
Another interesting Loop shape called "Bode Step" is described in cite:lurie02_system_archit_trades_using_bode.

#+begin_src latex :file open_loop_shaping_shape.pdf
  \begin{tikzpicture}
    % Phase Axis
    \draw[->] (-0.3, -0.5) -- ++(8, 0) node[above]{$\omega$}; \draw[<-] (0, 0)
    node[left]{$\angle L(j\omega)$} -- ++(0, -2.3);

    % Gain Axis
    \draw[->] (-0.3, 2) -- ++(8, 0) node[above]{$\omega$}; \draw[->] (0, 0.5) --
    ++(0, 3) node[left]{$\left|L(j\omega)\right|$};

    % Gain Slopes
    \draw[shift={(0,2)}] (0.5, 1.25) -- node[midway, above]{$-2$} (2, 0.5) --
    node[midway, above]{$-1$} (6, -0.5) -- node[midway, below left]{$-2$} (7.5,
    -1.25);

    % Forbiden region
    \path[shift={(0,1.8)}, fill=red!50!white] (0.5, 1.25) -- (2, 0.5) -| coordinate[near start](lfshaping) cycle;
    \path[shift={(0,2.2)}, fill=red!50!white] (6, -0.5) -- (7.5, -1.25) |- coordinate[near end](hfshaping) cycle;

    \draw[<-] (lfshaping) -- ++(0, -0.8) node[below, align=center]{{\scriptsize Ref. tracking}\\{\scriptsize Dist. rejection}};
    \draw[<-] (hfshaping) -- ++(0,  0.8)  node[above, align=center]{{\scriptsize Noise attenuation}};

    % Crossover frequency
    \node[below] (wc) at (4,2){$\omega_c$};
    \draw[<-] (wc.south) -- ++(0, -0.4) node[below, align=center]{{\scriptsize Bandwidth}};

    % Phase
    \draw[] (0.5, -2) -- (2, -2)[out=0, in=-180] to (4, -1.25)[out=0, in=-180] to
    (6, -2) -- (7.5, -2); \draw[] (0.5, -2) -- (2, -2)[out=0, in=-180] to (4,
    -1.25)[out=0, in=-180] to (6, -2) -- (7.5, -2);

    % Phase Margin
    \draw[->, dashed] (4, -2) -- (4, -1.25) node[above]{{\scriptsize Phase Margin}};
    \draw[dashed] (0, -2) node[left]{$-\pi$} -- (7.5, -2);
  \end{tikzpicture}
#+end_src

#+name: fig:open_loop_shaping_shape
#+caption: Typical Wanted Shape for the Loop Gain $L(s)$
#+RESULTS:
[[file:figs/open_loop_shaping_shape.png]]

The shaping of *closed-loop* transfer functions is obviously not as simple as they don't depend linearly on $K(s)$.
But this is were the $\mathcal{H}_\infty$ Synthesis will be useful!
More details on that in Sections [[sec:modern_interpretation_specification]] and [[sec:closed-loop-shaping]].

** Example of Manual Open Loop Shaping
<<sec:loop_shaping_example>>

#+begin_exampl
Let's take our example system described in Section [[sec:example_system]] and design a controller using the Open-Loop shaping synthesis approach.
The specifications are:
1. *Disturbance rejection*: Highest possible rejection below 1Hz
2. *Positioning speed*: Bandwidth of approximately 10Hz
3. *Noise attenuation*: Roll-off of -40dB/decade past 30Hz
4. *Robustness*: Gain margin > 3dB and Phase margin > 30 deg
#+end_exampl

#+begin_exercice
Using =SISOTOOL=, design a controller that fulfills the specifications.

#+begin_src matlab :eval no :tangle no
  sisotool(G)
#+end_src

#+HTML: <details><summary>Hint</summary>
You can follow this procedure:
1. In order to have good disturbance rejection at low frequency, add a simple or double *integrator*
2. In terms of the loop gain, the *bandwidth* can be defined at the frequency $\omega_c$ where $|l(j\omega_c)|$ first crosses 1 from above.
   Therefore, adjust the *gain* such that $L(j\omega)$ crosses 1 at around 10Hz
3. The roll-off at high frequency for noise attenuation should already be good enough.
   If not, add a *low pass filter*
4. Add a *Lead* centered around the crossover frequency (10 Hz) and tune it such that sufficient phase margin is added.
   Verify that the gain margin is good enough.
#+HTML: </details>
#+end_exercice

Let's say we came up with the following controller.
#+begin_src matlab
  K = 14e8 * ...              % Gain
      1/(s^2) * ...           % Double Integrator
      1/(1 + s/2/pi/40) * ... % Low Pass Filter
      (1 + s/(2*pi*10/sqrt(8)))/(1 + s/(2*pi*10*sqrt(8))); % Lead
#+end_src

The bode plot of the Loop Gain is shown in Figure [[fig:loop_gain_manual_afm]] and we can verify that we have the wanted stability margins using the =margin= command:
#+begin_src matlab
  [Gm, Pm, ~, Wc] = margin(G*K)
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
  data2orgtable([Gm; Pm; Wc/2/pi], {'Gain Margin $> 3$ [dB]', 'Phase Margin $> 30$ [deg]', 'Crossover $\approx 10$ [Hz]'}, {'Requirements', 'Manual Method'}, ' %.1f ');
#+end_src

#+RESULTS:
| Requirements                | Manual Method |
|-----------------------------+---------------|
| Gain Margin $> 3$ [dB]      |           3.1 |
| Phase Margin $> 30$ [deg]   |          35.4 |
| Crossover $\approx 10$ [Hz] |          10.1 |

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

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

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G*K, freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  hold off;
  ylim([1e-4, 1e1])

  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G*K, freqs, 'Hz')))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360); ylim([-360, 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/loop_gain_manual_afm.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:loop_gain_manual_afm
#+caption: Bode Plot of the obtained Loop Gain $L(s) = G(s) K(s)$
#+RESULTS:
[[file:figs/loop_gain_manual_afm.png]]

** $\mathcal{H}_\infty$ Loop Shaping Synthesis
<<sec:h_infinity_open_loop_shaping>>

The synthesis of controllers based on the Loop Shaping method can be automated using the $\mathcal{H}_\infty$ Synthesis.

Using Matlab, it can be easily performed using the =loopsyn= command:
#+begin_src matlab :eval no :tangle no
  K = loopsyn(G, Lw);
#+end_src
where:
- =G= is the (LTI) plant
- =Lw= is the wanted loop shape
- =K= is the synthesize controller

#+begin_seealso
  Matlab documentation of =loopsyn= ([[https://www.mathworks.com/help/robust/ref/loopsyn.html][link]]).
#+end_seealso

Therefore, by just providing the wanted loop shape and the plant model, the $\mathcal{H}_\infty$ Loop Shaping synthesis generates a /stabilizing/ controller such that the obtained loop gain $L(s)$ matches the specified one with an accuracy $\gamma$.

Even though we will not go into details and explain how such synthesis is working, an example is provided in the next section.

** Example of the $\mathcal{H}_\infty$ Loop Shaping Synthesis
<<sec:h_infinity_open_loop_shaping_example>>

To apply the $\mathcal{H}_\infty$ Loop Shaping Synthesis, the wanted shape of the loop gain should be determined from the specifications.
This is summarized in Table [[tab:open_loop_shaping_specifications]].

Such shape corresponds to the typical wanted Loop gain Shape shown in Figure [[fig:open_loop_shaping_shape]].

#+name: tab:open_loop_shaping_specifications
#+caption: Wanted Loop Shape corresponding to each specification
|                         | Specification                               | Corresponding Loop Shape                                        |
|-------------------------+---------------------------------------------+-----------------------------------------------------------------|
| *Disturbance Rejection* | Highest possible rejection below 1Hz        | Slope of -40dB/decade at low frequency to have a high loop gain |
| *Positioning Speed*     | Bandwidth of approximately 10Hz             | $L$ crosses 1 at 10Hz: $\vert L_w(j2 \pi 10)\vert = 1$          |
| *Noise Attenuation*     | Roll-off of -40dB/decade past 30Hz          | Roll-off of -40dB/decade past 30Hz                              |
| *Robustness*            | Gain margin > 3dB and Phase margin > 30 deg | Slope of -20dB/decade near the crossover                        |

Then, a (stable, minimum phase) transfer function $L_w(s)$ should be created that has the same gain as the wanted shape of the Loop gain.
For this example, a double integrator and a lead centered on 10Hz are used.
Then the gain is adjusted such that the $|L_w(j2 \pi 10)| = 1$.

Using Matlab, we have:
#+begin_src matlab
  Lw = 2.3e3 * ...
       1/(s^2) * ... % Double Integrator
       (1 + s/(2*pi*10/sqrt(3)))/(1 + s/(2*pi*10*sqrt(3))); % Lead
#+end_src

The $\mathcal{H}_\infty$ open loop shaping synthesis is then performed using the =loopsyn= command:
#+begin_src matlab
  [K, ~, GAM] = loopsyn(G, Lw);
#+end_src

The obtained Loop Gain is shown in Figure [[fig:open_loop_shaping_hinf_L]] and matches the specified one by a factor $\gamma \approx 2$.

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

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

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G*K, freqs, 'Hz'))), 'DisplayName', '$L(s)$');
  plot(freqs, abs(squeeze(freqresp(Lw, freqs, 'Hz'))), 'k--', 'DisplayName', '$L_w(s)$');
  plot(freqs, abs(squeeze(freqresp(Lw, freqs, 'Hz')))*GAM, 'k-.', 'DisplayName', '$L_w(s) / \gamma$, $L_w(s) \cdot \gamma$');
  plot(freqs, abs(squeeze(freqresp(Lw, freqs, 'Hz')))/GAM, 'k-.', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  hold off;
  legend('location', 'northeast');
  ylim([1e-4, 1e2]);

  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G*K, freqs, 'Hz')))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360); ylim([-360, 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/open_loop_shaping_hinf_L.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:open_loop_shaping_hinf_L
#+caption: Obtained Open Loop Gain $L(s) = G(s) K(s)$ and comparison with the wanted Loop gain $L_w$
#+RESULTS:
[[file:figs/open_loop_shaping_hinf_L.png]]


#+begin_important
When using the $\mathcal{H}_\infty$ Synthesis, it is usually recommended to analyze the obtained controller.

This is usually done by breaking down the controller into simple elements such as low pass filters, high pass filters, notches, leads, etc.
#+end_important

Let's briefly analyze the obtained controller which bode plot is shown in Figure [[fig:open_loop_shaping_hinf_K]]:
- two integrators are used at low frequency to have the wanted low frequency high gain
- a lead is added centered with the crossover frequency to increase the phase margin
- a notch is added at the resonance of the plant to increase the gain margin (this is very typical of $\mathcal{H}_\infty$ controllers, and can be an issue, more info on that latter)

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

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

  ax1 = nexttile([2,1]);
  plot(freqs, abs(squeeze(freqresp(K, freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  hold off;

  ax2 = nexttile;
  plot(freqs, 180/pi*angle(squeeze(freqresp(K, freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360); ylim([-180, 90]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/open_loop_shaping_hinf_K.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:open_loop_shaping_hinf_K
#+caption: Obtained controller $K$ using the open-loop $\mathcal{H}_\infty$ shaping
#+RESULTS:
[[file:figs/open_loop_shaping_hinf_K.png]]


Let's finally compare the obtained stability margins of the $\mathcal{H}_\infty$ controller and of the manually developed controller in Table [[tab:open_loop_shaping_compare]].

#+begin_src matlab :exports none
  [Gm_2, Pm_2, ~, Wc_2] = margin(G*K)
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
  data2orgtable([Gm, Gm_2; Pm, Pm_2; Wc/2/pi, Wc_2/2/pi], {'Gain Margin $> 3$ [dB]', 'Phase Margin $> 30$ [deg]', 'Crossover $\approx 10$ [Hz]'}, {'Specifications', 'Manual Method', '$\mathcal{H}_\infty$ Method'}, ' %.1f ');
#+end_src

#+name: tab:open_loop_shaping_compare
#+caption: Comparison of the characteristics obtained with the two methods
#+RESULTS:
| Specifications              | Manual Method | $\mathcal{H}_\infty$ Method |
|-----------------------------+---------------+-----------------------------|
| Gain Margin $> 3$ [dB]      |           3.1 |                        31.7 |
| Phase Margin $> 30$ [deg]   |          35.4 |                        54.7 |
| Crossover $\approx 10$ [Hz] |          10.1 |                         9.9 |

* A first Step into the $\mathcal{H}_\infty$ world
<<sec:h_infinity_introduction>>

** Introduction                                                      :ignore:

- Section [[sec:h_infinity_norm]]
- Section [[sec:h_infinity_synthesis]]
- Section [[sec:generalized_plant]]
- Section [[sec:h_infinity_general_synthesis]]
- Section [[sec:generalized_plant_derivation]]

** The $\mathcal{H}_\infty$ Norm
<<sec:h_infinity_norm>>

#+begin_definition
 The $\mathcal{H}_\infty$ norm is defined as the peak of the maximum singular value of the frequency response
 \begin{equation}
   \|G(s)\|_\infty = \max_\omega \bar{\sigma}\big( G(j\omega) \big)
 \end{equation}

 For a SISO system $G(s)$, it is simply the peak value of $|G(j\omega)|$ as a function of frequency:
 \begin{equation}
   \|G(s)\|_\infty = \max_{\omega} |G(j\omega)| \label{eq:hinf_norm_siso}
 \end{equation}
#+end_definition

#+begin_exampl
Let's compute the $\mathcal{H}_\infty$ norm of our test plant $G(s)$ using the =hinfnorm= function:
#+begin_src matlab :results value replace
  hinfnorm(G)
#+end_src

#+RESULTS:
: 7.9216e-06

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 [[fig:hinfinity_norm_siso_bode]].

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);
  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))), 'k-');
  plot([20, 100], [hinfnorm(G) hinfnorm(G)], 'k--');
  text(100, hinfnorm(G), '$\quad \|G\|_\infty$')
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Magnitude $|G(j\omega)|$');
  ylim([1e-8, 2e-5]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/hinfinity_norm_siso_bode.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:hinfinity_norm_siso_bode
#+caption: Example of the $\mathcal{H}_\infty$ norm of a SISO system
#+RESULTS:
[[file:figs/hinfinity_norm_siso_bode.png]]
#+end_exampl

** $\mathcal{H}_\infty$ Synthesis
<<sec:h_infinity_synthesis>>

#+begin_definition
  $\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.
#+end_definition

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


Then applying the $\mathcal{H}_\infty$ synthesis to a plant, the engineer work usually consists of the following steps
1. Write the problem as standard $\mathcal{H}_\infty$ problem
2. Translate the specifications as $\mathcal{H}_\infty$ norms of transfer functions
3. Make the synthesis and analyze the obtain controller
4. Reduce the order of the controller for implementation


Note that there are many ways to use the $\mathcal{H}_\infty$ Synthesis:
- Traditional $\mathcal{H}_\infty$ Synthesis (=hinfsyn= [[https://www.mathworks.com/help/robust/ref/hinfsyn.html][doc]])
- Open Loop Shaping $\mathcal{H}_\infty$ Synthesis (=loopsyn= [[https://www.mathworks.com/help/robust/ref/loopsyn.html][doc]])
- Mixed Sensitivity Loop Shaping (=mixsyn= [[https://www.mathworks.com/help/robust/ref/lti.mixsyn.html][doc]])
- Fixed-Structure $\mathcal{H}_\infty$ Synthesis (=hinfstruct= [[https://www.mathworks.com/help/robust/ref/lti.hinfstruct.html][doc]])
- Signal Based $\mathcal{H}_\infty$ Synthesis, and many more...

** The Generalized Plant
<<sec:generalized_plant>>

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

The generalized plant, usually noted $P(s)$, is shown in Figure [[fig:general_plant]].
It has two inputs and two outputs (both could contains many signals).
The meaning of the inputs and outputs are summarized in Table [[tab:notation_general]].

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

\begin{equation}
  \begin{bmatrix} z \\ v \end{bmatrix} = P \begin{bmatrix} w \\ u \end{bmatrix} = \begin{bmatrix} P_{11} & P_{12} \\ P_{21} & P_{22} \end{bmatrix} \begin{bmatrix} w \\ u \end{bmatrix}
\end{equation}

#+begin_src latex :file general_plant.pdf
  \begin{tikzpicture}
    \node[block={2.0cm}{2.0cm}] (P) {$P$};
    \node[above] at (P.north) {Generalized Plant};

    % Input and outputs coordinates
    \coordinate[] (inputw)  at ($(P.south west)!0.75!(P.north west)$);
    \coordinate[] (inputu)  at ($(P.south west)!0.25!(P.north west)$);
    \coordinate[] (outputz) at ($(P.south east)!0.75!(P.north east)$);
    \coordinate[] (outputv) at ($(P.south east)!0.25!(P.north east)$);

    % Connections and labels
    \draw[<-] (inputw) -- ++(-0.8, 0) node[above right]{$w$};
    \draw[<-] (inputu) -- ++(-0.8, 0) node[above right]{$u$};

    \draw[->] (outputz) -- ++(0.8, 0) node[above left]{$z$};
    \draw[->] (outputv) -- ++(0.8, 0) node[above left]{$v$};
  \end{tikzpicture}
#+end_src

#+name: fig:general_plant
#+caption: Inputs and Outputs of the generalized Plant
#+RESULTS:
[[file:figs/general_plant.png]]

#+begin_important
#+name: tab:notation_general
#+caption: Notations for the general configuration
| Notation | Meaning                                         |
|----------+-------------------------------------------------|
| $P$      | Generalized plant model                         |
| $w$      | Exogenous inputs: commands, disturbances, noise |
| $z$      | Exogenous outputs: signals to be minimized      |
| $v$      | Controller inputs: measurements                 |
| $u$      | Control signals                                 |
#+end_important

** The $\mathcal{H}_\infty$ Synthesis applied on the Generalized plant
<<sec:h_infinity_general_synthesis>>

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

#+begin_important
  - $\mathcal{H}_\infty$ Synthesis applied on the generalized plant ::
  Find a stabilizing controller $K$ that, using the sensed output $v$, generates a control signal $u$ such that the $\mathcal{H}_\infty$ norm of the closed-loop transfer function from $w$ to $z$ is minimized.

  After $K$ is found, the system is /robustified/ by adjusting the response around the unity gain frequency to increase stability margins.
#+end_important

#+begin_src latex :file general_control_names.pdf
  \begin{tikzpicture}

    % Blocs
    \node[block={2.0cm}{2.0cm}] (P) {$P$};
    \node[block={1.5cm}{1.5cm}, below=0.7 of P] (K) {$K$};

    % Input and outputs coordinates
    \coordinate[] (inputw)  at ($(P.south west)!0.75!(P.north west)$);
    \coordinate[] (inputu)  at ($(P.south west)!0.25!(P.north west)$);
    \coordinate[] (outputz) at ($(P.south east)!0.75!(P.north east)$);
    \coordinate[] (outputv) at ($(P.south east)!0.25!(P.north east)$);

    % Connections and labels
    \draw[<-] (inputw) node[above left, align=right]{(weighted)\\exogenous inputs\\$w$} -- ++(-1.5, 0);
    \draw[<-] (inputu) -- ++(-0.8, 0) |- node[left, near start, align=right]{control signals\\$u$} (K.west);

    \draw[->] (outputz) node[above right, align=left]{(weighted)\\exogenous outputs\\$z$} -- ++(1.5, 0);
    \draw[->] (outputv) -- ++(0.8, 0) |- node[right, near start, align=left]{sensed output\\$v$} (K.east);
  \end{tikzpicture}
#+end_src

#+name: fig:general_control_names
#+caption: General Control Configuration
#+RESULTS:
[[file:figs/general_control_names.png]]

Note that the closed-loop transfer function from $w$ to $z$ is:
\begin{equation}
  \frac{z}{w} = P_{11} + P_{12} K \big( I - P_{22} K \big)^{-1} P_{21} \triangleq F_l(P, K)
\end{equation}

Using Matlab, the $\mathcal{H}_\infty$ Synthesis applied on a Generalized plant can be applied using the =hinfsyn= command ([[https://www.mathworks.com/help/robust/ref/hinfsyn.html][documentation]]):
#+begin_src matlab :eval no :tangoe no
  K = hinfsyn(P, nmeas, ncont);
#+end_src
where:
- =P= is the generalized plant transfer function matrix
- =nmeas= is the number of sensed output (size of $v$)
- =ncont= is the number of control signals (size of $u$)
- =K= obtained controller that minimized the $\mathcal{H}_\infty$ norm from $w$ to $z$

** From a Classical Feedback Architecture to a Generalized Plant
<<sec:generalized_plant_derivation>>

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 shown in Figure [[fig:general_plant]]

#+begin_src latex :file classical_feedback_tracking.pdf
  \begin{tikzpicture}
    \node[addb={+}{}{}{}{-}] (addfb) at (0, 0){};
    \node[block, right=0.8 of addfb] (K){$K(s)$};
    \node[block, right=0.8 of K] (G){$G(s)$};

    \draw[<-] (addfb.west) -- ++(-0.8, 0) node[above right]{$r$};
    \draw[->] (addfb.east) -- (K.west) node[above left]{$\epsilon$};
    \draw[->] (K.east) -- (G.west) node[above left]{$u$};
    \draw[->] (G.east) -- ++(1.2, 0);
    \draw[->] ($(G.east) + (0.6, 0)$) node[branch]{} node[above]{$y$} -- ++(0, -0.8) -| (addfb.south);
  \end{tikzpicture}
#+end_src

#+begin_src latex :file mixed_sensitivity_ref_tracking.pdf
  \begin{tikzpicture}
    \node[block] (G) {$G(s)$};
    \node[addb={+}{-}{}{}{}, right=0.6 of G] (addw) {};
    \coordinate[above right=0.6 and 1.4 of addw] (u);
    \coordinate[above=0.6 of u] (epsilon);

    \coordinate[] (w) at ($(epsilon-|G.west)+(-1.4, 0)$);

    \node[block, below left=0.8 and 0 of addw] (K) {$K(s)$};

    % Connections
    \draw[->] (G.east) -- (addw.west);
    \draw[->] ($(addw.east)+(0.4, 0)$)node[branch]{} |- (epsilon) node[above left](z1){$\epsilon$};
    \draw[->] ($(G.west)+(-0.4, 0)$)node[branch](start){} |- (u) node[above left](z2){$u$};

    \draw[->] (addw.east) -- (addw-|z1) |- node[near start, right]{$v$} (K.east);
    \draw[->] (K.west) -| node[near end, left]{$u$} ($(G-|w)+(0.4, 0)$) -- (G.west);

    \draw[->] (w) node[above]{$w = r$} -| (addw.north);

    \draw [decoration={brace, raise=5pt}, decorate] (z1.north east) -- node[right=6pt]{$z$} (z2.south east);

    \begin{scope}[on background layer]
      \node[fit={(G.south-|start.west) ($(z1.north west)+(-0.4, 0)$)}, inner sep=6pt, draw, dashed, fill=black!20!white] (P) {};
      \node[below] at (P.north) {Generalized Plant $P(s)$};
    \end{scope}
  \end{tikzpicture}
#+end_src

#+begin_exercice
  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)
  [[file:figs/classical_feedback_tracking.png]]

#+HTML: <details><summary>Hint</summary>
  First, define the signals of the generalized plant:
  - Exogenous inputs: $w = r$
  - 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 as in Figure [[fig:general_plant]].
#+HTML: </details>

#+HTML: <details><summary>Answer</summary>
  The obtained generalized plant shown in Figure [[fig:mixed_sensitivity_ref_tracking]].

  #+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: </details>
#+end_exercice

* Modern Interpretation of Control Specifications
<<sec:modern_interpretation_specification>>

** Introduction                                                      :ignore:

- Section [[sec:closed_loop_tf]]
- Section [[sec:sensitivity_transfer_functions]]
- Section [[sec:module_margin]]
- Section [[sec:other_requirements]]

As shown in Section [[sec:open_loop_shaping]], 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)$.

However, the loop gain $L(s)$ does *not* directly give the performances of the closed-loop system, which are determined by the *closed-loop* transfer functions.

If we consider the feedback system shown in Figure [[fig:gang_of_four_feedback]], we can link to the following specifications to closed-loop transfer functions.
This is summarized in Table [[tab:spec_closed_loop_tf]].

#+name: tab:spec_closed_loop_tf
#+caption: Typical Specification and associated closed-loop transfer function
| Specification                  | Closed-Loop Transfer Function                 |
|--------------------------------+-----------------------------------------------|
| Reference Tracking             | From $r$ to $\epsilon$                        |
| Disturbance Rejection          | From $d$ to $y$                               |
| Measurement Noise Filtering    | From $n$ to $y$                               |
| Small Command Amplitude        | From $n,r,d$ to $u$                           |
| Stability                      | All closed-loop transfer function             |
| Robustness (stability margins) | Module margin (see Section [[sec:module_margin]]) |

#+begin_src latex :file gang_of_four_feedback.pdf
  \begin{tikzpicture}
    \node[addb={+}{}{}{}{-}] (addfb) at (0, 0){};
    \node[block, right=0.8 of addfb] (K){$K(s)$};
    \node[addb,  right=0.8 of K] (addd){};
    \node[block, right=0.8 of addd] (G){$G(s)$};
    \node[addb,  below right=0.4 and 0.2 of G] (addn){};

    \draw[<-] (addfb.west) -- ++(-0.8, 0) node[above right]{$r$};
    \draw[->] (addfb.east) -- (K.west) node[above left]{$\epsilon$};
    \draw[->] (K.east) -- (addd.west);
    \draw[<-] (addd.north) -- ++(0, 0.6) node[below right]{$d$};
    \draw[->] (addd.east) -- (G.west) node[above left]{$u$};
    \draw[->] (G.east) -- ++(1.6, 0) node[above left]{$y$};
    \draw[->] (G-|addn) node[branch]{} -- (addn.north);
    \draw[<-] (addn.east) -- ++(0.8, 0) node[above left]{$n$};
    \draw[->] (addn.west) -| (addfb.south);
  \end{tikzpicture}
#+end_src

#+name: fig:gang_of_four_feedback
#+caption: Simple Feedback Architecture
#+RESULTS:
[[file:figs/gang_of_four_feedback.png]]

** Closed Loop Transfer Functions
<<sec:closed_loop_tf>>

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

#+begin_exercice
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 [[fig:gang_of_four_feedback]].

#+HTML: <details><summary>Hint</summary>
Take one of the output (e.g. $y$), and write it as a function of the inputs $[d, r, n]$ going step by step around the loop:
\begin{aligned}
  y &= G u \\
    &= G (d + K \epsilon) \\
    &= G \big(d + K (r - n - y) \big) \\
    &= G d + GK r - GK n - GK y
\end{aligned}

Isolate $y$ at the right hand side, and finally obtain:
\[ y = \frac{GK}{1+ GK} r + \frac{G}{1 + GK} d - \frac{GK}{1 + GK} n \]

Do the same procedure for $u$ and $\epsilon$
#+HTML: </details>

#+HTML: <details><summary>Answer</summary>
The following equations should be obtained:
\begin{align}
  y        &= \frac{GK}{1 + GK} r + \frac{G}{1 + GK} d - \frac{GK}{1 + GK} n \\
  \epsilon &= \frac{1 }{1 + GK} r - \frac{G}{1 + GK} d - \frac{G }{1 + GK} n \\
  u        &= \frac{K }{1 + GK} r - \frac{1}{1 + GK} d - \frac{K }{1 + GK} n
\end{align}
#+HTML: </details>
#+end_exercice

#+begin_important
We can see that they are 4 different transfer functions describing the behavior of the system in Figure [[fig:gang_of_four_feedback]].
These called the *Gang of Four*:
\begin{align}
  S  &= \frac{1 }{1 + GK}, \quad \text{the sensitivity function} \\
  T  &= \frac{GK}{1 + GK}, \quad \text{the complementary sensitivity function} \\
  GS &= \frac{G }{1 + GK}, \quad \text{the load disturbance sensitivity function} \\
  KS &= \frac{K }{1 + GK}, \quad \text{the noise sensitivity function}
\end{align}
#+end_important

#+begin_seealso
If a feedforward controller is included, a *Gang of Six* transfer functions can be defined.
More on that in this [[https://www.youtube.com/watch?v=b_8v8scghh8][short video]].
#+end_seealso

And we have:
\begin{align}
  \epsilon &= S  r - GS d - GS n \\
  y        &= T  r + GS d - T  n \\
  u        &= KS r - S  d - KS n
\end{align}

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

** Sensitivity Function
<<sec:sensitivity_transfer_functions>>

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


Suppose we have developed a "/reference/" controller $K_r(s)$ and made three small changes to obtained three controllers $K_1(s)$, $K_2(s)$ and $K_3(s)$.
The obtained sensitivity functions are shown in Figure [[fig:sensitivity_shape_effect]] and the corresponding step responses are shown in Figure [[fig:sensitivity_shape_effect_step]].

The comparison of the sensitivity functions shapes and their effect on the step response is summarized in Table [[tab:compare_sensitivity_shapes]].

#+name: tab:compare_sensitivity_shapes
#+caption: Comparison of the sensitivity function shape and the corresponding step response for the three controller variations
| Controller | Sensitivity Function Shape                         | Change of the Step Response      |
|------------+----------------------------------------------------+----------------------------------|
| $K_1(s)$   | Larger bandwidth $\omega_b$                        | Faster rise time                 |
| $K_2(s)$   | Larger peak value $\Vert S\Vert_\infty$            | Large overshoot and oscillations |
| $K_3(s)$   | Larger low frequency gain $\vert S(j\cdot 0)\vert$ | Larger static error              |

#+begin_src matlab :exports none
  wc = 2*pi*1; L_w = 8; wi = 2*pi*0.02;
  Kr = 1/((s + wi)^2) * ... % Double Integrator
       (1 + s/(wc/sqrt(L_w)))/(1 + s/(wc*sqrt(L_w))); % Lead
  Kr = Kr/abs(evalfr(Kr*G, j*wc));

  wc = 2*pi*2; L_w = 8; wi = 2*pi*0.045;
  K1 = 1/((s + wi)^2) * ... % Double Integrator
       (1 + s/(wc/sqrt(L_w)))/(1 + s/(wc*sqrt(L_w))); % Lead
  K1 = K1/abs(evalfr(K1*G, j*wc));

  wc = 2*pi*1; L_w = 2; wi = 2*pi*0.03;
  K2 = 1/((s + wi)^2) * ... % Double Integrator
       (1 + s/(wc/sqrt(L_w)))/(1 + s/(wc*sqrt(L_w))); % Lead
  K2 = K2/abs(evalfr(K2*G, j*wc));

  wc = 2*pi*1; L_w = 8; wi = 2*pi*0.2;
  K3 = 1/((s + wi)^2) * ... % Double Integrator
       (1 + s/(wc/sqrt(L_w)))/(1 + s/(wc*sqrt(L_w))); % Lead
  K3 = K3/abs(evalfr(K3*G, j*wc));


  Sr = 1/(1 + Kr*G);
  S1 = 1/(1 + K1*G);
  S2 = 1/(1 + K2*G);
  S3 = 1/(1 + K3*G);

  Tr = Kr*G/(1 + Kr*G);
  T1 = K1*G/(1 + K1*G);
  T2 = K2*G/(1 + K2*G);
  T3 = K3*G/(1 + K3*G);
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(S1, freqs, 'Hz'))), 'DisplayName', '$K_1(s)$');
  plot(freqs, abs(squeeze(freqresp(S2, freqs, 'Hz'))), 'DisplayName', '$K_2(s)$');
  plot(freqs, abs(squeeze(freqresp(S3, freqs, 'Hz'))), 'DisplayName', '$K_3(s)$');
  plot(freqs, abs(squeeze(freqresp(Sr, freqs, 'Hz'))), 'k-', 'DisplayName', '$K_r(s)$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frquency [Hz]'); ylabel('Sensitivity Magnitude');
  hold off;
  legend('location', 'southeast');
  ylim([1e-3, 3]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/sensitivity_shape_effect.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:sensitivity_shape_effect
#+caption: Sensitivity function magnitude $|S(j\omega)|$ corresponding to the reference controller $K_r(s)$ and the three modified controllers $K_i(s)$
#+RESULTS:
[[file:figs/sensitivity_shape_effect.png]]

#+begin_src matlab :exports none
  t = linspace(0, 5, 1000);

  figure;
  hold on;
  plot(t, step(T1, t), 'DisplayName', '$K_1(s)$')
  plot(t, step(T2, t), 'DisplayName', '$K_2(s)$')
  plot(t, step(T3, t), 'DisplayName', '$K_3(s)$')
  plot(t, step(Tr, t), 'k-', 'DisplayName', '$K_r(s)$')
  hold off
  xlabel('Time [s]'); ylabel('Step Response');
  legend('location', 'northeast');
  ylim([0, 1.7])
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/sensitivity_shape_effect_step.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:sensitivity_shape_effect_step
#+caption: Step response (response from $r$ to $y$) for the different controllers
#+RESULTS:
[[file:figs/sensitivity_shape_effect_step.png]]

#+begin_src latex :file h-infinity-spec-S.pdf
  \begin{tikzpicture}
    \begin{axis}[%
      width=8cm,
      height=4cm,
      at={(0,0)},
      xmode=log,
      xmin=0.01,
      xmax=10000,
      ymin=-80,
      ymax=40,
      ylabel={Magnitude [dB]},
      xlabel={Frequency [Hz]},
      ytick={40, 20, 0, -20, -40, -60, -80},
      xminorgrids,
      yminorgrids
      ]
      \addplot [thick, color=black, forget plot]
      table[row sep=crcr]{%
        0.01 -60\\
        0.1  -60\\
        190   6\\
        10000 6\\
      };

      \draw[<-] (0.05, -60) -- (0.1, -70);
      \draw (0.1, -70) -- (2, -70) node[right, fill=white, draw]{\footnotesize{Small static error}};
      \draw[<-] (70, -3) -- (3, -3) node[left, fill=white, draw]{\footnotesize{Speed}};

      \draw[<-] (300, 6) -- (200, 20);
      \draw (200, 20) -- (10, 20) node[left, fill=white, draw]{\footnotesize{Robustness}};
    \end{axis}
  \end{tikzpicture}
#+end_src

#+begin_definition
- Closed-Loop Bandwidth ::
  The closed-loop bandwidth $\omega_b$ is the frequency where $|S(j\omega)|$ first crosses $1/\sqrt{2} = -3dB$ from below.

  In general, a large bandwidth corresponds to a faster rise time.
#+end_definition

#+begin_important
From the simple analysis above, we can draw a first estimation of the wanted shape for the sensitivity function in Figure [[fig:h-infinity-spec-S]].

The wanted characteristics on the magnitude of the sensitivity function are then:
- A small magnitude at low frequency to make the static errors small
- A wanted minimum closed-loop bandwidth in order to have fast rise time and good rejection of perturbations
- A small peak value in order to limit large overshoot and oscillations.
  This generally means higher robustness.
  This will become clear in the next section about the *module margin*.

#+name: fig:h-infinity-spec-S
#+caption: Typical wanted shape of the Sensitivity transfer function
[[file:figs/h-infinity-spec-S.png]]
#+end_important

** Robustness: Module Margin
<<sec:module_margin>>

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

#+begin_exampl
Let's consider the following plant $G_t(s)$:
#+begin_src matlab
  w0 = 2*pi*100;
  xi = 0.1;
  k = 1e7;

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

Let's say we have designed a controller $K_t(s)$ that gives the loop gain shown in Figure [[fig:phase_gain_margin_model_plant]].

#+begin_src matlab
  Kt = 1.2e6*(s + w0)/s;
#+end_src

The following characteristics can be determined from Figure [[fig:phase_gain_margin_model_plant]]:
- bandwidth of $\approx 10\, \text{Hz}$
- infinite gain margin (the phase of the loop-gain never reaches -180 degrees
- more than 90 degrees of phase margin

This might indicate very good robustness properties of the closed-loop system.

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

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

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(Gt*Kt, freqs, 'Hz'))), 'DisplayName', '$L(s)$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  hold off;
  legend('location', 'northeast');

  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gt*Kt, freqs, 'Hz')))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360); ylim([-200, 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/phase_gain_margin_model_plant.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:phase_gain_margin_model_plant
#+caption: Bode plot of the obtained Loop Gain $L(s)$
#+RESULTS:
[[file:figs/phase_gain_margin_model_plant.png]]


Now let's suppose the "real" plant $G_r(s)$ as a slightly lower damping factor:
#+begin_src matlab
  xi = 0.03;
#+end_src

#+begin_src matlab :exports none
  Gr = 1/k*(s/w0/4 + 1)/(s^2/w0^2 + 2*xi*s/w0 + 1);
#+end_src

The obtained "real" loop gain is shown in Figure [[fig:phase_gain_margin_real_plant]].
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.

It is confirmed by checking the stability of the closed loop system:
#+begin_src matlab :results value replace
  isstable(feedback(Gr,K))
#+end_src

#+RESULTS:
: 0

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

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

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(Gt*Kt, freqs, 'Hz'))), 'DisplayName', '$L(s)$');
  plot(freqs, abs(squeeze(freqresp(Gr*Kt, freqs, 'Hz'))), 'DisplayName', '$L_r(s)$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  hold off;
  legend('location', 'northeast');

  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gt*Kt, freqs, 'Hz')))));
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gr*Kt, freqs, 'Hz')))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360); ylim([-200, 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/phase_gain_margin_real_plant.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:phase_gain_margin_real_plant
#+caption: Bode plots of $L(s)$ (loop gain corresponding the nominal plant) and $L_r(s)$ (loop gain corresponding to the real plant)
#+RESULTS:
[[file:figs/phase_gain_margin_real_plant.png]]

Therefore, even a small change of the plant parameter makes the system unstable even though both the gain margin and the phase margin for the nominal plant are excellent.

This is due to the fact that the gain and phase margin are robustness indicators for a *pure* change or gain or a *pure* change of phase but not a combination of both.
#+end_exampl

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

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

- Nyquist Plot ::
  The Nyquist plot shows the evolution of $L(j\omega)$ in the complex plane from $\omega = 0 \to \infty$.
#+end_definition

#+begin_seealso
For more information about the /general/ Nyquist Stability Criteria, you may want to look at [[https://www.youtube.com/watch?v=sof3meN96MA][this]] video.
#+end_seealso

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

#+begin_definition
- 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.
#+end_definition

#+begin_exampl
A typical Nyquist plot is shown in Figure [[fig:module_margin_example]].
The gain, phase and module margins are graphically shown to have an idea of what they represent.

#+begin_src matlab :exports none
  % Example Plant
  k = 1e6; % Stiffness [N/m]
  c = 4e2; % Damping [N/(m/s)]
  m = 10; % Mass [kg]
  G = 1/(m*s^2 + c*s + k); % Plant

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

  L = G*K;
  L_resp = squeeze(freqresp(L, freqs, 'Hz'));

  % Module Margin
  Dm = min(abs(1 + L_resp));

  % Phase Gain Margin
  [Gm, Pm, Wcg, Wcp] = margin(L);

  freqs = logspace(0, 3, 1000);
  figure;
  hold on;
  % Gain Margin
  plot([-1, -1/Gm], [0, 0], '-', 'DisplayName', sprintf('$\\Delta G = %.1f$', Gm))

  % Phase Margin
  theta = -pi:0.01:-pi+Pm*pi/180;
  plot(cos(theta), sin(theta), '-', 'DisplayName', sprintf('$\\Delta \\phi = %.1f^o$', Pm));

  % Module Margin
  theta = 0 : 0.01 : 2*pi;
  plot(Dm*cos(theta)-1, Dm*sin(theta), '-', 'DisplayName', sprintf('$\\Delta M = %.1f$', Dm));

  % Nyquist Plot
  plot(real(L_resp),  imag(L_resp), 'k-', 'DisplayName', '$L(j\omega)$')

  plot(-1, 0, 'k*', 'HandleVisibility', 'off');

  hold off;
  xlabel('Real Axis'); ylabel('Imaginary Axis');
  xlim([-1.5, 0.5]); ylim([-1, 1]);
  axis equal;
  legend('location', 'southeast');
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/module_margin_example.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:module_margin_example
#+caption: Nyquist plot with visual indication of the Gain margin $\Delta G$, Phase margin $\Delta \phi$ and Module margin $\Delta M$
#+RESULTS:
[[file:figs/module_margin_example.png]]
#+end_exampl

As expected from Figure [[fig:module_margin_example]], 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}
  \Delta G \ge \frac{\Delta M}{\Delta M - 1}; \quad \Delta \phi \ge \frac{1}{\Delta M}
\end{equation}


Let's now try to express the Module margin $\Delta M$ as an $\mathcal{H}_\infty$ norm of a closed-loop transfer function:
\begin{align*}
\Delta M &= \text{minimum distance between } L(j\omega) \text{ and point } (-1) \\
         &= \min_\omega |L(j\omega) - (-1)| \\
         &= \min_\omega |1 + L(j\omega)| \\
         &= \frac{1}{\max_\omega \frac{1}{|1 + L(j\omega)|}} \\
         &= \frac{1}{\|S\|_\infty}
\end{align*}

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

\begin{equation}
  \Delta M = \frac{1}{\|S\|_\infty} \label{eq:module_margin_S}
\end{equation}

The wanted robustness of the closed-loop system can be specified by setting a maximum value on $\|S\|_\infty$.
#+end_important

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

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

#+begin_seealso
  To learn more about module/disk margin, you can check out [[https://www.youtube.com/watch?v=XazdN6eZF80][this]] video.
#+end_seealso

** TODO Other Requirements
<<sec:other_requirements>>

Interpretation of the $\mathcal{H}_\infty$ norm of systems:
- frequency by frequency attenuation / amplification

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

Consider an input sinusoidal signal $w(t) = \sin\left( \omega_0 t \right)$, then the output signal $z(t)$ will be equal to:
\[ z(t) = A \sin\left( \omega_0 t + \phi \right) \]
with:
- $A = |G_t(j\omega_0)|$ is the magnitude of $G_t(s)$ at $\omega_0$
- $\phi = \angle G_t(j\omega_0)$ is the phase of $G_t(s)$ at $\omega_0$


Noise Attenuation: typical wanted shape for $T$

#+name: tab:specification_modern
#+caption: Typical Specifications and corresponding wanted norms of open and closed loop tansfer functions
|                             | Open-Loop Shaping  | Closed-Loop Shaping                        |
|-----------------------------+--------------------+--------------------------------------------|
| Reference Tracking          | $L$ large          | $S$ small                                  |
| Disturbance Rejection       | $L$ large          | $GS$ small                                 |
| Measurement Noise Filtering | $L$ small          | $T$ small                                  |
| Small Command Amplitude     | $K$ and $L$ small  | $KS$ small                                 |
| Robustness                  | Phase/Gain margins | Module margin: $\Vert S\Vert_\infty$ small |

* $\mathcal{H}_\infty$ Shaping of closed-loop transfer functions
<<sec:closed-loop-shaping>>

** Introduction                                                      :ignore:

In the previous sections, we have seen that the performances of the system depends on the *shape* of the closed-loop transfer function.

Therefore, the synthesis problem is to design $K(s)$ such that closed-loop system is stable and such that various closed-loop transfer functions such as $S$, $KS$ and $T$ are shaped as wanted.
This is clearly not simple as these closed-loop transfer functions does not depend linearly on $K$.

But don't worry, the $\mathcal{H}_\infty$ synthesis will do this job for us!

This
Section [[sec:weighting_functions]]
Section [[sec:weighting_functions_design]]
Section [[sec:sensitivity_shaping_example]]
Section [[sec:shaping_multiple_tf]]

** How to Shape closed-loop transfer function? Using Weighting Functions!
<<sec:weighting_functions>>

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

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

#+begin_src latex :file loop_shaping_S_without_W.pdf
  \begin{tikzpicture}
    \node[block] (G) {$G(s)$};
    \node[addb={+}{-}{}{}{}, right=0.6 of G] (addw) {};
    \coordinate[above right=1.0 and 1.4 of addw] (epsilon);

    \coordinate[] (w) at ($(epsilon-|G.west)+(-1.0, 0)$);

    \node[block, below left=0.8 and 0 of addw] (K) {$K(s)$};

    % Connections
    \draw[->] (G.east) -- (addw.west);
    \draw[->] ($(addw.east)+(0.4, 0)$)node[branch]{} |- (epsilon) node[above](z1){$z = \epsilon$};

    \draw[->] (addw.east) -- (addw-|z1) |- node[near start, right]{$v$} (K.east);
    \draw[->] (K.west) -| node[near end, left]{$u$} ($(G-|w)+(0.1, 0)$) -- (G.west);

    \draw[->] (w) node[above]{$w = r$} -| (addw.north);

    \begin{scope}[on background layer]
      \node[fit={(G.south west) ($(z1.north west)+(-0.4, 0)$)}, inner sep=12pt, draw, dashed, fill=black!20!white] (P) {};
      \node[below] at (P.north) {Generalized Plant $P(s)$};
    \end{scope}
  \end{tikzpicture}
#+end_src

#+name: fig:loop_shaping_S_without_W
#+caption: Generalized Plant
#+RESULTS:
[[file:figs/loop_shaping_S_without_W.png]]


#+begin_important
The /trick/ is to include a *weighting function* $W_S(s)$ in the generalized plant as shown in Figure [[fig:loop_shaping_S_with_W]].

Now, the closed-loop transfer function from $w$ to $z$ is equal to $W_s(s)S(s)$ and applying the $\mathcal{H}_\infty$ synthesis to the /weighted/ generalized plant $\tilde{P}(s)$ will generate a controller $K(s)$ such that $\|W_s(s)S(s)\|_\infty$ is minimized.
#+end_important

Let's now show how this is equivalent as *shaping* the sensitivity function:
\begin{align}
                  & \left\| W_s(s) S(s) \right\|_\infty < 1\nonumber \\
  \Leftrightarrow & \left| W_s(j\omega) S(j\omega) \right| < 1 \quad \forall \omega\nonumber \\
  \Leftrightarrow & \left| S(j\omega) \right| < \frac{1}{\left| W_s(j\omega) \right|} \quad \forall \omega \label{eq:sensitivity_shaping}
\end{align}

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

  Therefore, the choice of the weighting function $W_s(s)$ is very important: its inverse magnitude will define the wanted *upper bound* of the sensitivity function magnitude.
#+end_important

#+begin_src latex :file loop_shaping_S_with_W.pdf
  \begin{tikzpicture}
    \node[block] (G) {$G(s)$};
    \node[addb={+}{-}{}{}{}, right=0.6 of G] (addw) {};
    \node[block, above right=1.0 and 1.0 of addw] (Ws) {$W_s(s)$};
    \coordinate[right=0.8 of Ws] (epsilon);

    \coordinate[] (w) at ($(epsilon-|G.west)+(-1.0, 0)$);

    \begin{scope}[on background layer]
      \node[fit={(G.south west) (Ws.north east)}, inner sep=8pt, draw, dashed, fill=black!20!white] (P) {};
      \node[above] at (P.north) {Weighted Generalized Plant $\tilde{P}(s)$};
    \end{scope}

    \node[block, below=0.4 of P] (K) {$K(s)$};

    % Connections
    \draw[->] (G.east) -- (addw.west);
    \draw[->] ($(addw.east)+(0.4, 0)$)node[branch]{} |- (Ws.west)node[above left]{$\epsilon$};
    \draw[->] (Ws.east) -- (epsilon) node[above](z1){$z = \tilde{\epsilon}$};

    \draw[->] (addw.east) -- (addw-|z1) |- node[near start, right]{$v$} (K.east);
    \draw[->] (K.west) -| node[near end, left]{$u$} ($(G-|w)+(0.2, 0)$) -- (G.west);

    \draw[->] (w) node[above]{$w = r$} -| (addw.north);
  \end{tikzpicture}
#+end_src

#+name: fig:loop_shaping_S_with_W
#+caption: Weighted Generalized Plant
#+RESULTS:
[[file:figs/loop_shaping_S_with_W.png]]

#+begin_exercice
Using matlab, compute the weighted generalized plant shown in Figure [[fig:first_order_weight]] as a function of $G(s)$ and $W_S(s)$.

#+HTML: <details><summary>Hint</summary>
The weighted generalized plant can be defined in Matlab using two techniques:
- by writing manually the 4 transfer functions from $[w, u]$ to $[\tilde{\epsilon}, v]$
- by pre-multiplying the (non-weighted) generalized plant by a block-diagonal transfer function matrix containing the weights for the outputs $z$ and =1= for the outputs $v$
#+HTML: </details>

#+HTML: <details><summary>Answer</summary>
The two solutions below can be used.

#+begin_src matlab :tangle no :eval no
  Pw = [Ws -Ws*G;
        1  -G];
#+end_src

#+begin_src matlab :tangle no :eval no
  Pw = blkdiag(Ws, 1)*P;
#+end_src

The second solution is however more general, and can also be used when weights are added at the inputs by post-multiplying instead of pre-multiplying.
#+HTML: </details>
#+end_exercice

** Design of Weighting Functions
<<sec:weighting_functions_design>>

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

#+begin_definition
- proper ::
  more poles than zeros, this implies $\lim_{\omega \to \infty} |W(j\omega)| < \infty$
- stable ::
  no poles in the right half plane
- minimum phase ::
  no zeros in the right half plane
#+end_definition

Matlab is providing the =makeweight= function that allows to design first-order weights by specifying the low frequency gain, high frequency gain, and the gain at a specific frequency:
#+begin_src matlab :tangle no :eval no
  W = makeweight(dcgain,[freq,mag],hfgain)
#+end_src
with:
- =dcgain=: low frequency gain
- =[freq,mag]=: frequency =freq= at which the gain is =mag=
- =hfgain=: high frequency gain

#+begin_exampl
The Matlab code below produces a weighting function with the following characteristics (Figure [[fig:first_order_weight]]):
- Low frequency gain of 100
- Gain of 1 at 10Hz
- High frequency gain of 0.5

#+begin_src matlab
  Ws = makeweight(1e2, [2*pi*10, 1], 1/2);
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(Ws, freqs, 'Hz'))), 'k-');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frquency [Hz]'); ylabel('Magnitude');
  hold off;
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/first_order_weight.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:first_order_weight
#+caption: Obtained Magnitude of the Weighting Function
#+RESULTS:
[[file:figs/first_order_weight.png]]
#+end_exampl

#+begin_seealso
Quite often, higher orders weights are required.

In such case, the following formula can be used:

\begin{equation}
  W(s) = \left( \frac{
           \frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_\infty}\right)^{\frac{2}{n}}}} s + \left(\frac{G_0}{G_c}\right)^{\frac{1}{n}}
         }{
           \left(\frac{1}{G_\infty}\right)^{\frac{1}{n}} \frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_\infty}\right)^{\frac{2}{n}}}} s + \left(\frac{1}{G_c}\right)^{\frac{1}{n}}
         }\right)^n \label{eq:weight_formula_advanced}
\end{equation}

The parameters permit to specify:
- the low frequency gain: $G_0 = lim_{\omega \to 0} |W(j\omega)|$
- the high frequency gain: $G_\infty = lim_{\omega \to \infty} |W(j\omega)|$
- the absolute gain at $\omega_0$: $G_c = |W(j\omega_0)|$
- the absolute slope between high and low frequency: $n$

A Matlab function implementing Equation eqref:eq:weight_formula_advanced is shown below:

#+name: lst:generateWeight
#+caption: Matlab Function that can be used to generate Weighting functions
#+begin_src matlab :tangle matlab/generateWeight.m :comments none :eval no
  function [W] = generateWeight(args)
      arguments
          args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1
          args.G1 (1,1) double {mustBeNumeric, mustBePositive} = 10
          args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1
          args.wc (1,1) double {mustBeNumeric, mustBePositive} = 2*pi
          args.n  (1,1) double {mustBeInteger, mustBePositive} = 1
      end

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

      s = zpk('s');

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

  end
#+end_src

Let's use this function to generate three weights with the same high and low frequency gains, but but different slopes.

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

The obtained shapes are shown in Figure [[fig:high_order_weight]].

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(W1, freqs, 'Hz'))), ...
       'DisplayName', '$n = 1$');
  plot(freqs, abs(squeeze(freqresp(W2, freqs, 'Hz'))), ...
       'DisplayName', '$n = 2$');
  plot(freqs, abs(squeeze(freqresp(W3, freqs, 'Hz'))), ...
       'DisplayName', '$n = 3$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frquency [Hz]'); ylabel('Magnitude');
  legend('location', 'northeast');
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/high_order_weight.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:high_order_weight
#+caption: Higher order weights using Equation eqref:eq:weight_formula_advanced
#+RESULTS:
[[file:figs/high_order_weight.png]]
#+end_seealso

** Shaping the Sensitivity Function
<<sec:sensitivity_shaping_example>>

Let's design a controller using the $\mathcal{H}_\infty$ synthesis that fulfils the following requirements:
1. Bandwidth of at least 10Hz
2. Small static errors for step responses
3. Robustness: Large module margin $\Delta M > 0.5$ ($\Rightarrow \Delta G > 2$ and $\Delta \phi > 29^o$)

As usual, the plant used is the one presented in Section [[sec:example_system]].

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

#+HTML: <details><summary>Hint</summary>
The typical wanted upper bound of the sensitivity function is shown in Figure [[fig:h-infinity-spec-S-bis]].

More precisely:
1. Recall that the closed-loop bandwidth is defined as the frequency $|S(j\omega)|$ first crosses $1/\sqrt{2} = -3dB$ from below
2. For the small static error, -60dB is usually enough as other factors (measurement noise, disturbances) will anyhow limit the performances
3. Recall that the module margin is equal to the inverse of the $\mathcal{H}_\infty$ norm of the sensitivity function:
   \[ \Delta M = \frac{1}{\|S\|_\infty} \]

Remember that the wanted upper bound of the sensitivity function is defined by the *inverse* magnitude of the weight.

#+name: fig:h-infinity-spec-S-bis
#+caption: Typical wanted shape of the Sensitivity transfer function
[[file:figs/h-infinity-spec-S.png]]
#+HTML: </details>

#+HTML: <details><summary>Answer</summary>
1. $|W_s(j \cdot 2 \pi 10)| = \sqrt{2}$
2. $|W_s(j \cdot 0)| = 10^3$
3. $\|W_s\|_\infty = 0.5$

Using Matlab, such weighting function can be generated using the =makeweight= function as shown below:
#+begin_src matlab :eval no :tangle no
  Ws = makeweight(1e3, [2*pi*10, sqrt(2)], 1/2);
#+end_src

Or using the =generateWeight= function:
#+begin_src matlab :eval no :tangle no
  Ws = generateWeight('G0', 1e3, ...
                      'G1', 1/2, ...
                      'Gc', sqrt(2), 'wc', 2*pi*10, ...
                      'n', 2);
#+end_src
#+HTML: </details>
#+end_exercice

Let's say we came up with the following weighting function:
#+begin_src matlab
  Ws = generateWeight('G0', 1e3, ...
                      'G1', 1/2, ...
                      'Gc', sqrt(2), 'wc', 2*pi*10, ...
                      'n', 2);
#+end_src

The weighting function is then added to the generalized plant.
#+begin_src matlab
  P = [1 -G;
       1 -G];
  Pw = blkdiag(Ws, 1)*P;
#+end_src

And the $\mathcal{H}_\infty$ synthesis is performed on the /weighted/ generalized plant.
#+begin_src matlab :results output replace
  K = hinfsyn(Pw, 1, 1, 'Display', 'on');
#+end_src

#+RESULTS:
#+begin_example
  Test bounds:  0.5 <=  gamma  <=  0.51

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

  Best performance (actual): 0.503
#+end_example

$\gamma \approx 0.5$ means that the $\mathcal{H}_\infty$ synthesis generated a controller $K(s)$ that stabilizes the closed-loop system, and such that:
\begin{aligned}
  & \| W_s(s) S(s) \|_\infty \approx 0.5 \\
  & \Leftrightarrow |S(j\omega)| < \frac{0.5}{|W_s(j\omega)|} \quad \forall \omega
\end{aligned}

This is indeed what we can see by comparing $|S|$ and $|W_S|$ in Figure [[fig:results_sensitivity_hinf]].

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

Having $\gamma$ slightly above one does not necessary means the obtained controller is not "good".
It just means that at some frequency, one of the closed-loop transfer functions is above the specified upper bound by a factor $\gamma$.
#+end_important

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(freqs, 1./abs(squeeze(freqresp(Ws, freqs, 'Hz'))), 'k--', 'DisplayName', '$|W_s|^{-1}$');
  plot(freqs, abs(squeeze(freqresp(1/(1 + K*G), freqs, 'Hz'))), 'k-', 'DisplayName', '$|S|$');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Magnitude');
  legend('location', 'southeast', 'FontSize', 8);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/results_sensitivity_hinf.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:results_sensitivity_hinf
#+caption: Weighting function and obtained closed-loop sensitivity
#+RESULTS:
[[file:figs/results_sensitivity_hinf.png]]

** Shaping multiple closed-loop transfer functions
<<sec:shaping_multiple_tf>>

As was shown in Section [[sec:modern_interpretation_specification]], depending on the specifications, up to four closed-loop transfer function may be shaped (the Gang of four).
This was summarized in Table [[tab:specification_modern]].

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

When multiple closed-loop transfer function are shaped at the same time, it is refereed to as "Mixed-Sensitivity $\mathcal{H}_\infty$ Control" and is the subject of Section [[sec:h_infinity_mixed_sensitivity]].

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

*** S KS                                                            :ignore:
#+HTML: <details><summary>Shaping of S and KS</summary>
#+begin_src latex :file general_conf_shaping_S_KS.pdf
  \begin{tikzpicture}
    % Blocs
    \node[block] (G) {$G$};

    \node[addb={+}{-}{}{}{}, right=0.6 of G] (addw) {};
    \node[block, above right=0.4 and 0.8 of addw] (W2) {$W_2$};
    \node[block, above=0.5 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.8 and 1.4 of G] (w);

    % Connections
    \draw[->] (G.east) -- (addw.west);
    \draw[->] ($(addw.east)+(0.2, 0)$)node[branch]{} |- (W1.west);
    \draw[->] (Gin)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[->] (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:general_conf_shaping_S_KS
#+caption: Generalized Plant to shape $S$ and $KS$
#+RESULTS:
[[file:figs/general_conf_shaping_S_KS.png]]

#+name: lst:general_plant_S_KS
#+caption: General Plant definition corresponding to Figure [[fig:general_conf_shaping_S_KS]]
#+begin_src matlab :eval no :tangle no
  P = [W1 -G*W1
       0   W2
       1  -G];
#+end_src

- $W_1(s)$ is used to shape $S$
- $W_2(s)$ is used to shape $KS$
#+HTML: </details>

*** S T                                                             :ignore:
#+HTML: <details><summary>Shaping of S and T</summary>
#+begin_src latex :file general_conf_shaping_S_T.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] (W2) {$W_2$};
    \node[block, above=0.5 of W2] (W1) {$W_1$};

    \begin{scope}[on background layer]
      \node[fit={(G.south west) (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.8 and 0.8 of G] (w);

    % Connections
    \draw[->] (G.east) -- (addw.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[->] (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:general_conf_shaping_S_T
#+caption: Generalized Plant to shape $S$ and $T$
#+RESULTS:
[[file:figs/general_conf_shaping_S_T.png]]

#+name: lst:general_plant_S_T
#+caption: General Plant definition corresponding to Figure [[fig:general_conf_shaping_S_T]]
#+begin_src matlab :eval no :tangle no
  P = [W1 -G*W1
       0   G*W2
       1   -G];
#+end_src

- $W_1$ is used to shape $S$
- $W_2$ is used to shape $T$
#+HTML: </details>

*** S GS                                                            :ignore:
#+HTML: <details><summary>Shaping of S and GS</summary>
#+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: </details>

*** S T KS                                                          :ignore:
#+HTML: <details><summary>Shaping of S, T and KS</summary>
#+begin_src latex :file general_conf_shaping_S_T_KS.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:general_conf_shaping_S_T_KS
#+caption: Generalized Plant to shape $S$, $T$ and $KS$
#+RESULTS:
[[file:figs/general_conf_shaping_S_T_KS.png]]

#+name: lst:general_plant_S_T_KS
#+caption: General Plant definition corresponding to Figure [[fig:general_conf_shaping_S_T_KS]]
#+begin_src matlab :eval no :tangle no
  P = [W1 -G*W1
       0   W2
       0   G*W3
       1   -G];
#+end_src

- $W_1$ is used to shape $S$
- $W_2$ is used to shape $KS$
- $W_3$ is used to shape $T$
#+HTML: </details>

*** S T GS                                                          :ignore:
#+HTML: <details><summary>Shaping of S, T and GS</summary>
#+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: </details>

*** S T KS GS                                                       :ignore:
#+HTML: <details><summary>Shaping of S, T, KS and GS</summary>
#+begin_src latex :file general_conf_shaping_S_T_KS_GS.pdf
  \begin{tikzpicture}
    % Blocs
    \node[block] (G) {$G$};

    \node[addb={+}{-}{}{}{}, right=0.6 of G] (addr) {};
    \node[addb, left=0.6 of G] (addd) {};
    \node[block, above right=0.4 and 0.8 of addr] (W2) {$W_2$};
    \node[block, above=0.5 of W2] (W1) {$W_1$};

    \node[block, above left=0.7 and 0.8 of addd] (W3) {$W_3$};
    \node[block, above=0.5 of W3] (W4) {$1$};

    \begin{scope}[on background layer]
      \node[fit={(W3.west|-G.south) (W4.north -| W2.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=1.0 of W1] (z);
    \coordinate[left=1.0 of W3] (w);

    % Connections
    \draw[->] (G.east) -- (addr.west);
    \draw[->] ($(addr.east)+( 0.2, 0)$)node[branch]{} |- (W1.west);
    \draw[->] ($(addd.west)+(-0.4, 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[->] (addr.east) -- (addw-|z) |- node[near start, right]{$v$} (K.east);
    \draw[->] (K.west) -| node[near end, left]{$u$} (addd-|w) -- (addd.west);
    \draw[->] (addd.east) -- (G.west);

    \draw[->] (W3.east) -| (addd.north);
    \draw[->] (W4.east) -| (addr.north);

    \draw[->] (w) node[above right]{$w_1$} -- (W3.west);
    \draw[->] (w|-W4) node[above right]{$w_2$} -- (W4.west);
  \end{tikzpicture}
#+end_src

#+name: fig:general_conf_shaping_S_T_KS_GS
#+caption: Generalized Plant to shape $S$, $T$, $KS$ and $GS$
#+RESULTS:
[[file:figs/general_conf_shaping_S_T_KS_GS.png]]

#+name: lst:general_plant_S_T_KS_GS
#+caption: General Plant definition corresponding to Figure [[fig:general_conf_shaping_S_T_KS_GS]]
#+begin_src matlab :eval no :tangle no
  P = [ W1  -W1*G*W3 -G*W1
        0    0        W2
        1   -G*W3    -G];
#+end_src

- $W_1$ is used to shape $S$
- $W_2$ is used to shape $KS$
- $W_1W_3$ is used to shape $GS$
- $W_2W_3$ is used to shape $T$
#+HTML: </details>



*** Limitation                                                      :ignore:

When shaping multiple closed-loop transfer functions, one should be verify careful about the three following points that are further discussed:
- The shaped closed-loop transfer functions are linked by mathematical relations and cannot be shaped
- Closed-loop transfer function can only be shaped in certain frequency range.
- The size of the obtained controller may be very large and not implementable in practice



#+begin_warning
  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:
  \begin{equation}
    S(s) + T(s) = 1
  \end{equation}

  This means that $|S(j\omega)|$ and $|T(j\omega)|$ cannot be made small at the same time!

  It is therefore *not* possible to shape the four closed-loop transfer functions independently.
  The weighting function should be carefully design such as these fundamental relations are not violated.
#+end_warning

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

\begin{align*}
  &|G(j\omega) K(j\omega)| \ll 1 \Longrightarrow |S(j\omega)| = \frac{1}{1 + |G(j\omega)K(j\omega)|} \approx 1 \\
  &|G(j\omega) K(j\omega)| \gg 1 \Longrightarrow |S(j\omega)| = \frac{1}{1 + |G(j\omega)K(j\omega)|} \approx \frac{1}{|G(j\omega)K(j\omega)|}
\end{align*}

Similar relationship can be found for $T$, $KS$ and $GS$.

#+begin_exercice
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$).

#+HTML: <details><summary>Hint</summary>
You can follows this procedure for $T$, $KS$ and $GS$:
1. Write the closed-loop transfer function $T(s)$ as a function of $K(s)$ and $G(s)$
2. Take $|K(j\omega)G(j\omega)| \gg 1$ and conclude on $|T(j\omega)|$
3. Take $|K(j\omega)G(j\omega)| \ll 1$ and conclude on $|T(j\omega)|$
#+HTML: </details>

#+HTML: <details><summary>Answer</summary>
The obtained constrains are shown in Figure [[fig:h-infinity-4-blocs-constrains]].
#+HTML: </details>
#+end_exercice

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

Therefore the weighting functions should only focus on certainty frequency range depending on the transfer function being shaped.
These regions are summarized in Figure [[fig:h-infinity-4-blocs-constrains]].

#+begin_src latex :file h-infinity-4-blocs-constrains.pdf
  \begin{tikzpicture}
    \begin{scope}[shift={(0, 0)}]
      \draw[fill=blue!20] (-0.2, -2.5) rectangle (1.4, 0.5);
      \draw[] (0.6, -0.5) node[]{$\sim GK^{-1}$};
      \draw[fill=red!20] (3.6, -2.5) rectangle (5.2, 0.5);
      \draw[] (4.5, -0.5) node[]{$\sim 1$};
      \draw[fill=red!20] (2.5, 0.15) circle (0.15);
      \draw[dashed] (-0.4, 0) -- (5.4, 0);
      \draw [] (0,-2) to[out=45,in=180+45] (2,0) to[out=45,in=180] (2.5,0.3) to[out=0,in=180] (3.5,0) to[out=0,in=180] (5, 0);
      \draw[dashed]  rectangle ;
      \begin{scope}[on background layer]
        \node[fit={(-0.5, -2.7) (5.5, 1.4)}, inner sep=0pt, draw, dashed, fill=black!20!white] (S) {};
        \node[below] at (S.north) {$S$};
      \end{scope}
    \end{scope}

    \begin{scope}[shift={(6.4, 0)}]
      \draw[fill=blue!20] (-0.2, -2.5) rectangle (1.4, 0.5);
      \draw[] (0.6, -0.5) node[]{$\sim K^{-1}$};
      \draw[fill=red!20] (3.6, -2.5) rectangle (5.2, 0.5);
      \draw[] (4.5, -0.5) node[]{$\sim G$};
      \draw[dashed] (-0.4, 0) -- (5.4, 0);
      \draw [] (0,-2) to[out=45,in=180+45] (1, -1) to[out=45, in=180] (2.5,-0.2) to[out=0,in=180-45] (4,-1) to[out=-45,in=180-45] (5, -2);
      \begin{scope}[on background layer]
        \node[fit={(-0.5, -2.7) (5.5, 1.4)}, inner sep=0pt, draw, dashed, fill=black!20!white] (GS) {};
        \node[below] at (GS.north) {$GS$};
      \end{scope}
    \end{scope}

    \begin{scope}[shift={(0, -4.4)}]
      \draw[fill=red!20] (-0.2, -2.5) rectangle (1.4, 0.5);
      \draw[] (0.6, -1.8) node[]{$\sim G^{-1}$};
      \draw[fill=blue!20] (3.6, -2.5) rectangle (5.2, 0.5);
      \draw[] (4.5, -0.3) node[]{$\sim K$};
      \draw[dashed] (-0.4, 0) -- (5.4, 0);
      \draw [] (0,-1.5) to[out=45,in=180+45] (1, -0.5) to[out=45, in=180] (2.5,0.3) to[out=0,in=180-45] (4,-0.5) to[out=-45,in=180-45] (5, -1.5);
      \begin{scope}[on background layer]
        \node[fit={(-0.5, -2.7) (5.5, 1.4)}, inner sep=0pt, draw, dashed, fill=black!20!white] (KS) {};
        \node[below] at (KS.north) {$KS$};
      \end{scope}
    \end{scope}

    \begin{scope}[shift={(6.4, -4.4)}]
      \draw[fill=red!20] (-0.2, -2.5) rectangle (1.4, 0.5);
      \draw[] (0.6, -0.5) node[]{$\sim 1$};
      \draw[fill=blue!20] (3.6, -2.5) rectangle (5.2, 0.5);
      \draw[] (4.5, -0.5) node[]{$\sim GK$};
      \draw[fill=red!20] (2.5, 0.15) circle (0.15);
      \draw[dashed] (-0.4, 0) -- (5.4, 0);
      \draw [] (0,0) to[out=0,in=180] (1.5,0) to[out=0,in=180] (2.5,0.3) to[out=0,in=-45] (3,0) to[out=-45,in=180-45] (5, -2);
      \begin{scope}[on background layer]
        \node[fit={(-0.5, -2.7) (5.5, 1.4)}, inner sep=0pt, draw, dashed, fill=black!20!white] (T) {};
        \node[below] at (T.north) {$T$};
      \end{scope}
    \end{scope}
  \end{tikzpicture}
#+end_src

#+name: fig:h-infinity-4-blocs-constrains
#+caption: Shaping the Gang of Four: Limitations
#+RESULTS:
[[file:figs/h-infinity-4-blocs-constrains.png]]

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

  Two approaches can be used to obtain controllers with reasonable order:
  1. use simple weights (usually first order)
  2. perform a model reduction on the obtained high order controller
#+end_warning

* Mixed-Sensitivity $\mathcal{H}_\infty$ Control - Example
<<sec:h_infinity_mixed_sensitivity>>

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


#+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$
<<sec:step_shaping_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$
<<sec:step_shaping_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$
<<sec:step_shaping_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$
<<sec:step_shaping_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$
<<sec:step_shaping_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$
<<sec:step_shaping_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$
<<sec:step_shaping_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
<<sec:conclusion>>

* Resources
:PROPERTIES:
:UNNUMBERED: notoc
:END:

yt:?listType=playlist&list=PLn8PRpmsu08qFLMfgTEzR8DxOPE7fBiin

yt:?listType=playlist&list=PLsjPUqcL7ZIFHCObUU_9xPUImZ203gB4o