Update all the analysis about robustness / uncertainty reduction

matlab/index.org

@@ -44,7 +44,7 @@ Then, three design methods for generating two complementary filters are proposed
 - in section [[sec:feedback_generate_comp_filters]], the classical feedback architecture is used
 - in section [[sec:analytical_formula_literature]], analytical formulas found in the literature are listed
 
-* Optimal Sensor Fusion for noise characteristics
+* Optimal Sensor Fusion - Minimize the Super Sensor Noise
   :PROPERTIES:
   :header-args:matlab+: :tangle matlab/optimal_comp_filters.m
   :header-args:matlab+: :comments org :mkdirp yes
@@ -625,7 +625,7 @@ The RMS value of the obtained super sensors are shown on table [[tab:rms_results
 ** Conclusion
 From the above complementary filter design with the $\mathcal{H}_2$ and $\mathcal{H}_\infty$ synthesis, it still seems that the $\mathcal{H}_2$ synthesis gives the complementary filters that permits to obtain the minimal super sensor noise (when measuring with the $\mathcal{H}_2$ norm).
 
-* Robustness to sensor dynamics uncertainty
+* Optimal Sensor Fusion - Minimize the Super Sensor Dynamical Uncertainty
   :PROPERTIES:
   :header-args:matlab+: :tangle matlab/comp_filter_robustness.m
   :header-args:matlab+: :comments org :mkdirp yes
@@ -633,28 +633,18 @@ From the above complementary filter design with the $\mathcal{H}_2$ and $\mathca
   <<sec:comp_filter_robustness>>
 
 ** Introduction                                                      :ignore:
-Let's first consider ideal sensors where $G_1 = 1$ and $G_2 = 1$ (figure [[fig:fusion_two_noisy_sensors_with_dyn_bis]]).
+We initially considered perfectly known sensor dynamics so that it can be perfectly inverted.
 
-#+name: fig:fusion_two_noisy_sensors_with_dyn_bis
-#+caption: Fusion of two sensors
-[[file:figs/fusion_two_noisy_sensors_with_dyn_bis.png]]
+We now take into account the fact that the sensor dynamics is only partially known.
+To do so, we model the uncertainty that we have on the sensor dynamics by multiplicative input uncertainty as shown in Fig. [[fig:sensor_fusion_dynamic_uncertainty]].
 
-# #+name: fig:fusion_two_signals
-# #+caption: Fusion of two noisy measurements of $x$
-# [[file:figs/fusion_two_signals.png]]
+#+name: fig:sensor_fusion_dynamic_uncertainty
+#+caption: Sensor fusion architecture with sensor dynamics uncertainty
+[[file:figs-tikz/sensor_fusion_dynamic_uncertainty.png]]
 
-We then have:
-\begin{align*}
-  \hat{x} &= (x + n_1) H_1 + (x + n_2) H_2 \\
-          &= x + n_1 H_1 + n_2 H_2
-\end{align*}
+The objective here is to design complementary filters $H_1(s)$ and $H_2(s)$ in order to minimize the dynamical uncertainty of the super sensor.
 
-So the estimation error is
-\[ \delta_x = \hat{x} - x = n_1 H_1 + n_2 H_2 \]
-
-And we see that the complementary filters are only shaping the noise and that they do not impact the transfer function from $x$ to $\hat{x}$ that is in the feedback path.
-
-** ZIP file containing the data and matlab files                     :ignore:
+** ZIP file containing the data and matlab files                    :ignore:
 #+begin_src bash :exports none :results none
   if [ matlab/comp_filter_robustness.m -nt data/comp_filter_robustness.zip ]; then
     cp matlab/comp_filter_robustness.m comp_filter_robustness.m;
@@ -677,38 +667,22 @@ And we see that the complementary filters are only shaping the noise and that th
   <<matlab-init>>
 #+end_src
 
-** Static Gain Mismatch                                            :noexport:
-Even though we here still consider that the two sensors have perfect dynamics, we consider gain mismatch for the two sensors:
-\begin{align*}
-  G_1(s) &= 1 + \delta_1(s) \\
-  G_2(s) &= 1 + \delta_2(s)
-\end{align*}
-
-Thus, we have:
-\begin{align*}
-  \hat{x} &= (x + n_1) (1 + \delta_1) H_1 + (x + n_2) (1 + \delta_2) H_2 \\
-          &= x (1 + \delta_1 H_1 + \delta_2 H_2) + n_1 (1 + \delta_1) H_1 + n_2(1 + \delta_2)  H_2
-\end{align*}
-
-So the transfer function from $x$ to $\hat{x}$ is:
-\begin{align*}
-  \frac{\hat{x}}{x} &= 1 + \delta_1 H_1 + \delta_2 H_2 \\
-                    &= 1 + \delta_1 H_1 + \delta_2 (1 - H_1) \\
-                    &= 1 + (\delta_1 - \delta_2) H_1 + \delta_2 \\
-\end{align*}
+#+begin_src matlab
+  freqs = logspace(-1, 3, 1000);
+#+end_src
 
 ** Unknown sensor dynamics dynamics
 In practical systems, the sensor dynamics has always some level of uncertainty.
-Let's represent that with multiplicative input uncertainty as shown on figure [[fig:fusion_gain_mismatch]].
+Let's represent that with multiplicative input uncertainty as shown on figure [[fig:sensor_fusion_dynamic_uncertainty]].
 
-#+name: fig:fusion_gain_mismatch
+#+name: fig:sensor_fusion_dynamic_uncertainty
 #+caption: Fusion of two sensors with input multiplicative uncertainty
-[[file:figs-tikz/fusion_gain_mismatch.png]]
+[[file:figs-tikz/sensor_fusion_dynamic_uncertainty.png]]
 
-We have:
+The dynamics of the super sensor is represented by
 \begin{align*}
-  \frac{\hat{x}}{x} &= (1 + W_1 \Delta_1) H_1 + (1 + W_2 \Delta_2) H_2 \\
-                    &= 1 + W_1 H_1 \Delta_1 + W_2 H_2 \Delta_2
+  \frac{\hat{x}}{x} &= (1 + w_1 \Delta_1) H_1 + (1 + w_2 \Delta_2) H_2 \\
+                    &= 1 + w_1 H_1 \Delta_1 + w_2 H_2 \Delta_2
 \end{align*}
 With $\Delta_i$ is any transfer function satisfying $\| \Delta_i \|_\infty < 1$.
 
@@ -741,31 +715,35 @@ We then have that the angle introduced by the super sensor is bounded by $\arcsi
 Thus, we choose should choose $\epsilon$ so that the maximum phase uncertainty introduced by the sensors is of an acceptable value.
 
 ** Design the complementary filters in order to limit the phase and gain uncertainty of the super sensor
-Let's say the two sensors dynamics $H_1$ and $H_2$ have been identified with the associated uncertainty weights $W_1$ and $W_2$.
+Let's say the two sensors dynamics have been identified with the associated uncertainty weights $w_1(s)$ and $w_2(s)$.
 
-If we want to have a maximum phase introduced by the sensors of 20 degrees, we have to design $H_1$ and $H_2$ such that:
+If we want to have a maximum phase introduced by the sensors of 20 degrees, we have to design $H_1(s)$ and $H_2(s)$ such that:
 \begin{align*}
-                      & arcsin(|H_1 W_1| + |H_2 W_2|) < 20 \text{ deg} \\
-  \Longleftrightarrow & |H_1 W_1| + |H_2 W_2| < 0.34
+                      & arcsin\Big( |H_1(j\omega) w_1(j\omega)| + |H_2(j\omega) w_2(j\omega)| \Big) < 20 \text{ deg} \quad \forall\omega \\
+  \Longleftrightarrow & |H_1 w_1| + |H_2 w_2| < 0.34 \quad \forall\omega
 \end{align*}
 
 We can do that with the $\mathcal{H}_\infty$ synthesis by setting upper bounds on the complementary filters using weights that corresponds to the sensor dynamics uncertainty.
 
-For simplicity, let's suppose $W_1(s) = W_2(s) = 0.1$ ($10\%$ uncertainty in the sensor gain).
-\[ |H_1 W_1| + |H_2 W_2| < 3.4 \]
-
+Basically, at frequencies where $|w_i(j\omega)|$ is large, $|H_i(j\omega)|$ has to be made small.
 Thus, by limiting the norm of the complementary filters, we can limit the maximum unwanted phase introduced by the uncertainty on the sensors dynamics.
 
 This is of primary importance in order to ensure the stability of the feedback loop using the super sensor signal.
 
+
 ** First Basic Example with gain mismatch
-Let's consider two ideal sensors except one sensor has not an expected gain of one but a gain of $0.6$.
+Let's consider two ideal sensors except one sensor has not an expected unity gain but a gain equal to $0.6$:
+\begin{align*}
+  G_1(s) &= 1 \\
+  G_2(s) &= 0.6
+\end{align*}
+
 #+begin_src matlab
   G1 = 1;
   G2 = 0.6;
 #+end_src
 
-Let's design two complementary filters as shown on figure [[fig:comp_filters_robustness_test]].
+Two pairs of complementary filters are designed and shown on figure [[fig:comp_filters_robustness_test]].
 The complementary filters shown in blue does not present a bump as the red ones but provides less sensor separation at high and low frequencies.
 
 #+begin_src matlab :exports none
@@ -822,7 +800,7 @@ The complementary filters shown in blue does not present a bump as the red ones
 #+CAPTION: The two complementary filters designed for the robustness test ([[./figs/comp_filters_robustness_test.png][png]], [[./figs/comp_filters_robustness_test.pdf][pdf]])
 [[file:figs/comp_filters_robustness_test.png]]
 
-We then compute the bode plot of the super sensor transfer function $H_1*G_1 + H_2*G_2$ for both complementary filters pair (figure [[fig:tf_super_sensor_comp]]).
+We then compute the bode plot of the super sensor transfer function $H_1(s)G_1(s) + H_2(s)G_2(s)$ for both complementary filters pair (figure [[fig:tf_super_sensor_comp]]).
 
 We see that the blue complementary filters with a lower maximum norm permits to limit the phase lag introduced by the gain mismatch.
 
@@ -863,7 +841,367 @@ We see that the blue complementary filters with a lower maximum norm permits to
 #+CAPTION: Comparison of the obtained super sensor transfer functions ([[./figs/tf_super_sensor_comp.png][png]], [[./figs/tf_super_sensor_comp.pdf][pdf]])
 [[file:figs/tf_super_sensor_comp.png]]
 
-** TODO More Complete example with model uncertainty
+** More Complete example with dynamical uncertainty
+We want to merge two sensors:
+- sensor 1 that has unknown dynamics above 10Hz: $|w_1(j\omega)| > 1$ for $\omega > 10\text{ Hz}$
+- sensor 2 that has unknown dynamics below 1Hz and above 1kHz $|w_2(j\omega)| > 1$ for $\omega < 1\text{ Hz}$ and $\omega > 1\text{ kHz}$
+
+*** Dynamical uncertainty of the individual sensors
+We define the weights that are used to characterize the dynamic uncertainty of the sensors.
+#+begin_src matlab
+  omegac = 100*2*pi; G0 = 0.1; Ginf = 10;
+  w1 = (Ginf*s/omegac + G0)/(s/omegac + 1);
+
+  omegac = 0.2*2*pi; G0 = 5; Ginf = 0.1;
+  w2 = (Ginf*s/omegac + G0)/(s/omegac + 1);
+  omegac = 5000*2*pi; G0 = 1; Ginf = 50;
+  w2 = w2*(Ginf*s/omegac + G0)/(s/omegac + 1);
+#+end_src
+
+From the weights, we define the uncertain transfer functions of the sensors. Some of the uncertain dynamics of both sensors are shown on Fig. [[fig:uncertainty_dynamics_sensors]] with the upper and lower bounds on the magnitude and on the phase.
+#+begin_src matlab
+  G1 = 1 + w1*ultidyn('Delta',[1 1]);
+  G2 = 1 + w2*ultidyn('Delta',[1 1]);
+#+end_src
+
+#+begin_src matlab :exports none
+  % Few random samples of the sensor dynamics are computed
+  G1s = usample(G1, 10);
+  G2s = usample(G2, 10);
+#+end_src
+
+#+begin_src matlab :exports none
+  % We here compute the maximum and minimum phase of both sensors
+  Dphi1 = 180/pi*asin(abs(squeeze(freqresp(w1, freqs, 'Hz'))));
+  Dphi2 = 180/pi*asin(abs(squeeze(freqresp(w2, freqs, 'Hz'))));
+  Dphi1(abs(squeeze(freqresp(w1, freqs, 'Hz'))) > 1) = 190;
+  Dphi2(abs(squeeze(freqresp(w2, freqs, 'Hz'))) > 1) = 190;
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  % Magnitude
+  ax1 = subaxis(2,1,1);
+  hold on;
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, 1 + abs(squeeze(freqresp(w1, freqs, 'Hz'))), '--');
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, max(1 - abs(squeeze(freqresp(w1, freqs, 'Hz'))), 0), '--');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs, 1 + abs(squeeze(freqresp(w2, freqs, 'Hz'))), '--');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs, max(1 - abs(squeeze(freqresp(w2, freqs, 'Hz'))), 0), '--');
+  for i = 1:length(G1s)
+    plot(freqs, abs(squeeze(freqresp(G1s(:, :, i, 1), freqs, 'Hz'))), '-', 'color', [0 0.4470 0.7410 0.4]);
+    plot(freqs, abs(squeeze(freqresp(G2s(:, :, i, 1), freqs, 'Hz'))), '-', 'color', [0.8500 0.3250 0.0980 0.4]);
+  end
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Magnitude');
+  ylim([1e-1, 10]);
+  hold off;
+
+  % Phase
+  ax2 = subaxis(2,1,2);
+  hold on;
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs,  Dphi1, '--');
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, -Dphi1, '--');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs,  Dphi2, '--');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs, -Dphi2, '--');
+  for i = 1:length(G1s)
+    plot(freqs, 180/pi*angle(squeeze(freqresp(G1s(:, :, i, 1), freqs, 'Hz'))), '-', 'color', [0 0.4470 0.7410 0.4]);
+    plot(freqs, 180/pi*angle(squeeze(freqresp(G2s(:, :, i, 1), freqs, 'Hz'))), '-', 'color', [0.8500 0.3250 0.0980 0.4]);
+  end
+  set(gca,'xscale','log');
+  yticks(-180:90:180);
+  ylim([-180 180]);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  hold off;
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/uncertainty_dynamics_sensors.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:uncertainty_dynamics_sensors
+#+CAPTION: Dynamic uncertainty of the two sensors ([[./figs/uncertainty_dynamics_sensors.png][png]], [[./figs/uncertainty_dynamics_sensors.pdf][pdf]])
+[[file:figs/uncertainty_dynamics_sensors.png]]
+
+*** Synthesis objective
+The uncertainty region of the super sensor dynamics is represented by a circle in the complex plane as shown in Fig. [[fig:uncertainty_gain_phase_variation]].
+
+At each frequency $\omega$, the radius of the circle is $|w_1(j\omega) H_1(j\omega)| + |w_2(j\omega) H_2(j\omega)|$.
+
+Thus, the phase shift $\Delta\phi(\omega)$ due to the super sensor uncertainty is bounded by:
+\[ |\Delta\phi(\omega)| \leq \arcsin\big( |w_1(j\omega) H_1(j\omega)| + |w_2(j\omega) H_2(j\omega)| \big) \]
+
+Let's define some allowed frequency depend phase shift $\Delta\phi_\text{max}(\omega) > 0$ such that:
+\[ |\Delta\phi(\omega)| < \Delta\phi_\text{max}(\omega), \quad \forall\omega \]
+
+
+If $H_1(s)$ and $H_2(s)$ are designed such that
+\[ |w_1(j\omega) H_1(j\omega)| + |w_2(j\omega) H_2(j\omega)| < \sin\big( \Delta\phi_\text{max}(\omega) \big) \]
+
+The maximum phase shift due to dynamic uncertainty at frequency $\omega$ will be $\Delta\phi_\text{max}(\omega)$.
+
+*** Requirements as an $\mathcal{H}_\infty$ norm
+We know try to express this requirement in terms of an $\mathcal{H}_\infty$ norm.
+
+Let define one weight $w_\phi(s)$ that represents the maximum wanted phase uncertainty:
+\[ |w_{\phi}(j\omega)|^{-1} \approx \sin(\Delta\phi_{\text{max}}(\omega)), \quad \forall\omega \]
+
+Then:
+\begin{align*}
+                      & |w_1(j\omega) H_1(j\omega)| + |w_2(j\omega) H_2(j\omega)| < \sin\big( \Delta\phi_\text{max}(\omega) \big), \quad \forall\omega \\
+  \Longleftrightarrow & |w_1(j\omega) H_1(j\omega)| + |w_2(j\omega) H_2(j\omega)| < |w_\phi(j\omega)|^{-1}, \quad \forall\omega \\
+  \Longleftrightarrow & \left| w_1(j\omega) H_1(j\omega) w_\phi(j\omega) \right| + \left| w_2(j\omega) H_2(j\omega) w_\phi(j\omega) \right| < 1, \quad \forall\omega
+\end{align*}
+
+Which is approximately equivalent to (with an approximation of maximum $\sqrt{2}$):
+#+name: eq:hinf_conf_phase_uncertainty
+\begin{equation}
+  \left\| \begin{matrix} w_1(s) w_\phi(s) H_1(s) \\ w_2(s) w_\phi(s) H_2(s) \end{matrix} \right\|_\infty < 1
+\end{equation}
+
+On should not forget that at frequency where both sensors has unknown dynamics ($|w_1(j\omega)| > 1$ and $|w_2(j\omega)| > 1$), the super sensor dynamics will also be unknown and the phase uncertainty cannot be bounded.
+Thus, at these frequencies, $|w_\phi|$ should be smaller than $1$.
+
+*** H-Infinity Synthesis
+Let's define $w_\phi(s)$ such that the maximum allowed phase uncertainty $\Delta\phi_\text{max}$ of the super sensor dynamics is $30\text{ deg}$ until frequency $\omega_0 = 500\text{ Hz}$
+
+#+begin_src matlab
+  Dphi = 20; % [deg]
+
+  n = 4; w0 = 2*pi*900; G0 = 1/sin(Dphi*pi/180); Ginf = 1/100; Gc = 1;
+  wphi = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/Ginf)^(2/n)))*s + (G0/Gc)^(1/n))/((1/Ginf)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/Ginf)^(2/n)))*s + (1/Gc)^(1/n)))^n;
+
+  W1 = w1*wphi;
+  W2 = w2*wphi;
+#+end_src
+
+#+begin_src matlab :exports none
+  % We here compute the wanted maximum and minimum phase of the super sensor
+  Dphimax = 180/pi*asin(1./abs(squeeze(freqresp(wphi, freqs, 'Hz'))));
+  Dphimax(1./abs(squeeze(freqresp(wphi, freqs, 'Hz'))) > 1) = 190;
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(freqs,  Dphimax, 'k--');
+  plot(freqs, -Dphimax, 'k--');
+  set(gca, 'XScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Magnitude');
+  hold off;
+  xlim([freqs(1), freqs(end)]);
+  ylim([-180 180]);
+  yticks(-180:45:180);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/maximum_wanted_phase_uncertainty.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:maximum_wanted_phase_uncertainty
+#+CAPTION: Maximum wanted phase uncertainty using this weight ([[./figs/maximum_wanted_phase_uncertainty.png][png]], [[./figs/maximum_wanted_phase_uncertainty.pdf][pdf]])
+[[file:figs/maximum_wanted_phase_uncertainty.png]]
+
+The obtained upper bounds on the complementary filters in order to limit the phase uncertainty of the super sensor are represented in Fig. [[fig:upper_bounds_comp_filter_max_phase_uncertainty]].
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '-', 'DisplayName', '$1/|w_1w_\phi|$');
+  plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '-', 'DisplayName', '$1/|w_2w_\phi|$');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Magnitude');
+  hold off;
+  xlim([freqs(1), freqs(end)]);
+  legend('location', 'northeast');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/upper_bounds_comp_filter_max_phase_uncertainty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:upper_bounds_comp_filter_max_phase_uncertainty
+#+CAPTION: Upper bounds on the complementary filters set in order to limit the maximum phase uncertainty of the super sensor to 30 degrees until 500Hz ([[./figs/upper_bounds_comp_filter_max_phase_uncertainty.png][png]], [[./figs/upper_bounds_comp_filter_max_phase_uncertainty.pdf][pdf]])
+[[file:figs/upper_bounds_comp_filter_max_phase_uncertainty.png]]
+
+The $\mathcal{H}_\infty$ synthesis is performed using the defined weights and the obtained complementary filters are shown in Fig. [[fig:comp_filter_hinf_uncertainty]].
+
+#+begin_src matlab
+  P = [W1 -W1;
+       0   W2;
+       1   0];
+#+end_src
+
+And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
+#+begin_src matlab :results output replace :exports both
+  [H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
+#+end_src
+
+#+RESULTS:
+#+begin_example
+[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
+Resetting value of Gamma min based on D_11, D_12, D_21 terms
+
+Test bounds:      0.0447 <  gamma  <=      1.3318
+
+  gamma    hamx_eig  xinf_eig  hamy_eig   yinf_eig   nrho_xy   p/f
+    1.332   1.3e+01 -1.0e-14   1.3e+00   -2.6e-18    0.0000    p
+    0.688   1.3e-11# ********   1.3e+00   -6.7e-15  ********    f
+    1.010   1.1e+01 -1.5e-14   1.3e+00   -2.5e-14    0.0000    p
+    0.849   6.9e-11# ********   1.3e+00   -2.3e-14  ********    f
+    0.930   5.2e-12# ********   1.3e+00   -6.1e-18  ********    f
+    0.970   5.6e-11# ********   1.3e+00   -2.3e-14  ********    f
+    0.990   5.0e-11# ********   1.3e+00   -1.7e-17  ********    f
+    1.000   2.1e-10# ********   1.3e+00    0.0e+00  ********    f
+    1.005   1.9e-10# ********   1.3e+00   -3.7e-14  ********    f
+    1.008   1.1e+01 -9.1e-15   1.3e+00    0.0e+00    0.0000    p
+    1.006   1.2e-09# ********   1.3e+00   -6.9e-16  ********    f
+    1.007   1.1e+01 -4.6e-15   1.3e+00   -1.8e-16    0.0000    p
+
+ Gamma value achieved:     1.0069
+#+end_example
+
+#+begin_src matlab
+  H1 = 1 - H2;
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+
+  ax1 = subplot(2,1,1);
+  hold on;
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$W_1$');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$W_2$');
+
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
+
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude');
+  set(gca, 'XTickLabel',[]);
+  legend('location', 'northeast');
+
+  ax2 = subplot(2,1,2);
+  hold on;
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-');
+  hold off;
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'XScale', 'log');
+  yticks([-360:90:360]);
+
+  linkaxes([ax1,ax2],'x');
+  xlim([freqs(1), freqs(end)]);
+  xticks([0.1, 1, 10, 100, 1000]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/comp_filter_hinf_uncertainty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:comp_filter_hinf_uncertainty
+#+CAPTION: Obtained complementary filters ([[./figs/comp_filter_hinf_uncertainty.png][png]], [[./figs/comp_filter_hinf_uncertainty.pdf][pdf]])
+[[file:figs/comp_filter_hinf_uncertainty.png]]
+
+*** Super sensor uncertainty
+We can now compute the uncertainty of the super sensor. The result is shown in Fig. [[fig:super_sensor_uncertainty_bode_plot]].
+
+#+begin_src matlab
+  Gss = G1*H1 + G2*H2;
+#+end_src
+
+#+begin_src matlab :exports none
+  Gsss = usample(Gss, 20);
+#+end_src
+
+#+begin_src matlab :exports none
+  % We here compute the maximum and minimum phase of the super sensor
+  Dphiss = 180/pi*asin(abs(squeeze(freqresp(w1*H1, freqs, 'Hz')))+abs(squeeze(freqresp(w2*H2, freqs, 'Hz'))));
+  Dphiss(abs(squeeze(freqresp(w1*H1, freqs, 'Hz')))+abs(squeeze(freqresp(w2*H2, freqs, 'Hz'))) > 1) = 190;
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  % Magnitude
+  ax1 = subaxis(2,1,1);
+  hold on;
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, 1 + abs(squeeze(freqresp(w1, freqs, 'Hz'))), '--');
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, max(1 - abs(squeeze(freqresp(w1, freqs, 'Hz'))), 0), '--');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs, 1 + abs(squeeze(freqresp(w2, freqs, 'Hz'))), '--');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs, max(1 - abs(squeeze(freqresp(w2, freqs, 'Hz'))), 0), '--');
+  plot(freqs, 1 + abs(squeeze(freqresp(w1*H1+w2*H2, freqs, 'Hz'))), 'k--');
+  plot(freqs, max(1 - abs(squeeze(freqresp(w1*H1+w2*H2, freqs, 'Hz'))), 0), 'k--');
+  for i = 1:length(Gsss)
+    plot(freqs, abs(squeeze(freqresp(Gsss(:, :, i, 1), freqs, 'Hz'))), '-', 'color', [0 0 0 0.2]);
+  end
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Magnitude');
+  ylim([1e-1, 10]);
+  hold off;
+
+  % Phase
+  ax2 = subaxis(2,1,2);
+  hold on;
+  % plot(freqs,  Dphimax, 'r-');
+  % plot(freqs, -Dphimax, 'r-');
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs,  Dphi1, '--');
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, -Dphi1, '--');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs,  Dphi2, '--');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs, -Dphi2, '--');
+  plot(freqs,  Dphiss, 'k--');
+  plot(freqs, -Dphiss, 'k--');
+  for i = 1:length(Gsss)
+    plot(freqs, 180/pi*angle(squeeze(freqresp(Gsss(:, :, i, 1), freqs, 'Hz'))), '-', 'color', [0 0 0 0.2]);
+  end
+  set(gca,'xscale','log');
+  yticks(-180:90:180);
+  ylim([-180 180]);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  hold off;
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/super_sensor_uncertainty_bode_plot.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:super_sensor_uncertainty_bode_plot
+#+CAPTION: Uncertainty on the dynamics of the super sensor ([[./figs/super_sensor_uncertainty_bode_plot.png][png]], [[./figs/super_sensor_uncertainty_bode_plot.pdf][pdf]])
+[[file:figs/super_sensor_uncertainty_bode_plot.png]]
+
+We here just used very wimple weights. We could shape the dynamical uncertainty of the super sensor by using more complex weights.
+We could for instance ask for less uncertainty at low frequency.
+
 * Equivalent Super Sensor
 
 The goal here is to find the parameters of a single sensor that would best represent a super sensor.