Add CVX optimization section

matlab/index.org

@@ -1502,8 +1502,7 @@ We see that the blue complementary filters with a lower maximum norm permits to
 :header-args:matlab+: :comments org :mkdirp yes
 :END:
 <<sec:mixed_synthesis_sensor_fusion>>
-
-** ZIP file containing the data and matlab files                     :ignore:
+** ZIP file containing the data and matlab files                    :ignore:
 #+begin_note
   The Matlab scripts is accessible [[file:matlab/mixed_synthesis_sensor_fusion.m][here]].
 #+end_note
@@ -1859,7 +1858,7 @@ The PSD and CPS of the super sensor's noise are shown in Fig. [[fig:psd_super_se
 [[file:figs/cps_super_sensor_mixed_syn.png]]
 
 ** Obtained Super Sensor's Uncertainty
-The uncertainty on the super sensor's dynamics is shown in Fig. [[]].
+The uncertainty on the super sensor's dynamics is shown in Fig. [[fig:super_sensor_dyn_uncertainty_mixed_syn]].
 
 #+begin_src matlab :exports none
   G1 = 1 + w1*ultidyn('Delta',[1 1]);
@@ -1943,6 +1942,461 @@ This synthesis methods allows both to:
 - limit the dynamical uncertainty of the super sensor
 - minimize the RMS value of the estimation
 
+* Mixed Synthesis - LMI Optimization
+** Introduction
+The following matlab scripts was written by Mohit.
+
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+  <<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  <<matlab-init>>
+#+end_src
+
+#+begin_src matlab
+  freqs = logspace(-1, 3, 1000);
+#+end_src
+
+** Noise characteristics and Uncertainty of the individual sensors
+We define the weights that are used to characterize the dynamic uncertainty of the sensors. This will be used for the $\mathcal{H}_\infty$ part of the synthesis.
+#+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
+
+We define the noise characteristics of the two sensors by choosing $N_1$ and $N_2$. This will be used for the $\mathcal{H}_2$ part of the synthesis.
+#+begin_src matlab
+  omegac = 100*2*pi; G0 = 1e-5; Ginf = 1e-4;
+  N1 = (Ginf*s/omegac + G0)/(s/omegac + 1)/(1 + s/2/pi/100);
+
+  omegac = 1*2*pi; G0 = 1e-3; Ginf = 1e-8;
+  N2 = ((sqrt(Ginf)*s/omegac + sqrt(G0))/(s/omegac + 1))^2/(1 + s/2/pi/4000)^2;
+#+end_src
+
+** Weights
+The weights for the $\mathcal{H}_2$ and $\mathcal{H}_\infty$ part are defined below.
+
+#+begin_src matlab
+  n = 4; w0 = 2*pi*900; G0 = 9; G1 = 1; Gc = 1.1;
+  H = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n;
+  wphi = 0.2*(s+3.142e04)/(s+628.3)/H;
+
+  W1u = ss(w1*wphi); W2u = ss(w2*wphi); % Weight on the uncertainty
+  W1n = ss(N1); W2n = ss(N2); % Weight on the noise
+#+end_src
+
+#+begin_src matlab
+  P = [W1u -W1u;
+       0    W2u;
+       W1n -W1n;
+       0    W2n;
+       1    0];
+#+end_src
+
+** LMI Optimization
+We are using the [[http://cvxr.com/cvx/][CVX toolbox]] to solve the optimization problem.
+
+We first put the generalized plant in a State-space form.
+#+begin_src matlab
+  A = P.A;
+  Bw = P.B(:,1);
+  Bu = P.B(:,2);
+  Cz1 = P.C(1:2,:); Dz1w = P.D(1:2,1); Dz1u = P.D(1:2,2); % Hinf
+  Cz2 = P.C(3:4,:); Dz2w = P.D(1:2,1); Dz2u = P.D(1:2,2); % H2
+  Cy = P.C(5,:); Dyw = P.D(5,1); Dyu = P.D(5,2);
+
+  n = size(P.A,1);
+  ny = 1; % number of measurements
+  nu = 1; % number of control inputs
+  nz = 2;
+  nw = 1;
+
+  Wtinf = 0;
+  Wt2 = 1;
+#+end_src
+
+We Define all the variables.
+#+begin_src matlab
+  cvx_startup;
+
+  cvx_begin sdp
+  cvx_quiet true
+  cvx_solver sedumi
+  variable X(n,n) symmetric;
+  variable Y(n,n) symmetric;
+  variable W(nz,nz) symmetric;
+  variable Ah(n,n);
+  variable Bh(n,ny);
+  variable Ch(nu,n);
+  variable Dh(nu,ny);
+  variable eta;
+  variable gam;
+#+end_src
+
+We define the minimization objective.
+#+begin_src matlab
+  minimize Wt2*eta+Wtinf*gam % mix objective
+  subject to:
+#+end_src
+
+The $\mathcal{H}_\infty$ constraint.
+#+begin_src matlab
+  gam<=1; % Keep the Hinf norm less than 1
+
+  [ X, eye(n,n) ;
+    eye(n,n), Y ] >= 0 ;
+
+  [ A*X + Bu*Ch + X*A' + Ch'*Bu', A+Bu*Dh*Cy+Ah', Bw+Bu*Dh*Dyw, X*Cz1' + Ch'*Dz1u' ;
+    (A+Bu*Dh*Cy+Ah')', Y*A + A'*Y + Bh*Cy + Cy'*Bh', Y*Bw + Bh*Dyw, (Cz1+Dz1u*Dh*Cy)' ;
+    (Bw+Bu*Dh*Dyw)', Bw'*Y + Dyw'*Bh', -eye(nw,nw), (Dz1w+Dz1u*Dh*Dyw)' ;
+    Cz1*X + Dz1u*Ch, Cz1+Dz1u*Dh*Cy, Dz1w+Dz1u*Dh*Dyw, -gam*eye(nz,nz)] <= 0 ;
+#+end_src
+
+The $\mathcal{H}_2$ constraint.
+#+begin_src matlab
+  trace(W) <= eta ;
+
+  [ W, Cz2*X+Dz2u*Ch, Cz2*X+Dz2u*Ch;
+    X*Cz2'+Ch'*Dz2u', X, eye(n,n) ;
+    (Cz2*X+Dz2u*Ch)', eye(n,n), Y ] >= 0 ;
+
+  [ A*X + Bu*Ch + X*A' + Ch'*Bu', A+Bu*Dh*Cy+Ah', Bw+Bu*Dh*Dyw ;
+    (A+Bu*Dh*Cy+Ah')', Y*A + A'*Y + Bh*Cy + Cy'*Bh', Y*Bw + Bh*Dyw ;
+    (Bw+Bu*Dh*Dyw)', Bw'*Y + Dyw'*Bh', -eye(nw,nw)] <= 0 ;
+#+end_src
+
+And we run the optimization.
+#+begin_src matlab
+  cvx_end
+  cvx_status
+#+end_src
+
+#+begin_src matlab :exports none
+  if(strcmp(cvx_status,'Inaccurate/Solved'))
+      display('The solver was unable to make a determination to within the default numerical tolerance.');
+      display('However, it determined that the results obtained satisfied a “relaxed” tolerance leve');
+      display('and therefore may still be suitable for further use.');
+  end
+#+end_src
+
+Finally, we can compute the obtained complementary filters.
+#+begin_src matlab
+  M = eye(n);
+  N = inv(M)*(eye(n,n)-Y*X);
+  Dk = Dh;
+  Ck = (Ch-Dk*Cy*X)*inv(M');
+  Bk = inv(N)*(Bh-Y*Bu*Dk);
+  Ak = inv(N)*(Ah-Y*(A+Bu*Dk*Cy)*X-N*Bk*Cy*X-Y*Bu*Ck*M')*inv(M');
+
+  H2 = tf(ss(Ak,Bk,Ck,Dk));
+  H1 = 1 - H2;
+#+end_src
+
+** Result
+The obtained complementary filters are compared with the required upper bounds on Fig. [[fig:LMI_obtained_comp_filters]].
+
+#+begin_src matlab :exports none
+  figure;
+
+  ax1 = subplot(2,1,1);
+  hold on;
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, 1./abs(squeeze(freqresp(W1u, freqs, 'Hz'))), '--', 'DisplayName', '$W_1$');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, 1./abs(squeeze(freqresp(W2u, 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',[]);
+  ylim([1e-3, 2]);
+  legend('location', 'southwest');
+
+  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]);
+  ylim([-180, 180]);
+
+  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/LMI_obtained_comp_filters.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:LMI_obtained_comp_filters
+#+CAPTION: Obtained complementary filters using the LMI optimization ([[./figs/LMI_obtained_comp_filters.png][png]], [[./figs/LMI_obtained_comp_filters.pdf][pdf]])
+[[file:figs/LMI_obtained_comp_filters.png]]
+
+** Comparison with the matlab Mixed Synthesis
+The Mixed $\mathcal{H}_2$ / $\mathcal{H}_\infty$ Synthesis is performed below.
+#+begin_src matlab
+  Nmeas = 1; Ncon = 1; Nz2 = 2;
+
+  [H2m,~,normz,~] = h2hinfsyn(P, Nmeas, Ncon, Nz2, [0, 1], 'HINFMAX', 1, 'H2MAX', Inf, 'DKMAX', 100, 'TOL', 0.01, 'DISPLAY', 'on');
+
+  H1m = 1 - H2m;
+#+end_src
+
+The obtained filters are compare with the one obtained using the CVX toolbox in Fig. [[]].
+
+#+begin_src matlab :exports none
+  figure;
+
+  ax1 = subplot(2,1,1);
+  hold on;
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, abs(squeeze(freqresp(H1m, freqs, 'Hz'))), '--', 'DisplayName', '$H_{1,\mathcal{H}_2/\mathcal{H}_\infty}$');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, abs(squeeze(freqresp(H2m, freqs, 'Hz'))), '--', 'DisplayName', '$H_{2,\mathcal{H}_2/\mathcal{H}_\infty}$');
+
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_{1, CVX}$');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_{2, CVX}$');
+
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude');
+  set(gca, 'XTickLabel',[]);
+  ylim([1e-3, 2]);
+  legend('location', 'southwest');
+
+  ax2 = subplot(2,1,2);
+  hold on;
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, 180/pi*phase(squeeze(freqresp(H1m, freqs, 'Hz'))), '--');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, 180/pi*phase(squeeze(freqresp(H2m, freqs, 'Hz'))), '--');
+
+  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]);
+  ylim([-180, 180]);
+
+  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/compare_cvx_h2hinf_comp_filters.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:compare_cvx_h2hinf_comp_filters
+#+CAPTION: Comparison between the complementary filters obtained with the CVX toolbox and with the =h2hinfsyn= command ([[./figs/compare_cvx_h2hinf_comp_filters.png][png]], [[./figs/compare_cvx_h2hinf_comp_filters.pdf][pdf]])
+[[file:figs/compare_cvx_h2hinf_comp_filters.png]]
+
+** H-Infinity Objective
+In terms of the $\mathcal{H}_\infty$ objective, both synthesis method are satisfying the requirements as shown in Fig. [[fig:comp_cvx_h2i_hinf_norm]].
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, 1./abs(squeeze(freqresp(W1u, freqs, 'Hz'))), '-.', 'DisplayName', '$1/W_{1u}$');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, 1./abs(squeeze(freqresp(W2u, freqs, 'Hz'))), '-.', 'DisplayName', '$1/W_{2u}$');
+
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, abs(squeeze(freqresp(H1m, freqs, 'Hz'))), '--', 'DisplayName', '$H_{1,\mathcal{H}_2/\mathcal{H}_\infty}$');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, abs(squeeze(freqresp(H2m, freqs, 'Hz'))), '--', 'DisplayName', '$H_{2,\mathcal{H}_2/\mathcal{H}_\infty}$');
+
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_{1, CVX}$');
+  set(gca,'ColorOrderIndex',2)
+  plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_{2, CVX}$');
+
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude');
+  set(gca, 'XTickLabel',[]);
+  ylim([1e-3, 2]);
+  legend('location', 'southwest');
+  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_cvx_h2i_hinf_norm.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:comp_cvx_h2i_hinf_norm
+#+CAPTION: H-Infinity norm requirement and results ([[./figs/comp_cvx_h2i_hinf_norm.png][png]], [[./figs/comp_cvx_h2i_hinf_norm.pdf][pdf]])
+[[file:figs/comp_cvx_h2i_hinf_norm.png]]
+
+** Obtained Super Sensor's noise
+The PSD and CPS of the super sensor's noise obtained with the CVX toolbox and =h2hinfsyn= command are compared in Fig. [[fig:psd_compare_cvx_h2i]] and [[fig:cps_compare_cvx_h2i]].
+
+#+begin_src matlab :exports none
+  PSD_S1 = abs(squeeze(freqresp(N1, freqs, 'Hz'))).^2;
+  PSD_S2 = abs(squeeze(freqresp(N2, freqs, 'Hz'))).^2;
+  PSD_cvx = abs(squeeze(freqresp(N1*H1, freqs, 'Hz'))).^2+abs(squeeze(freqresp(N2*H2, freqs, 'Hz'))).^2;
+  PSD_h2i = abs(squeeze(freqresp(N1*H1m, freqs, 'Hz'))).^2+abs(squeeze(freqresp(N2*H2m, freqs, 'Hz'))).^2;
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(freqs, PSD_S1, '-',  'DisplayName', '$\Phi_{\hat{x}_1}$');
+  plot(freqs, PSD_S2, '-',  'DisplayName', '$\Phi_{\hat{x}_2}$');
+  plot(freqs, PSD_cvx, 'k-', 'DisplayName', '$\Phi_{\hat{x}, CVX}$');
+  plot(freqs, PSD_h2i, 'k--', 'DisplayName', '$\Phi_{\hat{x}, \mathcal{H}_2/\mathcal{H}_\infty}$');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Power Spectral Density');
+  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/psd_compare_cvx_h2i.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:psd_compare_cvx_h2i
+#+CAPTION: Power Spectral Density of the Super Sensor obtained with the mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ synthesis ([[./figs/psd_compare_cvx_h2i.png][png]], [[./figs/psd_compare_cvx_h2i.pdf][pdf]])
+[[file:figs/psd_compare_cvx_h2i.png]]
+
+
+#+begin_src matlab :exports none
+  CPS_S1 = 1/pi*cumtrapz(2*pi*freqs, PSD_S1);
+  CPS_S2 = 1/pi*cumtrapz(2*pi*freqs, PSD_S2);
+  CPS_cvx = 1/pi*cumtrapz(2*pi*freqs, PSD_cvx);
+  CPS_h2i = 1/pi*cumtrapz(2*pi*freqs, PSD_h2i);
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(freqs, CPS_S1, '-', 'DisplayName', sprintf('$\\sigma_{\\hat{x}_1} = %.1e$', sqrt(CPS_S1(end))));
+  plot(freqs, CPS_S2, '-', 'DisplayName', sprintf('$\\sigma_{\\hat{x}_2} = %.1e$', sqrt(CPS_S2(end))));
+  plot(freqs, CPS_cvx, 'k-', 'DisplayName', sprintf('$\\sigma_{\\hat{x}_{CVX}} = %.1e$', sqrt(CPS_cvx(end))));
+  plot(freqs, CPS_h2i, 'k--', 'DisplayName', sprintf('$\\sigma_{\\hat{x}_{\\mathcal{H}_2/\\mathcal{H}_\\infty}} = %.1e$', sqrt(CPS_h2i(end))));
+  set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Cumulative Power Spectrum');
+  hold off;
+  xlim([2e-1, freqs(end)]);
+  ylim([1e-10 1e-5]);
+  legend('location', 'southeast');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/cps_compare_cvx_h2i.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:cps_compare_cvx_h2i
+#+CAPTION: Cumulative Power Spectrum of the Super Sensor obtained with the mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ synthesis ([[./figs/cps_compare_cvx_h2i.png][png]], [[./figs/cps_compare_cvx_h2i.pdf][pdf]])
+[[file:figs/cps_compare_cvx_h2i.png]]
+
+** Obtained Super Sensor's Uncertainty
+The uncertainty on the super sensor's dynamics is shown in Fig. [[]].
+
+#+begin_src matlab :exports none
+  G1 = 1 + w1*ultidyn('Delta',[1 1]);
+  G2 = 1 + w2*ultidyn('Delta',[1 1]);
+
+  % We here compute the maximum and minimum phase of the super sensor
+  Dphiss_cvx = 180/pi*asin(abs(squeeze(freqresp(w1*H1, freqs, 'Hz')))+abs(squeeze(freqresp(w2*H2, freqs, 'Hz'))));
+  Dphiss_cvx(abs(squeeze(freqresp(w1*H1, freqs, 'Hz')))+abs(squeeze(freqresp(w2*H2, freqs, 'Hz'))) > 1) = 190;
+
+  Dphiss_h2i = 180/pi*asin(abs(squeeze(freqresp(w1*H1m, freqs, 'Hz')))+abs(squeeze(freqresp(w2*H2m, freqs, 'Hz'))));
+  Dphiss_h2i(abs(squeeze(freqresp(w1*H1m, freqs, 'Hz')))+abs(squeeze(freqresp(w2*H2m, freqs, 'Hz'))) > 1) = 190;
+
+  % 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 = subplot(2,1,1);
+  hold on;
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, 1 + abs(squeeze(freqresp(w1, freqs, 'Hz'))), '--', 'DisplayName', 'Bounds - S1');
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, max(1 - abs(squeeze(freqresp(w1, freqs, 'Hz'))), 0), '--', 'HandleVisibility', 'off');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs, 1 + abs(squeeze(freqresp(w2, freqs, 'Hz'))), '--', 'DisplayName', 'Bounds - S2');
+  set(gca,'ColorOrderIndex',2);
+  plot(freqs, max(1 - abs(squeeze(freqresp(w2, freqs, 'Hz'))), 0), '--', 'HandleVisibility', 'off');
+  plot(freqs, 1 + abs(squeeze(freqresp(w1*H1+w2*H2, freqs, 'Hz'))), 'k--', 'DisplayName', 'Bounds - CVX');
+  plot(freqs, max(1 - abs(squeeze(freqresp(w1*H1+w2*H2, freqs, 'Hz'))), 0), 'k--', 'HandleVisibility', 'off');
+  plot(freqs, 1 + abs(squeeze(freqresp(w1*H1m+w2*H2m, freqs, 'Hz'))), 'k-', 'DisplayName', 'Bounds - $\mathcal{H}_2/\mathcal{H}_\infty$');
+  plot(freqs, max(1 - abs(squeeze(freqresp(w1*H1m+w2*H2m, freqs, 'Hz'))), 0), 'k-', 'HandleVisibility', 'off');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  legend('location', 'southwest');
+  ylabel('Magnitude');
+  ylim([5e-2, 10]);
+  hold off;
+
+  % Phase
+  ax2 = subplot(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, '--');
+  plot(freqs,  Dphiss_cvx, 'k--');
+  plot(freqs, -Dphiss_cvx, 'k--');
+  plot(freqs,  Dphiss_h2i, 'k-');
+  plot(freqs, -Dphiss_h2i, 'k-');
+  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_compare_cvx_h2i.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:super_sensor_uncertainty_compare_cvx_h2i
+#+CAPTION: Super Sensor Dynamical Uncertainty obtained with the mixed synthesis ([[./figs/super_sensor_uncertainty_compare_cvx_h2i.png][png]], [[./figs/super_sensor_uncertainty_compare_cvx_h2i.pdf][pdf]])
+[[file:figs/super_sensor_uncertainty_compare_cvx_h2i.png]]
+
 * H-Infinity synthesis to ensure both performance and robustness
 :PROPERTIES:
 :header-args:matlab+: :tangle matlab/hinf_syn_perf_robust.m

tikz/index.org

@@ -9,7 +9,7 @@
 #+HTML_HEAD: <script src="../js/bootstrap.min.js"></script>
 #+HTML_HEAD: <script src="../js/jquery.stickytableheaders.min.js"></script>
 #+HTML_HEAD: <script src="../js/readtheorg.js"></script>
-#+PROPERTY: header-args:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{/home/thomas/MEGA/These/Papers/dehaeze19_desig_compl_filte/tikz/}{config.tex}")
+#+PROPERTY: header-args:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{/home/thomas/Cloud/These/Papers/dehaeze19_desig_compl_filte/tikz/}{config.tex}")
 #+PROPERTY: header-args:latex+ :imagemagick t :fit yes
 #+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
 #+PROPERTY: header-args:latex+ :imoutoptions -quality 100