dehaeze20_virtu_senso_fusio/matlab/index.html

968 lines
46 KiB
HTML
Raw Permalink Normal View History

2020-10-08 10:52:20 +02:00
<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
2020-11-12 10:44:04 +01:00
<!-- 2020-11-12 jeu. 10:43 -->
2020-10-08 10:52:20 +02:00
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
2020-11-12 10:44:04 +01:00
<title>Sensor Fusion Paper - Matlab Computation</title>
2020-10-08 10:52:20 +02:00
<meta name="generator" content="Org mode" />
<meta name="author" content="Thomas Dehaeze" />
2020-11-12 10:44:04 +01:00
<link rel="stylesheet" type="text/css" href="https://research.tdehaeze.xyz/css/style.css"/>
<script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.js"></script>
<script>MathJax = {
tex: {
tags: 'ams',
macros: {bm: ["\\boldsymbol{#1}",1],}
}
};
</script>
<script type="text/javascript" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
2020-10-08 10:52:20 +02:00
</head>
<body>
2020-11-12 10:44:04 +01:00
<div id="org-div-home-and-up">
<a accesskey="h" href="../index.html"> UP </a>
|
<a accesskey="H" href="../index.html"> HOME </a>
</div><div id="content">
<h1 class="title">Sensor Fusion Paper - Matlab Computation</h1>
2020-10-08 10:52:20 +02:00
<div id="table-of-contents">
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
2020-11-12 10:44:04 +01:00
<li><a href="#org917402b">1. Definition of the plant</a></li>
<li><a href="#orgde878e2">2. Multiplicative input uncertainty</a></li>
<li><a href="#orgc908813">3. Specifications and performance weights</a></li>
<li><a href="#orgce7235f">4. Upper bounds on the norm of the complementary filters for NP, RS and RP</a></li>
<li><a href="#orgaf8f3b1">5. H-Infinity synthesis of complementary filters</a></li>
<li><a href="#orga721ae7">6. Complementary filters using analytical formula</a></li>
<li><a href="#orgaf265db">7. Comparison of complementary filters</a></li>
<li><a href="#orga1bd48e">8. Controller Analysis</a></li>
<li><a href="#org1eb8c07">9. Nominal Stability and Nominal Performance</a></li>
<li><a href="#org8c33736">10. Robust Stability and Robust Performance</a></li>
<li><a href="#org367d76d">11. Pre-filter</a></li>
<li><a href="#org846e5c4">12. Matlab Functions</a>
<ul>
<li><a href="#org2f50732">12.1. <code>createWeight</code></a></li>
<li><a href="#org8ef5bfb">12.2. <code>plotMagUncertainty</code></a></li>
<li><a href="#org3b855e9">12.3. <code>plotPhaseUncertainty</code></a></li>
</ul>
</li>
2020-10-08 10:52:20 +02:00
</ul>
</div>
</div>
<p>
2020-11-12 10:44:04 +01:00
The control architecture studied here is shown on figure <a href="#org16aa33e">1</a> where:
2020-10-08 10:52:20 +02:00
</p>
<ul class="org-ul">
<li>\(G^{\prime}\) is the plant to control</li>
<li>\(K = G^{-1} H_L^{-1}\) is the controller used with \(G\) a model of the plant</li>
<li>\(H_L\) and \(H_H\) are complementary filters (\(H_L + H_H = 1\))</li>
<li>\(K_r\) is a pre-filter that can be added</li>
</ul>
2020-11-12 10:44:04 +01:00
<div id="org16aa33e" class="figure">
<p><img src="figs-paper/sf_arch_class_prefilter.png" alt="sf_arch_class_prefilter.png" />
2020-10-08 10:52:20 +02:00
</p>
<p><span class="figure-number">Figure 1: </span>Control Architecture</p>
</div>
<p>
Here is the outline of the <code>matlab</code> analysis for this control architecture:
</p>
<ul class="org-ul">
2020-11-12 10:44:04 +01:00
<li>Section <a href="#orgfbfc3f4">1</a>:
The plant model \(G\) is defined</li>
<li>Section <a href="#orge119055">2</a>:
The plant uncertainty set \(\Pi_I\) is defined using the multiplicative input uncertainty: \(\Pi_I: \ G^\prime = G (1 + w_I \Delta)\).
Thus the weight \(w_I\) is defined such that the true system dynamics is included in the set \(\Pi_I\)</li>
<li>Section <a href="#org3dfa664">3</a>:
From the specifications on performance that are expressed in terms of upper bounds of \(S\) and \(T\), performance weights \(w_S\) and \(w_T\) are derived such that the goal is to obtain \(|S| < \frac{1}{|w_S|}\) and \(|T| < \frac{1}{|w_T|}, \ \forall \omega\)</li>
<li>Section <a href="#orgacb23fe">4</a>:
Upper bounds on the magnitude of the complementary filters \(|H_L|\) and \(|H_H|\) are defined in order to ensure Nominal Performance (NP), Robust Stability (RS) and Robust Performance (RP)</li>
<li>Section <a href="#org4f2901c">5</a>:
\(H_L\) and \(H_H\) are synthesize using the \(\mathcal{H}_\infty\) synthesis such that \(|H_L|\) and \(|H_H|\) are within the specified bounds and such that \(H_L + H_H = 1\) (complementary property)</li>
<li>Section <a href="#org1e02524">6</a>:
Analytical formulas of complementary filters are used</li>
<li>Section <a href="#orgcb66bdb">7</a>:
The obtained complementary filters using the \(\mathcal{H}_\infty\) synthesis and the analytical formula are compared</li>
<li>Section <a href="#org06b9c47">8</a>:
The obtain controller \(K = G^{-1} H_H^{-1}\) is analyzed</li>
<li>Section <a href="#org55c5c37">9</a>:
The Nominal Stability (NS) and Nominal Performance conditions are verified</li>
<li>Section <a href="#org4d2c1c6">10</a>:
Robust Stability and Robust Performance conditions are studied</li>
<li>Section <a href="#org7ea3388">11</a>:
A pre-filter that is used to limit the input usage due to the change of the reference is added</li>
2020-10-08 10:52:20 +02:00
</ul>
2020-11-12 10:44:04 +01:00
<div class="note" id="orgbb67345">
2020-10-08 10:52:20 +02:00
<p>
All the files (data and Matlab scripts) are accessible <a href="data/sensor_fusion.zip">here</a>.
</p>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-org917402b" class="outline-2">
<h2 id="org917402b"><span class="section-number-2">1</span> Definition of the plant</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-1">
<p>
2020-11-12 10:44:04 +01:00
<a id="orgfbfc3f4"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
The studied system consists of a solid positioned on top of a motorized uni-axial soft suspension.
</p>
<p>
The absolute position \(x\) of the solid is measured using an inertial sensor and a force \(F\) can be applied to the mass using a voice coil actuator.
</p>
<p>
2020-11-12 10:44:04 +01:00
The model of the system is represented on figure <a href="#orgf41ac95">2</a> where the mass of the solid is \(m = 10\ [kg]\), the stiffness of the suspension is \(k = 10^4\ [N/m]\) and the damping of the system is \(c = 10^2\ [N/(m/s)]\).
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="orgf41ac95" class="figure">
<p><img src="figs-paper/mech_sys_alone.png" alt="mech_sys_alone.png" />
2020-10-08 10:52:20 +02:00
</p>
<p><span class="figure-number">Figure 2: </span>One degree of freedom system</p>
</div>
<p>
2020-11-12 10:44:04 +01:00
The plant \(G\) is defined on matlab and its bode plot is shown on figure <a href="#org8998aea">3</a>.
2020-10-08 10:52:20 +02:00
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">m = 10; <span class="org-comment">% mass [kg]</span>
k = 1e4; <span class="org-comment">% stiffness [N/m]</span>
c = 1e2; <span class="org-comment">% damping [N/(m/s)]</span>
2020-10-08 10:52:20 +02:00
2020-11-12 10:44:04 +01:00
G = 1<span class="org-type">/</span>(m<span class="org-type">*</span>s<span class="org-type">^</span>2 <span class="org-type">+</span> c<span class="org-type">*</span>s <span class="org-type">+</span> k);
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<div id="org8998aea" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/bode_plot_mech_sys.png" alt="bode_plot_mech_sys.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 3: </span>Bode plot of \(G\)</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-orgde878e2" class="outline-2">
<h2 id="orgde878e2"><span class="section-number-2">2</span> Multiplicative input uncertainty</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-2">
<p>
2020-11-12 10:44:04 +01:00
<a id="orge119055"></a>
2020-10-08 10:52:20 +02:00
We choose to use the multiplicative input uncertainty to model the plant uncertainty:
</p>
2020-11-12 10:44:04 +01:00
\begin{equation}
\Pi_I: \ G^\prime(s) = G(s) (1 + w_I(s) \Delta(s)),\text{ with } |\Delta(j\omega)| < 1 \ \forall \omega
\end{equation}
2020-10-08 10:52:20 +02:00
<p>
2020-11-12 10:44:04 +01:00
The uncertainty weight is shown in figure <a href="#org4aeff4d">4</a>.
2020-10-08 10:52:20 +02:00
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">wI = createWeight(<span class="org-string">'n'</span>, 2, <span class="org-string">'w0'</span>, 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>80, <span class="org-string">'G0'</span>, 0.1, <span class="org-string">'G1'</span>, 10, <span class="org-string">'Gc'</span>, 1);
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<div id="org4aeff4d" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/bode_wi.png" alt="bode_wi.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 4: </span>Bode plot of \(w_I\)</p>
2020-10-08 10:52:20 +02:00
</div>
<p>
2020-11-12 10:44:04 +01:00
Elements in the uncertainty set \(\Pi_I\) are computed and their bode plot is shown on figure <a href="#org51c57f2">5</a>.
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="org51c57f2" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/plant_uncertainty_bode_plot.png" alt="plant_uncertainty_bode_plot.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 5: </span>Some elements in the uncertainty set \(\Pi_I\)</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-orgc908813" class="outline-2">
<h2 id="orgc908813"><span class="section-number-2">3</span> Specifications and performance weights</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-3">
<p>
2020-11-12 10:44:04 +01:00
<a id="org3dfa664"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
The control objective is to isolate the displacement \(x\) of the mass from the ground motion \(w\).
</p>
<p>
The specifications are described below:
</p>
<ul class="org-ul">
<li>at least a factor \(10\) of disturbance rejection at \(2\ \text{Hz}\) and with a slope of \(2\) below \(2\ \text{Hz}\) until a rejection of \(10^3\)</li>
<li>the noise attenuation should be at least \(10\) above \(100\ \text{Hz}\) and with a slope of \(-2\) above</li>
</ul>
<p>
2020-11-12 10:44:04 +01:00
These specifications can be represented as upper bounds on the closed loop transfer functions \(S\) and \(T\) (see figure <a href="#org8dcc955">6</a>).
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="org8dcc955" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/bode_requirements.png" alt="bode_requirements.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 6: </span>Upper bounds on \(S\) and \(T\)</p>
2020-10-08 10:52:20 +02:00
</div>
<p>
We now define two weights, \(w_S(s)\) and \(w_T(s)\) such that \(1/|w_S|\) and \(1/|w_T|\) are lower than the previously defined upper bounds.
Then, the performance specifications are satisfied if the following condition is valid:
\[ \big|S(j\omega)\big| < \frac{1}{|w_S(j\omega)|} ; \quad \big|T(j\omega)\big| < \frac{1}{|w_T(j\omega)|}, \quad \forall \omega \]
</p>
<p>
2020-11-12 10:44:04 +01:00
The weights are defined as follow. They magnitude is compared with the upper bounds on \(S\) and \(T\) on figure <a href="#org702c691">7</a>.
2020-10-08 10:52:20 +02:00
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">wS = 1600<span class="org-type">/</span>(s<span class="org-type">+</span>0.13)<span class="org-type">^</span>2;
wT = 1000<span class="org-type">*</span>((s<span class="org-type">/</span>(2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>1000)))<span class="org-type">^</span>2;
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<div id="org702c691" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/compare_weights_upper_bounds_S_T.png" alt="compare_weights_upper_bounds_S_T.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 7: </span>Weights \(w_S\) and \(w_T\) with the upper bounds on \(S\) and \(T\) obtained from the specifications</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-orgce7235f" class="outline-2">
<h2 id="orgce7235f"><span class="section-number-2">4</span> Upper bounds on the norm of the complementary filters for NP, RS and RP</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-4">
<p>
2020-11-12 10:44:04 +01:00
<a id="orgacb23fe"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
Now that we have defined \(w_I\), \(w_S\) and \(w_T\), we can derive conditions for Nominal Performance, Robust Stability and Robust Performance (\(j\omega\) is omitted here for readability):
</p>
\begin{align*}
\text{NP} &\Leftrightarrow |H_H| < \frac{1}{|w_S|} \text{ and } |H_L| < \frac{1}{|w_T|} \quad \forall \omega \\
\text{RS} &\Leftrightarrow |H_L| < \frac{1}{|w_I| (2 + |w_I|)} \quad \forall \omega \\
\text{RP for } S &\Leftarrow |H_H| < \frac{1 + |w_I|}{|w_S| (2 + |w_I|)} \quad \forall \omega \\
\text{RP for } T &\Leftrightarrow |H_L| < \frac{1}{|w_T| (1 + |w_I|) + |w_I|} \quad \forall \omega
\end{align*}
<p>
These conditions are upper bounds on the complementary filters used for control.
</p>
<p>
2020-11-12 10:44:04 +01:00
We plot these conditions on figure <a href="#org3739909">8</a>.
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="org3739909" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/weights_NP_RS_RP.png" alt="weights_NP_RS_RP.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 8: </span>Upper bounds on the norm of the complementary filters for NP, RS and RP</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-orgaf8f3b1" class="outline-2">
<h2 id="orgaf8f3b1"><span class="section-number-2">5</span> H-Infinity synthesis of complementary filters</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-5">
<p>
2020-11-12 10:44:04 +01:00
<a id="org4f2901c"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
We here synthesize the complementary filters using the \(\mathcal{H}_\infty\) synthesis.
The goal is to specify upper bounds on the norms of \(H_L\) and \(H_H\) while ensuring their complementary property (\(H_L + H_H = 1\)).
</p>
<p>
2020-11-12 10:44:04 +01:00
In order to do so, we use the generalized plant shown on figure <a href="#org3fbed5b">9</a> where \(w_L\) and \(w_H\) weighting transfer functions that will be used to shape \(H_L\) and \(H_H\) respectively.
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="org3fbed5b" class="figure">
<p><img src="figs-paper/sf_hinf_filters_plant_b.png" alt="sf_hinf_filters_plant_b.png" />
2020-10-08 10:52:20 +02:00
</p>
<p><span class="figure-number">Figure 9: </span>Generalized plant used for the \(\mathcal{H}_\infty\) synthesis of the complementary filters</p>
</div>
<p>
2020-11-12 10:44:04 +01:00
The \(\mathcal{H}_\infty\) synthesis applied on this generalized plant will give a transfer function \(H_L\) (figure <a href="#org3939c2e">10</a>) such that the \(\mathcal{H}_\infty\) norm of the transfer function from \(w\) to \([z_H,\ z_L]\) is less than one:
2020-10-08 10:52:20 +02:00
\[ \left\| \begin{array}{c} H_L w_L \\ (1 - H_L) w_H \end{array} \right\|_\infty < 1 \]
</p>
<p>
Thus, if the above condition is verified, we can define \(H_H = 1 - H_L\) and we have that:
\[ \left\| \begin{array}{c} H_L w_L \\ H_H w_H \end{array} \right\|_\infty < 1 \]
Which is almost (with an maximum error of \(\sqrt{2}\)) equivalent to:
</p>
\begin{align*}
|H_L| &< \frac{1}{|w_L|}, \quad \forall \omega \\
|H_H| &< \frac{1}{|w_H|}, \quad \forall \omega
\end{align*}
<p>
We then see that \(w_L\) and \(w_H\) can be used to shape both \(H_L\) and \(H_H\) while ensuring (by definition of \(H_H = 1 - H_L\)) their complementary property.
</p>
2020-11-12 10:44:04 +01:00
<div id="org3939c2e" class="figure">
<p><img src="figs-paper/sf_hinf_filters_b.png" alt="sf_hinf_filters_b.png" />
2020-10-08 10:52:20 +02:00
</p>
<p><span class="figure-number">Figure 10: </span>\(\mathcal{H}_\infty\) synthesis of the complementary filters</p>
</div>
<p>
2020-11-12 10:44:04 +01:00
Thus, if we choose \(w_L\) and \(w_H\) such that \(1/|w_L|\) and \(1/|w_H|\) lie below the upper bounds of figure <a href="#org3739909">8</a>, we will ensure the NP, RS and RP of the controlled system.
2020-10-08 10:52:20 +02:00
</p>
<p>
Depending if we are interested only in NP, RS or RP, we can adjust the weights \(w_L\) and \(w_H\).
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">omegab = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>9;
wH = (omegab)<span class="org-type">^</span>2<span class="org-type">/</span>(s <span class="org-type">+</span> omegab<span class="org-type">*</span>sqrt(1e<span class="org-type">-</span>5))<span class="org-type">^</span>2;
omegab = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>28;
wL = (s <span class="org-type">+</span> omegab<span class="org-type">/</span>(4.5)<span class="org-type">^</span>(1<span class="org-type">/</span>3))<span class="org-type">^</span>3<span class="org-type">/</span>(s<span class="org-type">*</span>(1e<span class="org-type">-</span>4)<span class="org-type">^</span>(1<span class="org-type">/</span>3) <span class="org-type">+</span> omegab)<span class="org-type">^</span>3;
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<div id="org90a3139" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/weights_wl_wh.png" alt="weights_wl_wh.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 11: </span>Weights on the complementary filters \(w_L\) and \(w_H\) and the associated performance weights</p>
2020-10-08 10:52:20 +02:00
</div>
<p>
We define the generalized plant \(P\) on matlab.
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">P = [0 wL;
2020-10-08 10:52:20 +02:00
wH <span class="org-type">-</span>wH;
2020-11-12 10:44:04 +01:00
1 0];
2020-10-08 10:52:20 +02:00
</pre>
</div>
<p>
And we do the \(\mathcal{H}_\infty\) synthesis using the <code>hinfsyn</code> command.
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">[Hl_hinf, <span class="org-type">~</span>, gamma, <span class="org-type">~</span>] = hinfsyn(P, 1, 1,<span class="org-string">'TOLGAM'</span>, 0.001, <span class="org-string">'METHOD'</span>, <span class="org-string">'ric'</span>, <span class="org-string">'DISPLAY'</span>, <span class="org-string">'on'</span>);
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<pre class="example" id="org67390b9">
2020-10-08 10:52:20 +02:00
[Hl_hinf, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
Test bounds: 0.0000 &lt; gamma &lt;= 1.7285
gamma hamx_eig xinf_eig hamy_eig yinf_eig nrho_xy p/f
1.729 4.1e+01 8.4e-12 1.8e-01 0.0e+00 0.0000 p
0.864 3.9e+01 -5.8e-02# 1.8e-01 0.0e+00 0.0000 f
1.296 4.0e+01 8.4e-12 1.8e-01 0.0e+00 0.0000 p
1.080 4.0e+01 8.5e-12 1.8e-01 0.0e+00 0.0000 p
0.972 3.9e+01 -4.2e-01# 1.8e-01 0.0e+00 0.0000 f
1.026 4.0e+01 8.5e-12 1.8e-01 0.0e+00 0.0000 p
0.999 3.9e+01 8.5e-12 1.8e-01 0.0e+00 0.0000 p
0.986 3.9e+01 -1.2e+00# 1.8e-01 0.0e+00 0.0000 f
0.993 3.9e+01 -8.2e+00# 1.8e-01 0.0e+00 0.0000 f
0.996 3.9e+01 8.5e-12 1.8e-01 0.0e+00 0.0000 p
0.994 3.9e+01 8.5e-12 1.8e-01 0.0e+00 0.0000 p
0.993 3.9e+01 -3.2e+01# 1.8e-01 0.0e+00 0.0000 f
Gamma value achieved: 0.9942
</pre>
<p>
2020-11-12 10:44:04 +01:00
We then define the high pass filter \(H_H = 1 - H_L\). The bode plot of both \(H_L\) and \(H_H\) is shown on figure <a href="#orgebb3cd3">12</a>.
2020-10-08 10:52:20 +02:00
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">Hh_hinf = 1 <span class="org-type">-</span> Hl_hinf;
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<div id="orgebb3cd3" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/hinf_filters_results.png" alt="hinf_filters_results.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 12: </span>Obtained complementary filters using \(\mathcal{H}_\infty\) synthesis</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-orga721ae7" class="outline-2">
<h2 id="orga721ae7"><span class="section-number-2">6</span> Complementary filters using analytical formula</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-6">
<p>
2020-11-12 10:44:04 +01:00
<a id="org1e02524"></a>
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<ul class="org-ul">
<li class="off"><code>[&#xa0;]</code> The problem here is that the High pass filter gain goes to zero at high frequency, and when inverted, it will go to infinity.</li>
</ul>
2020-10-08 10:52:20 +02:00
<p>
We here use analytical formula for the complementary filters \(H_L\) and \(H_H\).
</p>
<p>
The first two formulas that are used to generate complementary filters are:
</p>
\begin{align*}
H_L(s) &= \frac{(1+\alpha) (\frac{s}{\omega_0})+1}{\left((\frac{s}{\omega_0})+1\right) \left((\frac{s}{\omega_0})^2 + \alpha (\frac{s}{\omega_0}) + 1\right)}\\
H_H(s) &= \frac{(\frac{s}{\omega_0})^2 \left((\frac{s}{\omega_0})+1+\alpha\right)}{\left((\frac{s}{\omega_0})+1\right) \left((\frac{s}{\omega_0})^2 + \alpha (\frac{s}{\omega_0}) + 1\right)}
\end{align*}
<p>
where:
</p>
<ul class="org-ul">
<li>\(\omega_0\) is the blending frequency in rad/s.</li>
<li>\(\alpha\) is used to change the shape of the filters:
<ul class="org-ul">
<li>Small values for \(\alpha\) will produce high magnitude of the filters \(|H_L(j\omega)|\) and \(|H_H(j\omega)|\) near \(\omega_0\) but smaller value for \(|H_L(j\omega)|\) above \(\approx 1.5 \omega_0\) and for \(|H_H(j\omega)|\) below \(\approx 0.7 \omega_0\)</li>
<li>A large \(\alpha\) will do the opposite</li>
</ul></li>
</ul>
<p>
2020-11-12 10:44:04 +01:00
This is illustrated on figure <a href="#org72b4c79">13</a>.
2020-10-08 10:52:20 +02:00
As it is usually wanted to have the \(\| S \|_\infty < 2\), \(\alpha\) between \(0.5\) and \(1\) gives a good trade-off between the performance and the robustness.
The slope of those filters at high and low frequencies is \(-2\) and \(2\) respectively for \(H_L\) and \(H_H\).
</p>
2020-11-12 10:44:04 +01:00
<div id="org72b4c79" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/comp_filters_param_alpha.png" alt="comp_filters_param_alpha.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 13: </span>Effect of the parameter \(\alpha\) on the shape of the generated second order complementary filters</p>
2020-10-08 10:52:20 +02:00
</div>
<p>
The parameters \(\alpha\) and \(\omega_0\) are chosen in order to have that the complementary filters stay below the defined upper bounds.
</p>
<p>
2020-11-12 10:44:04 +01:00
The obtained complementary filters are shown on figure <a href="#org8e27bc0">14</a>.
2020-10-08 10:52:20 +02:00
The Robust Performance is not fulfilled for \(T\), and we see that the RP condition as a slop of \(-3\). We thus have to use different formula for the complementary filters here.
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">w0 = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>13;
alpha = 0.8;
2020-10-08 10:52:20 +02:00
2020-11-12 10:44:04 +01:00
Hh2_ana = (s<span class="org-type">/</span>w0)<span class="org-type">^</span>2<span class="org-type">*</span>((s<span class="org-type">/</span>w0)<span class="org-type">+</span>1<span class="org-type">+</span>alpha)<span class="org-type">/</span>(((s<span class="org-type">/</span>w0)<span class="org-type">+</span>1)<span class="org-type">*</span>((s<span class="org-type">/</span>w0)<span class="org-type">^</span>2 <span class="org-type">+</span> alpha<span class="org-type">*</span>(s<span class="org-type">/</span>w0) <span class="org-type">+</span> 1));
Hl2_ana = ((1<span class="org-type">+</span>alpha)<span class="org-type">*</span>(s<span class="org-type">/</span>w0)<span class="org-type">+</span>1)<span class="org-type">/</span>(((s<span class="org-type">/</span>w0)<span class="org-type">+</span>1)<span class="org-type">*</span>((s<span class="org-type">/</span>w0)<span class="org-type">^</span>2 <span class="org-type">+</span> alpha<span class="org-type">*</span>(s<span class="org-type">/</span>w0) <span class="org-type">+</span> 1));
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<div id="org8e27bc0" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/complementary_filters_second_order.png" alt="complementary_filters_second_order.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 14: </span>Second order complementary filters using the analytical formula</p>
2020-10-08 10:52:20 +02:00
</div>
<p>
The following formula gives complementary filters with slopes of \(-3\) and \(3\):
</p>
\begin{align*}
H_L(s) &= \frac{\left(1+(\alpha+1)(\beta+1)\right) (\frac{s}{\omega_0})^2 + (1+\alpha+\beta)(\frac{s}{\omega_0}) + 1}{\left(\frac{s}{\omega_0} + 1\right) \left( (\frac{s}{\omega_0})^2 + \alpha (\frac{s}{\omega_0}) + 1 \right) \left( (\frac{s}{\omega_0})^2 + \beta (\frac{s}{\omega_0}) + 1 \right)}\\
H_H(s) &= \frac{(\frac{s}{\omega_0})^3 \left( (\frac{s}{\omega_0})^2 + (1+\alpha+\beta) (\frac{s}{\omega_0}) + (1+(\alpha+1)(\beta+1)) \right)}{\left(\frac{s}{\omega_0} + 1\right) \left( (\frac{s}{\omega_0})^2 + \alpha (\frac{s}{\omega_0}) + 1 \right) \left( (\frac{s}{\omega_0})^2 + \beta (\frac{s}{\omega_0}) + 1 \right)}
\end{align*}
<p>
The parameters are:
</p>
<ul class="org-ul">
<li>\(\omega_0\) is the blending frequency in rad/s</li>
<li>\(\alpha\) and \(\beta\) that are used to change the shape of the filters similarly to the parameter \(\alpha\) for the second order complementary filters</li>
</ul>
<p>
2020-11-12 10:44:04 +01:00
The filters are defined below and the result is shown on figure <a href="#org9e9b7f8">15</a> where we can see that the complementary filters are below the defined upper bounds.
2020-10-08 10:52:20 +02:00
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">alpha = 1;
beta = 10;
w0 = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>14;
2020-10-08 10:52:20 +02:00
2020-11-12 10:44:04 +01:00
Hh3_ana = (s<span class="org-type">/</span>w0)<span class="org-type">^</span>3 <span class="org-type">*</span> ((s<span class="org-type">/</span>w0)<span class="org-type">^</span>2 <span class="org-type">+</span> (1<span class="org-type">+</span>alpha<span class="org-type">+</span>beta)<span class="org-type">*</span>(s<span class="org-type">/</span>w0) <span class="org-type">+</span> (1<span class="org-type">+</span>(alpha<span class="org-type">+</span>1)<span class="org-type">*</span>(beta<span class="org-type">+</span>1)))<span class="org-type">/</span>((s<span class="org-type">/</span>w0 <span class="org-type">+</span> 1)<span class="org-type">*</span>((s<span class="org-type">/</span>w0)<span class="org-type">^</span>2<span class="org-type">+</span>alpha<span class="org-type">*</span>(s<span class="org-type">/</span>w0)<span class="org-type">+</span>1)<span class="org-type">*</span>((s<span class="org-type">/</span>w0)<span class="org-type">^</span>2<span class="org-type">+</span>beta<span class="org-type">*</span>(s<span class="org-type">/</span>w0)<span class="org-type">+</span>1));
Hl3_ana = ((1<span class="org-type">+</span>(alpha<span class="org-type">+</span>1)<span class="org-type">*</span>(beta<span class="org-type">+</span>1))<span class="org-type">*</span>(s<span class="org-type">/</span>w0)<span class="org-type">^</span>2 <span class="org-type">+</span> (1<span class="org-type">+</span>alpha<span class="org-type">+</span>beta)<span class="org-type">*</span>(s<span class="org-type">/</span>w0) <span class="org-type">+</span> 1)<span class="org-type">/</span>((s<span class="org-type">/</span>w0 <span class="org-type">+</span> 1)<span class="org-type">*</span>((s<span class="org-type">/</span>w0)<span class="org-type">^</span>2<span class="org-type">+</span>alpha<span class="org-type">*</span>(s<span class="org-type">/</span>w0)<span class="org-type">+</span>1)<span class="org-type">*</span>((s<span class="org-type">/</span>w0)<span class="org-type">^</span>2<span class="org-type">+</span>beta<span class="org-type">*</span>(s<span class="org-type">/</span>w0)<span class="org-type">+</span>1));
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<div id="org9e9b7f8" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/complementary_filters_third_order.png" alt="complementary_filters_third_order.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 15: </span>Third order complementary filters using the analytical formula</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-orgaf265db" class="outline-2">
<h2 id="orgaf265db"><span class="section-number-2">7</span> Comparison of complementary filters</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-7">
<p>
2020-11-12 10:44:04 +01:00
<a id="orgcb66bdb"></a>
The generated complementary filters using \(\mathcal{H}_\infty\) and the analytical formulas are compared on figure <a href="#org0e60c80">16</a>.
2020-10-08 10:52:20 +02:00
</p>
<p>
Although they are very close to each other, there is some difference to note here:
</p>
<ul class="org-ul">
<li>the analytical formula provides a very simple way to generate the complementary filters (and thus the controller), they could even be used to tune the controller online using the parameters \(\alpha\) and \(\omega_0\). However, these formula have the property that \(|H_H|\) and \(|H_L|\) are symmetrical with the frequency \(\omega_0\) which may not be desirable.</li>
<li>while the \(\mathcal{H}_\infty\) synthesis of the complementary filters is not as straightforward as using the analytical formula, it provides a more optimized procedure to obtain the complementary filters</li>
</ul>
<p>
The complementary filters obtained with the \(\mathcal{H}_\infty\) will be used for further analysis.
</p>
2020-11-12 10:44:04 +01:00
<div id="org0e60c80" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/comp_hinf_analytical.png" alt="comp_hinf_analytical.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 16: </span>Comparison of the complementary filters obtained with \(\mathcal{H}_\infty\) synthesis and with the analytical formula</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-orga1bd48e" class="outline-2">
<h2 id="orga1bd48e"><span class="section-number-2">8</span> Controller Analysis</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-8">
<p>
2020-11-12 10:44:04 +01:00
<a id="org06b9c47"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
The controller \(K\) is computed from the plant model \(G\) and the low pass filter \(H_H\):
\[ K = G^{-1} H_H^{-1} \]
</p>
<p>
As this is not proper and thus realizable, a second order low pass filter is added with a crossover frequency much larger than the control bandwidth.
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">omega = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>500;
K = 1<span class="org-type">/</span>(Hh_hinf<span class="org-type">*</span>G) <span class="org-type">*</span> 1<span class="org-type">/</span>((1<span class="org-type">+</span>s<span class="org-type">/</span>omega)<span class="org-type">*</span>(1<span class="org-type">+</span>s<span class="org-type">/</span>omega<span class="org-type">+</span>(s<span class="org-type">/</span>omega)<span class="org-type">^</span>2));
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<pre class="example" id="org4380854">
2020-10-08 10:52:20 +02:00
zpk(K)
ans =
4.961e12 (s+9.915e04) (s^2 + 5s + 500) (s^2 + 284.6s + 2.135e04) (s^2 + 130.5s + 9887)
--------------------------------------------------------------------------------------------------
(s+9.914e04) (s+6283) (s^2 + 0.3576s + 0.03198) (s^2 + 413.8s + 6.398e04) (s^2 + 6283s + 3.948e07)
Continuous-time zero/pole/gain model.
</pre>
<p>
2020-11-12 10:44:04 +01:00
The bode plot of the controller is shown on figure <a href="#orga448a8a">17</a>:
2020-10-08 10:52:20 +02:00
</p>
<ul class="org-ul">
<li>two integrator are present at low frequency</li>
<li>the resonance of the plant at \(3.5\ \text{Hz}\) is inverted (notched)</li>
<li>a lead is added at \(10\ \text{Hz}\)</li>
</ul>
2020-11-12 10:44:04 +01:00
<div id="orga448a8a" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/bode_plot_controller.png" alt="bode_plot_controller.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 17: </span>Bode plot of the obtained controller \(K\)</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-org1eb8c07" class="outline-2">
<h2 id="org1eb8c07"><span class="section-number-2">9</span> Nominal Stability and Nominal Performance</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-9">
<p>
2020-11-12 10:44:04 +01:00
<a id="org55c5c37"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
The nominal stability of the system is first checked with the <code>allmargin</code> matlab command.
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">allmargin(K<span class="org-type">*</span>G<span class="org-type">*</span>Hl_hinf)
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<pre class="example" id="orgab687e2">
2020-10-08 10:52:20 +02:00
allmargin(K*G*Hl_hinf)
ans =
struct with fields:
2020-11-12 10:44:04 +01:00
GainMargin: 3.86196691340257
GMFrequency: 225.992831144248
PhaseMargin: 34.0890006269983
PMFrequency: 88.3632737129296
DelayMargin: 0.0067331740287082
DMFrequency: 88.3632737129296
2020-10-08 10:52:20 +02:00
Stable: 1
</pre>
<p>
The system is stable and the stability margins are good.
</p>
<p>
2020-11-12 10:44:04 +01:00
The bode plot of the loop gain \(L = K*G*H_L\) is shown on figure <a href="#orgcce877a">18</a>.
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="orgcce877a" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/bode_plot_loop_gain.png" alt="bode_plot_loop_gain.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 18: </span>Bode Plot of the Loop Gain \(L = K G H_L\)</p>
2020-10-08 10:52:20 +02:00
</div>
<p>
In order to check the Nominal Performance of the system, we compute the sensibility and the complementary sensibility transfer functions.
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">S = 1<span class="org-type">/</span>(K<span class="org-type">*</span>G<span class="org-type">*</span>Hl_hinf <span class="org-type">+</span> 1);
T = K<span class="org-type">*</span>G<span class="org-type">*</span>Hl_hinf<span class="org-type">/</span>(K<span class="org-type">*</span>G<span class="org-type">*</span>Hl_hinf <span class="org-type">+</span> 1);
2020-10-08 10:52:20 +02:00
</pre>
</div>
<p>
2020-11-12 10:44:04 +01:00
We then compare their norms with the upper bounds on the performance of the system (figure <a href="#orgbbc226a">19</a>).
2020-10-08 10:52:20 +02:00
As expected, we guarantee the Nominal Performance of the system.
</p>
2020-11-12 10:44:04 +01:00
<div id="orgbbc226a" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/verification_NP.png" alt="verification_NP.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 19: </span>Bode plot of \(S\) and \(T\) in order to verify the nominal performance of the system</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-org8c33736" class="outline-2">
<h2 id="org8c33736"><span class="section-number-2">10</span> Robust Stability and Robust Performance</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-10">
<p>
2020-11-12 10:44:04 +01:00
<a id="org4d2c1c6"></a>
2020-10-08 10:52:20 +02:00
In order to verify the Robust stability of the system, we can use the following equivalence:
\[ \text{RS} \Leftrightarrow \left| w_I T \right| < 1 \quad \forall \omega \]
</p>
<p>
2020-11-12 10:44:04 +01:00
This is shown on figure <a href="#orga0177a7">20</a>.
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="orga0177a7" class="figure">
<p><img src="figs/robust_stability.png" alt="robust_stability.png" />
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 20: </span>Robust Stability Check: \(|w_I T| < 1, \quad \forall \omega\)</p>
2020-10-08 10:52:20 +02:00
</div>
<p>
2020-11-12 10:44:04 +01:00
To check Robust Stability, we can also look at the loop gain of the uncertain system (figure <a href="#orgfb25f29">21</a>) or the Nyquist plot (figure <a href="#orgf81e346">22</a>).
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="orgfb25f29" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/loop_gain_robustness.png" alt="loop_gain_robustness.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 21: </span>Loop Gain of the uncertain system</p>
2020-10-08 10:52:20 +02:00
</div>
2020-11-12 10:44:04 +01:00
<div id="orgf81e346" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/nyquist_robustness.png" alt="nyquist_robustness.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 22: </span>Nyquist plot of the uncertain system</p>
2020-10-08 10:52:20 +02:00
</div>
<p>
The Robust Performance is verified by plotting \(|S|\) and \(|T|\) for the uncertain system along side the upper bounds defined for performance.
2020-11-12 10:44:04 +01:00
This is shown on figure <a href="#org58bb0a6">23</a> and we can indeed confirmed that the robust performance property of the system is valid.
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div id="org58bb0a6" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/robust_performance_result.png" alt="robust_performance_result.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 23: </span>Verification of the Robust Performance</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-org367d76d" class="outline-2">
<h2 id="org367d76d"><span class="section-number-2">11</span> Pre-filter</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-11">
<p>
2020-11-12 10:44:04 +01:00
<a id="org7ea3388"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
For now, we have not specified any performance requirements on the input usage due to change of reference.
2020-11-12 10:44:04 +01:00
Do limit the input usage due to change of reference, we can use a pre-filter \(K_r\) as shown on figure <a href="#org16aa33e">1</a>.
2020-10-08 10:52:20 +02:00
</p>
<p>
If we want a response time of 0.01 second, we can choose a first order low pass filter with a crossover frequency at \(1/0.01 = 100\ \text{Hz}\) as defined below.
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">Kr = 1<span class="org-type">/</span>(1<span class="org-type">+</span>s<span class="org-type">/</span>2<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">/</span>100);
2020-10-08 10:52:20 +02:00
</pre>
</div>
<p>
2020-11-12 10:44:04 +01:00
We then run a simulation for a step of \(1\mu m\) for the system without and with the pre-filter \(K_r\) (figure <a href="#orgb70a4d3">24</a>).
2020-10-08 10:52:20 +02:00
This confirms that a pre-filter can be used to limit the input usage due to change of reference using this architecture.
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab">t = linspace(0, 0.01, 500);
2020-10-08 10:52:20 +02:00
opts = stepDataOptions;
2020-11-12 10:44:04 +01:00
opts.StepAmplitude = 1e<span class="org-type">-</span>6;
2020-10-08 10:52:20 +02:00
2020-11-12 10:44:04 +01:00
u = step((K)<span class="org-type">/</span>(1<span class="org-type">+</span>G<span class="org-type">*</span>K<span class="org-type">*</span>Hl_hinf), t, opts);
uf = step((Kr<span class="org-type">*</span>K)<span class="org-type">/</span>(1<span class="org-type">+</span>G<span class="org-type">*</span>K<span class="org-type">*</span>Hl_hinf), t, opts);
y = step((K<span class="org-type">*</span>G)<span class="org-type">/</span>(1<span class="org-type">+</span>G<span class="org-type">*</span>K<span class="org-type">*</span>Hl_hinf), t, opts);
yf = step((Kr<span class="org-type">*</span>G<span class="org-type">*</span>K)<span class="org-type">/</span>(1<span class="org-type">+</span>G<span class="org-type">*</span>K<span class="org-type">*</span>Hl_hinf), t, opts);
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
<div id="orgb70a4d3" class="figure">
2020-10-08 10:52:20 +02:00
<p><img src="figs/u_and_y_with_Kr.png" alt="u_and_y_with_Kr.png" />
</p>
2020-11-12 10:44:04 +01:00
<p><span class="figure-number">Figure 24: </span>Input usage and response due to a step change of reference when using a pre-filter \(K_r\)</p>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-org846e5c4" class="outline-2">
<h2 id="org846e5c4"><span class="section-number-2">12</span> Matlab Functions</h2>
2020-10-08 10:52:20 +02:00
<div class="outline-text-2" id="text-12">
<p>
2020-11-12 10:44:04 +01:00
<a id="org4fb0718"></a>
</p>
</div>
<div id="outline-container-org2f50732" class="outline-3">
<h3 id="org2f50732"><span class="section-number-3">12.1</span> <code>createWeight</code></h3>
<div class="outline-text-3" id="text-12-1">
<p>
<a id="org01d0043"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
2020-11-12 10:44:04 +01:00
This Matlab function is accessible <a href="src/createWeight.m">here</a>.
2020-10-08 10:52:20 +02:00
</p>
<div class="org-src-container">
2020-11-12 10:44:04 +01:00
<pre class="src src-matlab"><span class="org-keyword">function</span> <span class="org-variable-name">[W]</span> = <span class="org-function-name">createWeight</span>(<span class="org-variable-name">args</span>)
<span class="org-comment">% createWeight -</span>
<span class="org-comment">%</span>
<span class="org-comment">% Syntax: [in_data] = createWeight(in_data)</span>
<span class="org-comment">%</span>
<span class="org-comment">% Inputs:</span>
<span class="org-comment">% - n - Weight Order</span>
<span class="org-comment">% - G0 - Low frequency Gain</span>
<span class="org-comment">% - G1 - High frequency Gain</span>
<span class="org-comment">% - Gc - Gain of W at frequency w0</span>
<span class="org-comment">% - w0 - Frequency at which |W(j w0)| = Gc</span>
<span class="org-comment">%</span>
<span class="org-comment">% Outputs:</span>
<span class="org-comment">% - W - Generated Weight</span>
arguments
args.n (1,1) double {mustBeInteger, mustBePositive} = 1
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.w0 (1,1) double {mustBeNumeric, mustBePositive} = 1
<span class="org-keyword">end</span>
mustBeBetween(args.G0, args.Gc, args.G1);
s = tf(<span class="org-string">'s'</span>);
W = (((1<span class="org-type">/</span>args.w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(args.G0<span class="org-type">/</span>args.Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>args.n))<span class="org-type">/</span>(1<span class="org-type">-</span>(args.Gc<span class="org-type">/</span>args.G1)<span class="org-type">^</span>(2<span class="org-type">/</span>args.n)))<span class="org-type">*</span>s <span class="org-type">+</span> (args.G0<span class="org-type">/</span>args.Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>args.n))<span class="org-type">/</span>((1<span class="org-type">/</span>args.G1)<span class="org-type">^</span>(1<span class="org-type">/</span>args.n)<span class="org-type">*</span>(1<span class="org-type">/</span>args.w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(args.G0<span class="org-type">/</span>args.Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>args.n))<span class="org-type">/</span>(1<span class="org-type">-</span>(args.Gc<span class="org-type">/</span>args.G1)<span class="org-type">^</span>(2<span class="org-type">/</span>args.n)))<span class="org-type">*</span>s <span class="org-type">+</span> (1<span class="org-type">/</span>args.Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>args.n)))<span class="org-type">^</span>args.n;
<span class="org-keyword">end</span>
<span class="org-comment">% Custom validation function</span>
<span class="org-keyword">function</span> <span class="org-function-name">mustBeBetween</span>(<span class="org-variable-name">a</span>,<span class="org-variable-name">b</span>,<span class="org-variable-name">c</span>)
<span class="org-keyword">if</span> <span class="org-type">~</span>((a <span class="org-type">&gt;</span> b <span class="org-type">&amp;&amp;</span> b <span class="org-type">&gt;</span> c) <span class="org-type">||</span> (c <span class="org-type">&gt;</span> b <span class="org-type">&amp;&amp;</span> b <span class="org-type">&gt;</span> a))
eid = <span class="org-string">'createWeight:inputError'</span>;
msg = <span class="org-string">'Gc should be between G0 and G1.'</span>;
throwAsCaller(MException(eid,msg))
<span class="org-keyword">end</span>
<span class="org-keyword">end</span>
2020-10-08 10:52:20 +02:00
</pre>
</div>
2020-11-12 10:44:04 +01:00
</div>
</div>
2020-10-08 10:52:20 +02:00
2020-11-12 10:44:04 +01:00
<div id="outline-container-org8ef5bfb" class="outline-3">
<h3 id="org8ef5bfb"><span class="section-number-3">12.2</span> <code>plotMagUncertainty</code></h3>
<div class="outline-text-3" id="text-12-2">
2020-10-08 10:52:20 +02:00
<p>
2020-11-12 10:44:04 +01:00
<a id="orgc93a4bb"></a>
2020-10-08 10:52:20 +02:00
</p>
<p>
2020-11-12 10:44:04 +01:00
This Matlab function is accessible <a href="src/plotMagUncertainty.m">here</a>.
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-keyword">function</span> <span class="org-variable-name">[p]</span> = <span class="org-function-name">plotMagUncertainty</span>(<span class="org-variable-name">W</span>, <span class="org-variable-name">freqs</span>, <span class="org-variable-name">args</span>)
<span class="org-comment">% plotMagUncertainty -</span>
<span class="org-comment">%</span>
<span class="org-comment">% Syntax: [p] = plotMagUncertainty(W, freqs, args)</span>
<span class="org-comment">%</span>
<span class="org-comment">% Inputs:</span>
<span class="org-comment">% - W - Multiplicative Uncertainty Weight</span>
<span class="org-comment">% - freqs - Frequency Vector [Hz]</span>
<span class="org-comment">% - args - Optional Arguments:</span>
<span class="org-comment">% - G</span>
<span class="org-comment">% - color_i</span>
<span class="org-comment">% - opacity</span>
<span class="org-comment">%</span>
<span class="org-comment">% Outputs:</span>
<span class="org-comment">% - p - Plot Handle</span>
arguments
W
freqs double {mustBeNumeric, mustBeNonnegative}
args.G = tf(1)
args.color_i (1,1) double {mustBeInteger, mustBeNonnegative} = 0
args.opacity (1,1) double {mustBeNumeric, mustBeNonnegative} = 0.3
args.DisplayName char = <span class="org-string">''</span>
<span class="org-keyword">end</span>
<span class="org-comment">% Get defaults colors</span>
colors = <span class="org-type">get</span>(<span class="org-variable-name">groot</span>, <span class="org-string">'defaultAxesColorOrder'</span>);
p = <span class="org-type">patch</span>([freqs flip(freqs)], ...
[abs(squeeze(freqresp(args.G, freqs, <span class="org-string">'Hz'</span>)))<span class="org-type">.*</span>(1 <span class="org-type">+</span> abs(squeeze(freqresp(W, freqs, <span class="org-string">'Hz'</span>)))); ...
flip(abs(squeeze(freqresp(args.G, freqs, <span class="org-string">'Hz'</span>)))<span class="org-type">.*</span>max(1 <span class="org-type">-</span> abs(squeeze(freqresp(W, freqs, <span class="org-string">'Hz'</span>))), 1e<span class="org-type">-</span>6))], <span class="org-string">'w'</span>, ...
<span class="org-string">'DisplayName'</span>, args.DisplayName);
<span class="org-keyword">if</span> args.color_i <span class="org-type">==</span> 0
p.FaceColor = [0; 0; 0];
<span class="org-keyword">else</span>
p.FaceColor = colors(args.color_i, <span class="org-type">:</span>);
<span class="org-keyword">end</span>
p.EdgeColor = <span class="org-string">'none'</span>;
p.FaceAlpha = args.opacity;
<span class="org-keyword">end</span>
</pre>
</div>
</div>
2020-10-08 10:52:20 +02:00
</div>
2020-11-12 10:44:04 +01:00
<div id="outline-container-org3b855e9" class="outline-3">
<h3 id="org3b855e9"><span class="section-number-3">12.3</span> <code>plotPhaseUncertainty</code></h3>
<div class="outline-text-3" id="text-12-3">
2020-10-08 10:52:20 +02:00
<p>
2020-11-12 10:44:04 +01:00
<a id="org9d5c371"></a>
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<p>
This Matlab function is accessible <a href="src/plotPhaseUncertainty.m">here</a>.
2020-10-08 10:52:20 +02:00
</p>
2020-11-12 10:44:04 +01:00
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-keyword">function</span> <span class="org-variable-name">[p]</span> = <span class="org-function-name">plotPhaseUncertainty</span>(<span class="org-variable-name">W</span>, <span class="org-variable-name">freqs</span>, <span class="org-variable-name">args</span>)
<span class="org-comment">% plotPhaseUncertainty -</span>
<span class="org-comment">%</span>
<span class="org-comment">% Syntax: [p] = plotPhaseUncertainty(W, freqs, args)</span>
<span class="org-comment">%</span>
<span class="org-comment">% Inputs:</span>
<span class="org-comment">% - W - Multiplicative Uncertainty Weight</span>
<span class="org-comment">% - freqs - Frequency Vector [Hz]</span>
<span class="org-comment">% - args - Optional Arguments:</span>
<span class="org-comment">% - G</span>
<span class="org-comment">% - color_i</span>
<span class="org-comment">% - opacity</span>
<span class="org-comment">%</span>
<span class="org-comment">% Outputs:</span>
<span class="org-comment">% - p - Plot Handle</span>
arguments
W
freqs double {mustBeNumeric, mustBeNonnegative}
args.G = tf(1)
args.color_i (1,1) double {mustBeInteger, mustBeNonnegative} = 0
args.opacity (1,1) double {mustBeNumeric, mustBePositive} = 0.3
args.DisplayName char = <span class="org-string">''</span>
<span class="org-keyword">end</span>
<span class="org-comment">% Get defaults colors</span>
colors = <span class="org-type">get</span>(<span class="org-variable-name">groot</span>, <span class="org-string">'defaultAxesColorOrder'</span>);
<span class="org-comment">% Compute Phase Uncertainty</span>
Dphi = 180<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">*</span>asin(abs(squeeze(freqresp(W, freqs, <span class="org-string">'Hz'</span>))));
Dphi(abs(squeeze(freqresp(W, freqs, <span class="org-string">'Hz'</span>))) <span class="org-type">&gt;</span> 1) = 360;
<span class="org-comment">% Compute Plant Phase</span>
G_ang = 180<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">*</span>angle(squeeze(freqresp(args.G, freqs, <span class="org-string">'Hz'</span>)));
p = <span class="org-type">patch</span>([freqs flip(freqs)], [G_ang<span class="org-type">+</span>Dphi; flip(G_ang<span class="org-type">-</span>Dphi)], <span class="org-string">'w'</span>, ...
<span class="org-string">'DisplayName'</span>, args.DisplayName);
<span class="org-keyword">if</span> args.color_i <span class="org-type">==</span> 0
p.FaceColor = [0; 0; 0];
<span class="org-keyword">else</span>
p.FaceColor = colors(args.color_i, <span class="org-type">:</span>);
<span class="org-keyword">end</span>
p.EdgeColor = <span class="org-string">'none'</span>;
p.FaceAlpha = args.opacity;
<span class="org-keyword">end</span>
</pre>
</div>
2020-10-08 10:52:20 +02:00
</div>
</div>
</div>
</div>
<div id="postamble" class="status">
<p class="author">Author: Thomas Dehaeze</p>
2020-11-12 10:44:04 +01:00
<p class="date">Created: 2020-11-12 jeu. 10:43</p>
2020-10-08 10:52:20 +02:00
</div>
</body>
</html>