dehaeze20_virtu_senso_fusio/matlab/index.html

968 lines
46 KiB
HTML

<?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 jeu. 10:43 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<title>Sensor Fusion Paper - Matlab Computation</title>
<meta name="generator" content="Org mode" />
<meta name="author" content="Thomas Dehaeze" />
<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>
</head>
<body>
<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>
<div id="table-of-contents">
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<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>
</ul>
</div>
</div>
<p>
The control architecture studied here is shown on figure <a href="#org16aa33e">1</a> where:
</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>
<div id="org16aa33e" class="figure">
<p><img src="figs-paper/sf_arch_class_prefilter.png" alt="sf_arch_class_prefilter.png" />
</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">
<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>
</ul>
<div class="note" id="orgbb67345">
<p>
All the files (data and Matlab scripts) are accessible <a href="data/sensor_fusion.zip">here</a>.
</p>
</div>
<div id="outline-container-org917402b" class="outline-2">
<h2 id="org917402b"><span class="section-number-2">1</span> Definition of the plant</h2>
<div class="outline-text-2" id="text-1">
<p>
<a id="orgfbfc3f4"></a>
</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>
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)]\).
</p>
<div id="orgf41ac95" class="figure">
<p><img src="figs-paper/mech_sys_alone.png" alt="mech_sys_alone.png" />
</p>
<p><span class="figure-number">Figure 2: </span>One degree of freedom system</p>
</div>
<p>
The plant \(G\) is defined on matlab and its bode plot is shown on figure <a href="#org8998aea">3</a>.
</p>
<div class="org-src-container">
<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>
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);
</pre>
</div>
<div id="org8998aea" class="figure">
<p><img src="figs/bode_plot_mech_sys.png" alt="bode_plot_mech_sys.png" />
</p>
<p><span class="figure-number">Figure 3: </span>Bode plot of \(G\)</p>
</div>
</div>
</div>
<div id="outline-container-orgde878e2" class="outline-2">
<h2 id="orgde878e2"><span class="section-number-2">2</span> Multiplicative input uncertainty</h2>
<div class="outline-text-2" id="text-2">
<p>
<a id="orge119055"></a>
We choose to use the multiplicative input uncertainty to model the plant uncertainty:
</p>
\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}
<p>
The uncertainty weight is shown in figure <a href="#org4aeff4d">4</a>.
</p>
<div class="org-src-container">
<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);
</pre>
</div>
<div id="org4aeff4d" class="figure">
<p><img src="figs/bode_wi.png" alt="bode_wi.png" />
</p>
<p><span class="figure-number">Figure 4: </span>Bode plot of \(w_I\)</p>
</div>
<p>
Elements in the uncertainty set \(\Pi_I\) are computed and their bode plot is shown on figure <a href="#org51c57f2">5</a>.
</p>
<div id="org51c57f2" class="figure">
<p><img src="figs/plant_uncertainty_bode_plot.png" alt="plant_uncertainty_bode_plot.png" />
</p>
<p><span class="figure-number">Figure 5: </span>Some elements in the uncertainty set \(\Pi_I\)</p>
</div>
</div>
</div>
<div id="outline-container-orgc908813" class="outline-2">
<h2 id="orgc908813"><span class="section-number-2">3</span> Specifications and performance weights</h2>
<div class="outline-text-2" id="text-3">
<p>
<a id="org3dfa664"></a>
</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>
These specifications can be represented as upper bounds on the closed loop transfer functions \(S\) and \(T\) (see figure <a href="#org8dcc955">6</a>).
</p>
<div id="org8dcc955" class="figure">
<p><img src="figs/bode_requirements.png" alt="bode_requirements.png" />
</p>
<p><span class="figure-number">Figure 6: </span>Upper bounds on \(S\) and \(T\)</p>
</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>
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>.
</p>
<div class="org-src-container">
<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;
</pre>
</div>
<div id="org702c691" class="figure">
<p><img src="figs/compare_weights_upper_bounds_S_T.png" alt="compare_weights_upper_bounds_S_T.png" />
</p>
<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>
</div>
</div>
</div>
<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>
<div class="outline-text-2" id="text-4">
<p>
<a id="orgacb23fe"></a>
</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>
We plot these conditions on figure <a href="#org3739909">8</a>.
</p>
<div id="org3739909" class="figure">
<p><img src="figs/weights_NP_RS_RP.png" alt="weights_NP_RS_RP.png" />
</p>
<p><span class="figure-number">Figure 8: </span>Upper bounds on the norm of the complementary filters for NP, RS and RP</p>
</div>
</div>
</div>
<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>
<div class="outline-text-2" id="text-5">
<p>
<a id="org4f2901c"></a>
</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>
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.
</p>
<div id="org3fbed5b" class="figure">
<p><img src="figs-paper/sf_hinf_filters_plant_b.png" alt="sf_hinf_filters_plant_b.png" />
</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>
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:
\[ \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>
<div id="org3939c2e" class="figure">
<p><img src="figs-paper/sf_hinf_filters_b.png" alt="sf_hinf_filters_b.png" />
</p>
<p><span class="figure-number">Figure 10: </span>\(\mathcal{H}_\infty\) synthesis of the complementary filters</p>
</div>
<p>
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.
</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">
<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;
</pre>
</div>
<div id="org90a3139" class="figure">
<p><img src="figs/weights_wl_wh.png" alt="weights_wl_wh.png" />
</p>
<p><span class="figure-number">Figure 11: </span>Weights on the complementary filters \(w_L\) and \(w_H\) and the associated performance weights</p>
</div>
<p>
We define the generalized plant \(P\) on matlab.
</p>
<div class="org-src-container">
<pre class="src src-matlab">P = [0 wL;
wH <span class="org-type">-</span>wH;
1 0];
</pre>
</div>
<p>
And we do the \(\mathcal{H}_\infty\) synthesis using the <code>hinfsyn</code> command.
</p>
<div class="org-src-container">
<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>);
</pre>
</div>
<pre class="example" id="org67390b9">
[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>
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>.
</p>
<div class="org-src-container">
<pre class="src src-matlab">Hh_hinf = 1 <span class="org-type">-</span> Hl_hinf;
</pre>
</div>
<div id="orgebb3cd3" class="figure">
<p><img src="figs/hinf_filters_results.png" alt="hinf_filters_results.png" />
</p>
<p><span class="figure-number">Figure 12: </span>Obtained complementary filters using \(\mathcal{H}_\infty\) synthesis</p>
</div>
</div>
</div>
<div id="outline-container-orga721ae7" class="outline-2">
<h2 id="orga721ae7"><span class="section-number-2">6</span> Complementary filters using analytical formula</h2>
<div class="outline-text-2" id="text-6">
<p>
<a id="org1e02524"></a>
</p>
<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>
<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>
This is illustrated on figure <a href="#org72b4c79">13</a>.
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>
<div id="org72b4c79" class="figure">
<p><img src="figs/comp_filters_param_alpha.png" alt="comp_filters_param_alpha.png" />
</p>
<p><span class="figure-number">Figure 13: </span>Effect of the parameter \(\alpha\) on the shape of the generated second order complementary filters</p>
</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>
The obtained complementary filters are shown on figure <a href="#org8e27bc0">14</a>.
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">
<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;
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));
</pre>
</div>
<div id="org8e27bc0" class="figure">
<p><img src="figs/complementary_filters_second_order.png" alt="complementary_filters_second_order.png" />
</p>
<p><span class="figure-number">Figure 14: </span>Second order complementary filters using the analytical formula</p>
</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>
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.
</p>
<div class="org-src-container">
<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;
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));
</pre>
</div>
<div id="org9e9b7f8" class="figure">
<p><img src="figs/complementary_filters_third_order.png" alt="complementary_filters_third_order.png" />
</p>
<p><span class="figure-number">Figure 15: </span>Third order complementary filters using the analytical formula</p>
</div>
</div>
</div>
<div id="outline-container-orgaf265db" class="outline-2">
<h2 id="orgaf265db"><span class="section-number-2">7</span> Comparison of complementary filters</h2>
<div class="outline-text-2" id="text-7">
<p>
<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>.
</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>
<div id="org0e60c80" class="figure">
<p><img src="figs/comp_hinf_analytical.png" alt="comp_hinf_analytical.png" />
</p>
<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>
</div>
</div>
</div>
<div id="outline-container-orga1bd48e" class="outline-2">
<h2 id="orga1bd48e"><span class="section-number-2">8</span> Controller Analysis</h2>
<div class="outline-text-2" id="text-8">
<p>
<a id="org06b9c47"></a>
</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">
<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));
</pre>
</div>
<pre class="example" id="org4380854">
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>
The bode plot of the controller is shown on figure <a href="#orga448a8a">17</a>:
</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>
<div id="orga448a8a" class="figure">
<p><img src="figs/bode_plot_controller.png" alt="bode_plot_controller.png" />
</p>
<p><span class="figure-number">Figure 17: </span>Bode plot of the obtained controller \(K\)</p>
</div>
</div>
</div>
<div id="outline-container-org1eb8c07" class="outline-2">
<h2 id="org1eb8c07"><span class="section-number-2">9</span> Nominal Stability and Nominal Performance</h2>
<div class="outline-text-2" id="text-9">
<p>
<a id="org55c5c37"></a>
</p>
<p>
The nominal stability of the system is first checked with the <code>allmargin</code> matlab command.
</p>
<div class="org-src-container">
<pre class="src src-matlab">allmargin(K<span class="org-type">*</span>G<span class="org-type">*</span>Hl_hinf)
</pre>
</div>
<pre class="example" id="orgab687e2">
allmargin(K*G*Hl_hinf)
ans =
struct with fields:
GainMargin: 3.86196691340257
GMFrequency: 225.992831144248
PhaseMargin: 34.0890006269983
PMFrequency: 88.3632737129296
DelayMargin: 0.0067331740287082
DMFrequency: 88.3632737129296
Stable: 1
</pre>
<p>
The system is stable and the stability margins are good.
</p>
<p>
The bode plot of the loop gain \(L = K*G*H_L\) is shown on figure <a href="#orgcce877a">18</a>.
</p>
<div id="orgcce877a" class="figure">
<p><img src="figs/bode_plot_loop_gain.png" alt="bode_plot_loop_gain.png" />
</p>
<p><span class="figure-number">Figure 18: </span>Bode Plot of the Loop Gain \(L = K G H_L\)</p>
</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">
<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);
</pre>
</div>
<p>
We then compare their norms with the upper bounds on the performance of the system (figure <a href="#orgbbc226a">19</a>).
As expected, we guarantee the Nominal Performance of the system.
</p>
<div id="orgbbc226a" class="figure">
<p><img src="figs/verification_NP.png" alt="verification_NP.png" />
</p>
<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>
</div>
</div>
</div>
<div id="outline-container-org8c33736" class="outline-2">
<h2 id="org8c33736"><span class="section-number-2">10</span> Robust Stability and Robust Performance</h2>
<div class="outline-text-2" id="text-10">
<p>
<a id="org4d2c1c6"></a>
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>
This is shown on figure <a href="#orga0177a7">20</a>.
</p>
<div id="orga0177a7" class="figure">
<p><img src="figs/robust_stability.png" alt="robust_stability.png" />
</p>
<p><span class="figure-number">Figure 20: </span>Robust Stability Check: \(|w_I T| < 1, \quad \forall \omega\)</p>
</div>
<p>
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>).
</p>
<div id="orgfb25f29" class="figure">
<p><img src="figs/loop_gain_robustness.png" alt="loop_gain_robustness.png" />
</p>
<p><span class="figure-number">Figure 21: </span>Loop Gain of the uncertain system</p>
</div>
<div id="orgf81e346" class="figure">
<p><img src="figs/nyquist_robustness.png" alt="nyquist_robustness.png" />
</p>
<p><span class="figure-number">Figure 22: </span>Nyquist plot of the uncertain system</p>
</div>
<p>
The Robust Performance is verified by plotting \(|S|\) and \(|T|\) for the uncertain system along side the upper bounds defined for performance.
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.
</p>
<div id="org58bb0a6" class="figure">
<p><img src="figs/robust_performance_result.png" alt="robust_performance_result.png" />
</p>
<p><span class="figure-number">Figure 23: </span>Verification of the Robust Performance</p>
</div>
</div>
</div>
<div id="outline-container-org367d76d" class="outline-2">
<h2 id="org367d76d"><span class="section-number-2">11</span> Pre-filter</h2>
<div class="outline-text-2" id="text-11">
<p>
<a id="org7ea3388"></a>
</p>
<p>
For now, we have not specified any performance requirements on the input usage due to change of reference.
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>.
</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">
<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);
</pre>
</div>
<p>
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>).
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">
<pre class="src src-matlab">t = linspace(0, 0.01, 500);
opts = stepDataOptions;
opts.StepAmplitude = 1e<span class="org-type">-</span>6;
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);
</pre>
</div>
<div id="orgb70a4d3" class="figure">
<p><img src="figs/u_and_y_with_Kr.png" alt="u_and_y_with_Kr.png" />
</p>
<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>
</div>
</div>
</div>
<div id="outline-container-org846e5c4" class="outline-2">
<h2 id="org846e5c4"><span class="section-number-2">12</span> Matlab Functions</h2>
<div class="outline-text-2" id="text-12">
<p>
<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>
</p>
<p>
This Matlab function is accessible <a href="src/createWeight.m">here</a>.
</p>
<div class="org-src-container">
<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>
</pre>
</div>
</div>
</div>
<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">
<p>
<a id="orgc93a4bb"></a>
</p>
<p>
This Matlab function is accessible <a href="src/plotMagUncertainty.m">here</a>.
</p>
<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>
</div>
<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">
<p>
<a id="org9d5c371"></a>
</p>
<p>
This Matlab function is accessible <a href="src/plotPhaseUncertainty.m">here</a>.
</p>
<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>
</div>
</div>
</div>
</div>
<div id="postamble" class="status">
<p class="author">Author: Thomas Dehaeze</p>
<p class="date">Created: 2020-11-12 jeu. 10:43</p>
</div>
</body>
</html>