lecture-h-infinity/index.org

1451 lines
58 KiB
Org Mode
Raw Normal View History

2020-11-30 17:44:13 +01:00
#+TITLE: A brief and practical introduction to $\mathcal{H}_\infty$ Control
2020-11-25 19:35:11 +01:00
:DRAWER:
#+STARTUP: overview
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas
#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html
2020-11-25 19:37:47 +01:00
#+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>
2020-11-25 19:35:11 +01:00
#+HTML_MATHJAX: align: center tagside: right font: TeX
#+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
2020-11-30 17:44:13 +01:00
#+PROPERTY: header-args:latex+ :tangle no
2020-11-25 19:35:11 +01:00
#+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:
* TODO Introduction :ignore:
2020-11-25 19:35:11 +01:00
* 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
2020-11-30 17:44:13 +01:00
#+begin_src matlab :tangle no
addpath('matlab')
#+end_src
2020-11-25 19:35:11 +01:00
* Introduction to the Control Methodology - Model Based Control
2020-11-30 17:44:13 +01:00
<<sec:model_based_control>>
** Model Based Control - Methodology
<<sec:model_based_control_methodology>>
2020-11-25 19:35:11 +01:00
The typical methodology when applying Model Based Control to a plant is schematically shown in Figure [[fig:control-procedure]].
It consists of three steps:
1. *Identification or modeling*: $\Longrightarrow$ mathematical model
2. *Translate the specifications into mathematical criteria*:
- _Specifications_: Response Time, Noise Rejection, Maximum input amplitude, Robustness, ...
- _Mathematical Criteria_: Cost Function, Shape of TF
# - Cost Function, Needed Bandwidth, Roll-off, ...
2020-11-25 19:38:32 +01:00
# - $\Longrightarrow$ We will use the $\mathcal{H}_\infty$ Norm
2020-11-25 19:35:11 +01:00
3. *Synthesis*: research of $K$ 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 mainly focus on steps 2 and 3.
2020-11-27 18:26:06 +01:00
** Some Background: From Classical Control to Robust Control
2020-11-30 17:44:13 +01:00
<<sec:comp_classical_modern_robust_control>>
2020-11-25 19:35:11 +01:00
#+name: tab:comparison_control_methods
#+caption: Table summurazing the main differences between classical, modern and robust control
| | *Classical Control* | *Modern Control* | *Robust Control* |
| <l> | <c> | <c> | <c> |
|-------------------------+------------------------------------+--------------------------------------+-------------------------------------------------------------------------|
| *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 |
|-------------------------+------------------------------------+--------------------------------------+-------------------------------------------------------------------------|
| *Control Architectures* | Proportional, Integral, Derivative | Full State Feedback | General Control Configuration |
| | Leads, Lags | LQR, LQG | |
| | | Kalman Filters | |
|-------------------------+------------------------------------+--------------------------------------+-------------------------------------------------------------------------|
| *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 |
2020-11-27 18:26:06 +01:00
** Example System
2020-11-30 17:44:13 +01:00
<<sec:example_system>>
2020-11-27 18:26:06 +01:00
Let's consider the model shown in Figure [[fig:mech_sys_1dof_inertial_contr]].
It could represent a suspension system with a payload to position or isolate using an force actuator and an inertial sensor.
2020-11-27 18:26:06 +01:00
The notations used are listed 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)] |
2020-11-27 18:26:06 +01:00
| $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:
2020-11-27 18:26:06 +01:00
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>
2020-11-27 18:26:06 +01:00
#+end_exercice
2020-11-30 17:44:13 +01:00
Hi Musa,
Thank you very much for sharing this awesome package.
For a long time, I am dreaming of being abble to export source blocks to HTML tha are surounded by <details> blocks.
For now, I am manually adding #+HTML: <details><summary>Code</summary> and #+HTML: </details> around the source blocks I want to hide...
This is a very simple solution, but not so elegent nor practical.
Do you have any idea if it would be easy to extend to org-mode export of source blocks to add such functionallity?
Similarly, I would love to be able to export a <span> block with the name of the file corresponding to the source block.
For instance, if a particular source block is tangled to script.sh, it would be so nice to display the filename when exporting!
Thanks in advance
2020-11-27 18:26:06 +01:00
Having obtained $G(s)$ and $G_d(s)$, we can transform the system shown in Figure [[fig:mech_sys_1dof_inertial_contr]] into a classical feedback form as shown in Figure [[fig:open_loop_shaping]].
#+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
#+RESULTS:
[[file:figs/classical_feedback_test_system.png]]
Let's define the system parameters on Matlab.
#+begin_src matlab +n
2020-11-27 18:26:06 +01:00
k = 1e6; % Stiffness [N/m]
c = 4e2; % Damping [N/(m/s)]
m = 10; % Mass [kg]
2020-11-27 18:26:06 +01:00
#+end_src
And now the system dynamics $G(s)$ and $G_d(s)$ (their bode plots are shown in Figures [[fig:bode_plot_example_afm]] and [[fig:bode_plot_example_Gd]]).
#+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
2020-11-27 18:26:06 +01:00
#+end_src
#+begin_src matlab :exports none
freqs = logspace(0, 3, 1000);
2020-11-27 18:26:06 +01:00
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;
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
2020-11-30 17:44:13 +01:00
<<sec:open_loop_shaping>>
** Introduction to Loop Shaping
<<sec:open_loop_shaping_introduction>>
#+begin_definition
2020-11-30 17:44:13 +01:00
*Loop Shaping* refers to a design procedure that involves explicitly shaping the magnitude of the *Loop Transfer Function* $L(s)$.
#+end_definition
#+begin_definition
The *Loop Gain* $L(s)$ usually refers to as the product of the controller and the plant ("Gain around the loop", see Figure [[fig:open_loop_shaping]]):
\begin{equation}
L(s) = G(s) \cdot K(s) \label{eq:loop_gain}
\end{equation}
#+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
2020-11-30 17:44:13 +01:00
This synthesis method is widely used as many characteristics of the closed-loop system depend on the shape of the open loop gain $L(s)$ such as:
- *Performance*: $L$ large
- *Good disturbance rejection*: $L$ large
- *Limitation of measurement noise on plant output*: $L$ small
- *Small magnitude of input signal*: $K$ and $L$ small
- *Nominal stability*: $L$ small (RHP zeros and time delays)
- *Robust stability*: $L$ small (neglected dynamics)
The Open Loop shape is usually done manually has the loop gain $L(s)$ depends linearly on $K(s)$ eqref:eq:loop_gain.
2020-11-30 17:44:13 +01:00
$K(s)$ then consists of a combination of leads, lags, notches, etc. such that $L(s)$ has the wanted shape (an example is shown in Figure [[fig:open_loop_shaping_shape]]).
#+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]{Reference\\Tracking};
\draw[<-] (hfshaping) -- ++(0, 0.8) node[above, align=center]{Noise\\Rejection};
% Crossover frequency
\node[below] (wc) at (4,2){$\omega_c$};
\draw[<-] (wc.south) -- ++(0, -0.4) node[below, align=center]{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]{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]]
** Example of Open Loop Shaping
2020-11-30 17:44:13 +01:00
<<sec:loop_shaping_example>>
#+begin_exampl
Let's take our example system and try to apply the Open-Loop shaping strategy to design a controller that fulfils the following specifications:
- *Performance*: Bandwidth of approximately 10Hz
- *Noise Attenuation*: Roll-off of -40dB/decade past 30Hz
- *Robustness*: Gain margin > 3dB and Phase margin > 30 deg
#+end_exampl
2020-11-27 18:26:06 +01:00
#+begin_exercice
Using =SISOTOOL=, design a controller that fulfill the specifications.
#+begin_src matlab :eval no
sisotool(G)
#+end_src
#+end_exercice
In order to have the wanted Roll-off, two integrators are used, a lead is also added to have sufficient phase margin.
The obtained controller is shown below, and the bode plot of the Loop Gain is shown in Figure [[fig:loop_gain_manual_afm]].
#+begin_src matlab
K = 14e8 * ... % Gain
2020-11-27 18:26:06 +01:00
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
2020-11-27 18:26:06 +01:00
#+end_src
#+begin_src matlab :exports none
freqs = logspace(0, 3, 1000);
2020-11-27 18:26:06 +01:00
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])
2020-11-27 18:26:06 +01:00
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');
2020-11-27 18:26:06 +01:00
#+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]]
And we can verify that we have the wanted stability margins:
#+begin_src matlab
2020-11-27 18:26:06 +01:00
[Gm, Pm, ~, Wc] = margin(G*K)
#+end_src
#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
2020-11-30 17:44:13 +01:00
data2orgtable([Gm; Pm; Wc/2/pi], {'Gain Margin $> 3$ [dB]', 'Phase Margin $> 30$ [deg]', 'Crossover $\approx 10$ [Hz]'}, {'Requirements', 'Manual Method'}, ' %.1f ');
2020-11-27 18:26:06 +01:00
#+end_src
#+RESULTS:
2020-11-30 17:44:13 +01:00
| Requirements | Manual Method |
|-----------------------------+---------------|
| Gain Margin $> 3$ [dB] | 3.1 |
| Phase Margin $> 30$ [deg] | 35.4 |
| Crossover $\approx 10$ [Hz] | 10.1 |
2020-11-27 18:26:06 +01:00
** $\mathcal{H}_\infty$ Loop Shaping Synthesis
2020-11-30 17:44:13 +01:00
<<sec:h_infinity_open_loop_shaping>>
2020-11-27 18:26:06 +01:00
The Open Loop Shaping synthesis can be performed using the $\mathcal{H}_\infty$ Synthesis.
Even though we will not go into details, we will provide one example.
2020-11-30 17:44:13 +01:00
Using Matlab, the $\mathcal{H}_\infty$ Loop Shaping Synthesis can be performed using the =loopsyn= command:
2020-11-27 18:26:06 +01:00
#+begin_src matlab :eval no
K = loopsyn(G, Gd);
#+end_src
where:
- =G= is the (LTI) plant
- =Gd= 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
** Example of the $\mathcal{H}_\infty$ Loop Shaping Synthesis
2020-11-30 17:44:13 +01:00
<<sec:h_infinity_open_loop_shaping_example>>
2020-11-27 18:26:06 +01:00
Let's reuse the previous plant.
2020-11-27 18:26:06 +01:00
Translate the specification into the wanted shape of the open loop gain.
- *Performance*: Bandwidth of approximately 10Hz: $|L_w(j2 \pi 10)| = 1$
- *Noise Attenuation*: Roll-off of -40dB/decade past 30Hz
- *Robustness*: Gain margin > 3dB and Phase margin > 30 deg
2020-11-27 18:26:06 +01:00
#+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
2020-11-27 18:26:06 +01:00
2020-11-30 17:44:13 +01:00
The $\mathcal{H}_\infty$ optimal open loop shaping synthesis is performed using the =loopsyn= command:
#+begin_src matlab
[K, ~, GAM] = loopsyn(G, Lw);
2020-11-27 18:26:06 +01:00
#+end_src
The Bode plot of the obtained controller is shown in Figure [[fig:open_loop_shaping_hinf_K]].
#+begin_important
It is always important to analyze the controller after the synthesis is performed.
In the end, a synthesize controller is just a combination of low pass filters, high pass filters, notches, leads, etc.
#+end_important
Let's briefly analyze this controller:
- two integrators are used at low frequency to have the wanted low frequency high gain
- a lead is added centered with the crossover frequency to increase the phase margin
- a notch is added at the resonance of the plant to increase the gain margin (this is very typical of $\mathcal{H}_\infty$ controllers, and can be an issue, more info on that latter)
2020-11-27 18:26:06 +01:00
#+begin_src matlab :exports none
freqs = logspace(0, 3, 1000);
2020-11-27 18:26:06 +01:00
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)]);
2020-11-27 18:26:06 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/open_loop_shaping_hinf_K.pdf', 'width', 'wide', 'height', 'normal');
2020-11-27 18:26:06 +01:00
#+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]]
The obtained Loop Gain is shown in Figure [[fig:open_loop_shaping_hinf_L]].
2020-11-27 18:26:06 +01:00
#+begin_src matlab :exports none
freqs = logspace(0, 3, 1000);
2020-11-27 18:26:06 +01:00
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');
2020-11-27 18:26:06 +01:00
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
hold off;
legend('location', 'northeast');
ylim([1e-4, 1e2]);
2020-11-27 18:26:06 +01:00
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]]
Let's now 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 |
2020-11-27 18:26:06 +01:00
2020-11-30 17:44:13 +01:00
* First Steps in the $\mathcal{H}_\infty$ world
<<sec:h_infinity_introduction>>
** The $\mathcal{H}_\infty$ Norm
2020-11-30 17:44:13 +01:00
<<sec:h_infinity_norm>>
2020-11-25 19:35:11 +01:00
#+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
2020-11-30 17:44:13 +01:00
Let's compute the $\mathcal{H}_\infty$ norm of our test plant $G(s)$ using the =hinfnorm= function:
2020-11-25 19:35:11 +01:00
#+begin_src matlab :results value replace
hinfnorm(G)
#+end_src
#+RESULTS:
: 7.9216e-06
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
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]].
2020-11-25 19:35:11 +01:00
#+begin_src matlab :exports none
freqs = logspace(0, 3, 1000);
2020-11-25 19:35:11 +01:00
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$')
2020-11-25 19:35:11 +01:00
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
2020-11-30 17:44:13 +01:00
<<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.
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
Then applying the $\mathcal{H}_\infty$ synthesis to a plant, the engineer work usually consists of the following steps
2020-11-25 19:38:32 +01:00
1. Write the problem as standard $\mathcal{H}_\infty$ problem
2020-11-30 17:44:13 +01:00
2. Translate the specifications as $\mathcal{H}_\infty$ norms of transfer functions
2020-11-25 19:35:11 +01:00
3. Make the synthesis and analyze the obtain controller
4. Reduce the order of the controller for implementation
2020-11-30 17:44:13 +01:00
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]])
2020-11-25 19:38:32 +01:00
- Signal Based $\mathcal{H}_\infty$ Synthesis
2020-11-25 19:35:11 +01:00
** The Generalized Plant
2020-11-30 17:44:13 +01:00
<<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}
2020-11-25 19:35:11 +01:00
#+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
2020-11-30 17:44:13 +01:00
#+name: fig:general_plant
#+caption: Inputs and Outputs of the generalized Plant
2020-11-25 19:35:11 +01:00
#+RESULTS:
[[file:figs/general_plant.png]]
2020-11-30 17:44:13 +01:00
#+begin_important
2020-11-25 19:35:11 +01:00
#+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 |
2020-11-30 17:44:13 +01:00
#+end_important
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
** The General Synthesis Problem Formulation
<<sec:h_infinity_general_synthesis>>
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
Once the generalized plant is obtained, the $\mathcal{H}_\infty$ synthesis problem can be stated as follows:
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
#+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
2020-11-25 19:35:11 +01:00
\begin{tikzpicture}
2020-11-30 17:44:13 +01:00
% 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);
2020-11-25 19:35:11 +01:00
\end{tikzpicture}
#+end_src
2020-11-30 17:44:13 +01:00
#+name: fig:general_control_names
#+caption: General Control Configuration
2020-11-25 19:35:11 +01:00
#+RESULTS:
2020-11-30 17:44:13 +01:00
[[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
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
2020-11-25 19:35:11 +01:00
#+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)$};
2020-11-25 19:35:11 +01:00
\end{scope}
\end{tikzpicture}
#+end_src
#+begin_exercice
2020-11-30 17:44:13 +01:00
Compute the Generalized plant of corresponding to the tracking control architecture shown in Figure [[fig:classical_feedback_tracking]]
2020-11-25 19:35:11 +01:00
#+name: fig:classical_feedback_tracking
#+caption: Classical Feedback Control Architecture (Tracking)
[[file:figs/classical_feedback_tracking.png]]
2020-11-30 17:44:13 +01:00
#+HTML: <details><summary>Hint</summary>
2020-11-25 19:35:11 +01:00
First, define the signals of the generalized plant:
- Exogenous inputs: $w = r$
- Signals to be minimized: $z_1 = \epsilon$, $z_2 = u$
- Control signals: $v = y$
- Control inputs: $u$
Then, Remove $K$ and rearrange the inputs and outputs.
2020-11-30 17:44:13 +01:00
#+HTML: </details>
#+HTML: <details><summary>Answer</summary>
The obtained generalized plant shown in Figure [[fig:mixed_sensitivity_ref_tracking]].
2020-11-25 19:35:11 +01:00
#+name: fig:mixed_sensitivity_ref_tracking
#+caption: Generalized plant of the Classical Feedback Control Architecture (Tracking)
[[file:figs/mixed_sensitivity_ref_tracking.png]]
2020-11-30 17:44:13 +01:00
#+HTML: </details>
#+end_exercice
Using Matlab, the generalized plant can be defined as follows:
#+begin_src matlab :tangle no :eval no
P = [1 -G;
0 1;
1 -G]
P.InputName = {'w', 'u'};
P.OutputName = {'e', 'u', 'v'};
#+end_src
* Modern Interpretation of the Control Specifications
<<sec:modern_interpretation_specification>>
** Introduction
The
#+name: fig:gang_of_four_feedback
#+caption: Simple Feedback Architecture
#+RESULTS:
[[file:figs/gang_of_four_feedback.png]]
- *Reference tracking* Overshoot, Static error, Settling time
- From $r$ to $\epsilon$
- *Disturbances rejection*
- From $d$ to $y$
- $G(s) S(s) = T_{d \rightarrow \epsilon}$
- *Measurement noise filtering*
- From $n$ to $y$
- $T(s) = T_{n \rightarrow \epsilon}$
- *Small command amplitude*
- From $n, r, d$ to $u$
- $K(s) S(s) = T_{r \rightarrow u}$
- *Stability*
- $S(s)$, $T(s)$, $K(s)S(s)$, $G(s)S(s)$
- *Robustness to plant uncertainty* (stability margins)
- *Controller implementation*
** Closed Loop Transfer Functions
<<sec:losed_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_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]]
#+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>
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
#+HTML: <details><summary>Anwser</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
2020-11-30 17:44:13 +01:00
We can see that their are 4 different transfer functions describing the behavior of the system in Figure [[fig:gang_of_four_feedback]].
They are 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}
2020-11-30 17:44:13 +01:00
#+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
2020-11-30 17:44:13 +01:00
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}
** Sensitivity Transfer Function
<<sec:sensitivity_transfer_functions>>
#+begin_src matlab
K1 = 14e8 * ... % Gain
1/(s^2) * ... % Double Integrator
(1 + s/(2*pi*10/sqrt(8)))/(1 + s/(2*pi*10*sqrt(8))); % Lead
K2 = 1e8 * ... % Gain
1/(s^2) * ... % Double Integrator
(1 + s/(2*pi*1/sqrt(8)))/(1 + s/(2*pi*1*sqrt(8))); % Lead
K3 = 1e8 * ... % Gain
1/(s^2) * ... % Double Integrator
(1 + s/(2*pi*1/sqrt(2)))/(1 + s/(2*pi*1*sqrt(2))); % Lead
S1 = 1/(1 + K1*G);
S2 = 1/(1 + K2*G);
S3 = 1/(1 + K3*G);
T1 = K1*G/(1 + K1*G);
T2 = K2*G/(1 + K2*G);
T3 = K3*G/(1 + K3*G);
bodeFig({S1, S2, S3})
#+end_src
#+begin_src matlab
freqs = logspace(-1, 2, 1000);
figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');
ax1 = nexttile;
hold on;
plot(freqs, abs(squeeze(freqresp(S1, freqs, 'Hz'))), 'DisplayName', '$L(s)$');
plot(freqs, abs(squeeze(freqresp(S2, freqs, 'Hz'))), 'DisplayName', '$L_w(s)$');
plot(freqs, abs(squeeze(freqresp(S3, freqs, 'Hz'))), 'DisplayName', '$L_w(s) / \gamma$, $L_w(s) \cdot \gamma$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frquency [Hz]'); ylabel('Sensitivity Magnitude');
hold off;
ax2 = nexttile;
t = linspace(0, 1, 1000);
y1 = step(T1, t);
y2 = step(T2, t);
y3 = step(T3, t);
hold on;
plot(t, y1)
plot(t, y2)
plot(t, y3)
hold off
xlabel('Time [s]'); ylabel('Step Response');
#+end_src
#+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]{$\SI{-60}{\decibel} \rightarrow$ \footnotesize{Static error}};
\draw[<-] (1, -40) -- (10, -40) node[right, fill=white, draw]{$\SI{20}{\decibel/dec} \rightarrow$ \footnotesize{Ref. track.}};
\draw[<-] (100, 0) -- (3, 0) node[left, fill=white, draw]{$\omega_c \rightarrow$ \footnotesize{Speed}};
\draw[<-] (300, 6) -- (200, 20);
\draw (200, 20) -- (10, 20) node[left, fill=white, draw]{$\SI{6}{\decibel} \rightarrow$ \footnotesize{Module margin}};
\end{axis}
\end{tikzpicture}
#+end_src
#+name: fig:h-infinity-spec-S
#+caption: Typical wanted shape of the Sensitivity transfer function
#+RESULTS:
[[file:figs/h-infinity-spec-S.png]]
** Robustness: Module Margin
<<sec:module_margin>>
- [ ] Definition of Module margin
- [ ] Why it represents robustness
- [ ] Example
\[ M_S < 2 \Rightarrow \text{GM} > 2 \text{ and } \text{PM} > 29^o \]
** How to *Shape* transfer function? Using of Weighting Functions!
<<sec:weighting_functions>>
Let's say we want to shape the sensitivity transfer function corresponding to the transfer function from $r$ to $\epsilon$ of the control architecture shown in Figure [[fig:loop_shaping_S_without_W]].
#+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 left](z1){$\epsilon$};
\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);
\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]]
If the $\mathcal{H}_\infty$ synthesis is directly applied on the generalized plant $P(s)$ shown in Figure [[fig:loop_shaping_S_without_W]], if will minimize the $\mathcal{H}_\infty$ norm of transfer function from $r$ to $\epsilon$ (the sensitivity transfer function).
However, as the $\mathcal{H}_\infty$ norm is the maximum peak value of the transfer function's magnitude, it does not allow to *shape* the norm over all frequencies.
2020-11-30 17:44:13 +01:00
A /trick/ is to include a *weighting function* in the generalized plant as shown in Figure [[fig:loop_shaping_S_with_W]].
Applying the $\mathcal{H}_\infty$ synthesis to the /weighted/ generalized plant $\tilde{P}(s)$ (Figure [[fig:loop_shaping_S_with_W]]) will generate a controller $K(s)$ that minimizes the $\mathcal{H}_\infty$ norm between $r$ and $\tilde{\epsilon}$:
\begin{equation}
\begin{aligned}
& \left\| \frac{\tilde{\epsilon}}{r} \right\|_\infty < \gamma (=1) \\
\Leftrightarrow & \left\| W_s(s) S(s) \right\|_\infty < 1 \\
\Leftrightarrow & \left| W_s(j\omega) S(j\omega) \right| < 1 \quad \forall \omega \\
\Leftrightarrow & \left| S(j\omega) \right| < \frac{1}{\left| W_s(j\omega) \right|} \quad \forall \omega
\end{aligned}\label{eq:sensitivity_shaping}
\end{equation}
#+begin_important
As shown in Equation eqref:eq:sensitivity_shaping, the $\mathcal{H}_\infty$ synthesis 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 frequency dependent upper bound of the sensitivity transfer function magnitude.
#+end_important
2020-11-30 17:44:13 +01:00
#+begin_src latex :file loop_shaping_S_with_W.pdf
\begin{tikzpicture}
2020-11-30 17:44:13 +01:00
\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);
2020-11-30 17:44:13 +01:00
\coordinate[] (w) at ($(epsilon-|G.west)+(-1.0, 0)$);
2020-11-30 17:44:13 +01:00
\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}
2020-11-30 17:44:13 +01:00
\node[block, below=0.4 of P] (K) {$K(s)$};
2020-11-30 17:44:13 +01:00
% 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 left](z1){$\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.4, 0)$) -- (G.west);
\draw[->] (w) node[above]{$w = r$} -| (addw.north);
\end{tikzpicture}
#+end_src
2020-11-30 17:44:13 +01:00
#+name: fig:loop_shaping_S_with_W
#+caption: Weighted Generalized Plant
#+RESULTS:
2020-11-30 17:44:13 +01:00
[[file:figs/loop_shaping_S_with_W.png]]
2020-11-30 17:44:13 +01:00
Once the weighting function is designed, it should be added to the generalized plant as shown in Figure [[fig:loop_shaping_S_with_W]].
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
The weighted generalized plant can be defined in Matlab by either re-defining all the inputs or by pre-multiplying the (non-weighted) generalized plant by a block-diagonal MIMO transfer function containing the weights for the outputs $z$ and =1= for the outputs $v$.
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
#+begin_src matlab :tangle no :eval no
Pw = [Ws -Ws*G;
1 -G]
% Alternative
Pw = blkdiag(Ws, 1)*P;
#+end_src
** Design of Weighting Functions
<<sec:weighting_functions_design>>
Weighting function used must be *proper*, *stable* and *minimum phase* transfer functions.
- proper ::
more poles than zeros, this implies $\lim_{\omega \to \infty} |W(j\omega)| < \infty$
- stable ::
no poles in the right half plane
- minimum phase ::
no zeros in the right half plane
Matlab is providing the =makeweight= function that creates a first-order weights by specifying the low frequency gain, high frequency gain, and a gain at a specific frequency:
#+begin_src matlab :tangle no :eval no
W = makeweight(dcgain,[freq,mag],hfgain)
#+end_src
with:
- =dcgain=
- =freq=
- =mag=
- =hfgain=
#+begin_exampl
The Matlab code below produces a weighting function with a magnitude shape shown in Figure [[fig:first_order_weight]].
#+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 the design of these weights:
\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
** Sensitivity Function Shaping - Example
<<sec:sensitivity_shaping_example>>
- Robustness: Module margin > 2 ($\Rightarrow \text{GM} > 2 \text{ and } \text{PM} > 29^o$)
- Bandwidth:
- Slope of -2
First, the weighting functions is generated.
#+begin_src matlab
Ws = generateWeight('G0', 1e3, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 2);
#+end_src
It is then added to the generalized plant.
#+begin_src matlab
Pw = blkdiag(Ws, 1)*P;
#+end_src
And the $\mathcal{H}_\infty$ synthesis is performed.
#+begin_src matlab :results output replace
K = hinfsyn(Pw, 1, 1, 'Display', 'on');
#+end_src
#+RESULTS:
#+begin_example
K = hinfsyn(Pw, 1, 1, 'Display', 'on');
Test bounds: 0.5 <= gamma <= 0.51
gamma X>=0 Y>=0 rho(XY)<1 p/f
5.05e-01 0.0e+00 0.0e+00 4.497e-28 p
Limiting gains...
5.05e-01 0.0e+00 0.0e+00 0.000e+00 p
5.05e-01 -1.8e+01 # -2.9e-15 1.514e-15 f
Best performance (actual): 0.504
#+end_example
The obtained $\gamma \approx 0.5$ means that it found a controller $K(s)$ that stabilize the closed-loop system, and such that:
\begin{aligned}
& \| W_s(s) S(s) \|_\infty < 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_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]]
* $\mathcal{H}_\infty$ Mixed-Sensitivity Synthesis
<<sec:h_infinity_mixed_sensitivity>>
** Problem
** Typical Procedure
** Step 1 - Shaping of the Sensitivity Function
** Step 2 - Shaping of
2020-11-25 19:35:11 +01:00
2020-11-30 17:44:13 +01:00
* Conclusion
2020-11-25 19:35:11 +01:00
* Resources
yt:?listType=playlist&list=PLn8PRpmsu08qFLMfgTEzR8DxOPE7fBiin
yt:?listType=playlist&list=PLsjPUqcL7ZIFHCObUU_9xPUImZ203gB4o