2020-10-26 21:35:47 +01:00
#+TITLE : Complementary Filters Shaping Using $\mathcal{H}_\infty$ Synthesis - Matlab Computation
:DRAWER:
#+HTML_LINK_HOME : ../index.html
2020-11-12 10:44:59 +01:00
#+HTML_LINK_UP : ../index.html
#+HTML_HEAD : <link rel="stylesheet" type="text/css" href="https://research.tdehaeze.xyz/css/style.css"/>
#+HTML_HEAD : <script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.js"></script>
2020-10-26 21:35:47 +01:00
2021-04-28 15:59:35 +02:00
#+LaTeX_CLASS : scrreprt
#+LaTeX_HEADER_EXTRA : \input{/home/thomas/Cloud/org-theme/preamble.tex}
2020-10-26 21:35:47 +01:00
#+PROPERTY : header-args:matlab :session *MATLAB*
#+PROPERTY : header-args:matlab+ :tangle matlab/comp_filters_design.m
#+PROPERTY : header-args:matlab+ :comments org
#+PROPERTY : header-args:matlab+ :exports both
#+PROPERTY : header-args:matlab+ :results none
#+PROPERTY : header-args:matlab+ :eval no-export
#+PROPERTY : header-args:matlab+ :noweb yes
#+PROPERTY : header-args:matlab+ :mkdirp yes
#+PROPERTY : header-args:matlab+ :output-dir figs
2021-04-28 15:59:35 +02:00
#+PROPERTY : header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/tikz/org/}{config.tex}")
#+PROPERTY : header-args:latex+ :imagemagick t :fit yes
#+PROPERTY : header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY : header-args:latex+ :imoutoptions -quality 100
#+PROPERTY : header-args:latex+ :results file raw replace
#+PROPERTY : header-args:latex+ :buffer no
#+PROPERTY : header-args:latex+ :tangle no
#+PROPERTY : header-args:latex+ :eval no-export
#+PROPERTY : header-args:latex+ :exports results
#+PROPERTY : header-args:latex+ :mkdirp yes
#+PROPERTY : header-args:latex+ :output-dir figs
#+PROPERTY : header-args:latex+ :post pdf2svg(file=*this*, ext="png")
2020-10-26 21:35:47 +01:00
:END:
* Introduction :ignore:
In this document, the design of complementary filters is studied.
One use of complementary filter is described below:
#+begin_quote
The basic idea of a complementary filter involves taking two or more sensors, filtering out unreliable frequencies for each sensor, and combining the filtered outputs to get a better estimate throughout the entire bandwidth of the system.
To achieve this, the sensors included in the filter should complement one another by performing better over specific parts of the system bandwidth.
#+end_quote
This document is divided into several sections:
2021-04-30 11:18:08 +02:00
- in section [[#sec:h_inf_synthesis_complementary_filters ]], the $\mathcal{H}_\infty$ synthesis is used for generating two complementary filters
2020-10-26 21:35:47 +01:00
- in section [[sec:three_comp_filters ]], a method using the $\mathcal{H}_\infty$ synthesis is proposed to shape three of more complementary filters
- in section [[sec:comp_filters_ligo ]], the $\mathcal{H}_\infty$ synthesis is used and compared with FIR complementary filters used for LIGO
#+begin_note
Add the Matlab code use to obtain the results presented in the paper are accessible [[file:matlab.zip ][here ]] and presented below.
#+end_note
* H-Infinity synthesis of complementary filters
:PROPERTIES:
:header-args:matlab+: :tangle matlab/h_inf_synthesis_complementary_filters.m
:header-args:matlab+: :comments org :mkdirp yes
2021-04-30 11:18:08 +02:00
:CUSTOM_ID: sec:h_inf_synthesis_complementary_filters
2020-10-26 21:35:47 +01:00
:END:
** Introduction :ignore:
#+begin_note
The Matlab file corresponding to this section is accessible [[file:matlab/h_inf_synthesis_complementary_filters.m ][here ]].
#+end_note
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
2021-04-29 15:03:59 +02:00
<<matlab-dir >>
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
2021-04-29 15:03:59 +02:00
<<matlab-init >>
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab
2021-04-29 15:03:59 +02:00
freqs = logspace(-1, 3, 1000);
#+end_src
#+begin_src matlab :tangle no
addpath('./matlab');
addpath('./matlab/src');
#+end_src
#+begin_src matlab :exec no
addpath('./src');
2020-10-26 21:35:47 +01:00
#+end_src
** Synthesis Architecture
We here synthesize two complementary filters using the $\mathcal{H}_\infty$ synthesis.
The goal is to specify upper bounds on the norms of the two complementary filters $H_1(s)$ and $H_2(s)$ while ensuring their complementary property ($H_1(s) + H_2(s) = 1$).
In order to do so, we use the generalized plant shown on figure [[fig:h_infinity_robst_fusion ]] where $W_1(s)$ and $W_2(s)$ are weighting transfer functions that will be used to shape $H_1(s)$ and $H_2(s)$ respectively.
#+name : fig:h_infinity_robst_fusion
#+caption : $\mathcal{H}_\infty$ synthesis of the complementary filters
2021-04-30 11:18:08 +02:00
[[file:figs/h_infinity_robust_fusion.png ]]
2020-10-26 21:35:47 +01:00
The $\mathcal{H}_\infty$ synthesis applied on this generalized plant will give a transfer function $H_2$ (figure [[fig:h_infinity_robst_fusion ]]) such that the $\mathcal{H}_\infty$ norm of the transfer function from $w$ to $[z_1,\ z_2]$ is less than one:
\[ \left\| \begin{array}{c} (1 - H_2(s)) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_ \infty < 1 \]
Thus, if the above condition is verified, we can define $H_1(s) = 1 - H_2(s)$ and we have that:
\[ \left\| \begin{array}{c} H_1(s) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_ \infty < 1 \]
Which is almost (with an maximum error of $\sqrt{2}$) equivalent to:
\begin{align*}
|H_1(j\omega)| &< \frac{1}{|W_1(j\omega)|}, \quad \forall \omega \\
|H_2(j\omega)| &< \frac{1}{|W_2(j\omega)|}, \quad \forall \omega
\end{align*}
We then see that $W_1(s)$ and $W_2(s)$ can be used to shape both $H_1(s)$ and $H_2(s)$ while ensuring their complementary property by the definition of $H_1(s) = 1 - H_2(s)$.
** Design of Weighting Function
A formula is proposed to help the design of the weighting functions:
\begin{equation}
W(s) = \left( \frac{
\frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_ \infty}\right)^{\frac{2}{n}}}} s + \left(\frac{G_0}{G_c}\right)^{\frac{1}{n}}
}{
\left(\frac{1}{G_\infty}\right)^{\frac{1}{n}} \frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_ \infty}\right)^{\frac{2}{n}}}} s + \left(\frac{1}{G_c}\right)^{\frac{1}{n}}
}\right)^n
\end{equation}
The parameters permits to specify:
- the low frequency gain: $G_0 = lim_ {\omega \to 0} |W(j\omega)|$
- the high frequency gain: $G_\infty = lim_ {\omega \to \infty} |W(j\omega)|$
- the absolute gain at $\omega_0$: $G_c = |W(j\omega_0)|$
- the absolute slope between high and low frequency: $n$
The general shape of a weighting function generated using the formula is shown in figure [[fig:weight_formula ]].
2021-04-28 15:59:35 +02:00
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
n = 3; w0 = 2*pi*10; G0 = 1e-3; G1 = 1e1; Gc = 2;
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
W = (((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;
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(W, freqs, 'Hz'))), 'k-');
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
plot([1e-3 1e0], [G0 G0], 'k--', 'LineWidth', 1)
text(1e0, G0, '$\quad G_0$')
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
plot([1e1 1e3], [G1 G1], 'k--', 'LineWidth', 1)
text(1e1,G1,'$G_{\infty}\quad$','HorizontalAlignment', 'right')
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
plot([w0/2/pi w0/2/pi], [1 2*Gc], 'k--', 'LineWidth', 1)
text(w0/2/pi,1,'$\omega_c$','VerticalAlignment', 'top', 'HorizontalAlignment', 'center')
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
plot([w0/2/pi/2 2*w0/2/pi], [Gc Gc], 'k--', 'LineWidth', 1)
text(w0/2/pi/2, Gc, '$G_c \quad$','HorizontalAlignment', 'right')
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
text(w0/5/pi/2, abs(evalfr(W, j*w0/5)), 'Slope: $n \quad$', 'HorizontalAlignment', 'right')
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
text(w0/2/pi, abs(evalfr(W, j*w0)), '$\bullet$', 'HorizontalAlignment', 'center')
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([5e-4, 20]);
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/weight_formula.pdf', 'width', 'wide', 'height', 'normal');
2021-04-28 15:59:35 +02:00
#+end_src
2020-10-26 21:35:47 +01:00
#+name : fig:weight_formula
2021-04-28 15:59:35 +02:00
#+caption : Gain of the Weighting Function formula
#+RESULTS :
[[file:figs/weight_formula.png ]]
2020-10-26 21:35:47 +01:00
#+begin_src matlab
2021-04-29 15:03:59 +02:00
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = (((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;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
2021-05-03 17:47:19 +02:00
ylim([1e-4, 20]);
2021-04-29 15:03:59 +02:00
xticks([0.1, 1, 10, 100, 1000]);
2021-05-03 17:47:19 +02:00
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/weights_W1_W2.pdf', 'width', 'wide', 'height', 'normal');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:weights_W1_W2
#+caption : Weights on the complementary filters $W_1$ and $W_2$ and the associated performance weights
#+RESULTS :
[[file:figs/weights_W1_W2.png ]]
** H-Infinity Synthesis
We define the generalized plant $P$ on matlab.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P = [W1 -W1;
0 W2;
1 0];
2020-10-26 21:35:47 +01:00
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
#+begin_src matlab :results output replace :exports both
2021-04-29 15:03:59 +02:00
[H2, ~, gamma, ~ ] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
2020-10-26 21:35:47 +01:00
#+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.1000 < gamma <= 1050.0000
gamma hamx_eig xinf_eig hamy_eig yinf_eig nrho_xy p/f
1.050e+03 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
525.050 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
262.575 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
131.337 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
65.719 2.8e+01 2.4e-07 4.1e+00 -9.5e-14 0.0000 p
32.909 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
16.505 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
8.302 2.8e+01 2.4e-07 4.1e+00 -7.2e-14 0.0000 p
4.201 2.8e+01 2.4e-07 4.1e+00 -2.5e-25 0.0000 p
2.151 2.7e+01 2.4e-07 4.1e+00 -3.8e-14 0.0000 p
1.125 2.6e+01 2.4e-07 4.1e+00 -5.4e-24 0.0000 p
0.613 2.3e+01 -3.7e+01# 4.1e+00 0.0e+00 0.0000 f
0.869 2.6e+01 -3.7e+02# 4.1e+00 0.0e+00 0.0000 f
0.997 2.6e+01 -1.1e+04# 4.1e+00 0.0e+00 0.0000 f
1.061 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.029 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.013 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.005 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.001 2.6e+01 -3.1e+04# 4.1e+00 -3.8e-14 0.0000 f
1.003 2.6e+01 -2.8e+05# 4.1e+00 0.0e+00 0.0000 f
1.004 2.6e+01 2.4e-07 4.1e+00 -5.8e-24 0.0000 p
1.004 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
Gamma value achieved: 1.0036
#+end_example
We then define the high pass filter $H_1 = 1 - H_2$. The bode plot of both $H_1$ and $H_2$ is shown on figure [[fig:hinf_filters_results ]].
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H1 = 1 - H2;
2020-10-26 21:35:47 +01:00
#+end_src
** Obtained Complementary Filters
The obtained complementary filters are shown on figure [[fig:hinf_filters_results ]].
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 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$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([1e-4, 20]);
yticks([1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1]);
2021-05-03 17:47:19 +02:00
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 18;
2021-04-29 15:03:59 +02:00
% Phase
ax2 = nexttile;
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([-180:90:180]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/hinf_filters_results.pdf', 'width', 'wide', 'height', 'tall');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:hinf_filters_results
#+caption : Obtained complementary filters using $\mathcal{H}_\infty$ synthesis
#+RESULTS :
[[file:figs/hinf_filters_results.png ]]
* Generating 3 complementary filters
:PROPERTIES:
:header-args:matlab+: :tangle matlab/three_comp_filters.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:three_comp_filters >>
** Introduction :ignore:
#+begin_note
The Matlab file corresponding to this section is accessible [[file:matlab/three_comp_filters.m ][here ]].
#+end_note
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
2021-04-29 15:03:59 +02:00
<<matlab-dir >>
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
2021-04-29 15:03:59 +02:00
<<matlab-init >>
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab
2021-04-29 15:03:59 +02:00
freqs = logspace(-2, 4, 1000);
#+end_src
#+begin_src matlab :tangle no
addpath('./matlab');
addpath('./matlab/src');
#+end_src
#+begin_src matlab :exec no
addpath('./src');
2020-10-26 21:35:47 +01:00
#+end_src
** Theory
We want:
\begin{align*}
& |H_1(j\omega)| < 1/|W_1(j\omega)|, \quad \forall\omega\\
& |H_2(j\omega)| < 1/|W_2(j\omega)|, \quad \forall\omega\\
& |H_3(j\omega)| < 1/|W_3(j\omega)|, \quad \forall\omega\\
& H_1(s) + H_2(s) + H_3(s) = 1
\end{align*}
For that, we use the $\mathcal{H}_\infty$ synthesis with the architecture shown on figure [[fig:comp_filter_three_hinf ]].
#+name : fig:comp_filter_three_hinf
#+caption : Generalized architecture for generating 3 complementary filters
2021-04-30 11:18:08 +02:00
[[file:figs/comp_filter_three_hinf.png ]]
2020-10-26 21:35:47 +01:00
The $\mathcal{H}_\infty$ objective is:
\begin{align*}
& |(1 - H_2(j\omega) - H_3(j\omega)) W_1(j\omega)| < 1, \quad \forall\omega\\
& |H_2(j\omega) W_2(j\omega)| < 1, \quad \forall\omega\\
& |H_3(j\omega) W_3(j\omega)| < 1, \quad \forall\omega\\
\end{align*}
And thus if we choose $H_1 = 1 - H_2 - H_3$ we have solved the problem.
** Weights
First we define the weights.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
n = 2; w0 = 2*pi*1; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
W2 = 0.22*(1 + s/2/pi/1)^2/(sqrt(1e-4) + s/2/pi/1)^2* (1 + s/2/pi/10)^2/(1 + s/2/pi/1000)^2;
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W3 = (((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;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca,'ColorOrderIndex',3)
plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
2021-05-03 17:47:19 +02:00
xlim([freqs(1), freqs(end)]); ylim([2e-4, 1.3e1])
leg = legend('location', 'northeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/three_weighting_functions.pdf', 'width', 'wide', 'height', 'normal');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:three_weighting_functions
#+caption : Three weighting functions used for the $\mathcal{H}_\infty$ synthesis of the complementary filters
#+RESULTS :
[[file:figs/three_weighting_functions.png ]]
** H-Infinity Synthesis
Then we create the generalized plant =P= .
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P = [W1 -W1 -W1;
0 W2 0 ;
0 0 W3;
1 0 0];
2020-10-26 21:35:47 +01:00
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis.
#+begin_src matlab :results output replace :exports both
2021-04-29 15:03:59 +02:00
[H, ~, gamma, ~ ] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
2020-10-26 21:35:47 +01:00
#+end_src
#+RESULTS :
#+begin_example
[H, ~, gamma, ~ ] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
Resetting value of Gamma min based on D_11, D_12, D_21 terms
Test bounds: 0.1000 < gamma <= 1050.0000
gamma hamx_eig xinf_eig hamy_eig yinf_eig nrho_xy p/f
1.050e+03 3.2e+00 4.5e-13 6.3e-02 -1.2e-11 0.0000 p
525.050 3.2e+00 1.3e-13 6.3e-02 0.0e+00 0.0000 p
262.575 3.2e+00 2.1e-12 6.3e-02 -1.5e-13 0.0000 p
131.337 3.2e+00 1.1e-12 6.3e-02 -7.2e-29 0.0000 p
65.719 3.2e+00 2.0e-12 6.3e-02 0.0e+00 0.0000 p
32.909 3.2e+00 7.4e-13 6.3e-02 -5.9e-13 0.0000 p
16.505 3.2e+00 1.4e-12 6.3e-02 0.0e+00 0.0000 p
8.302 3.2e+00 1.6e-12 6.3e-02 0.0e+00 0.0000 p
4.201 3.2e+00 1.6e-12 6.3e-02 0.0e+00 0.0000 p
2.151 3.2e+00 1.6e-12 6.3e-02 0.0e+00 0.0000 p
1.125 3.2e+00 2.8e-12 6.3e-02 0.0e+00 0.0000 p
0.613 3.0e+00 -2.5e+03# 6.3e-02 0.0e+00 0.0000 f
0.869 3.1e+00 -2.9e+01# 6.3e-02 0.0e+00 0.0000 f
0.997 3.2e+00 1.9e-12 6.3e-02 0.0e+00 0.0000 p
0.933 3.1e+00 -6.9e+02# 6.3e-02 0.0e+00 0.0000 f
0.965 3.1e+00 -3.0e+03# 6.3e-02 0.0e+00 0.0000 f
0.981 3.1e+00 -8.6e+03# 6.3e-02 0.0e+00 0.0000 f
0.989 3.2e+00 -2.7e+04# 6.3e-02 0.0e+00 0.0000 f
0.993 3.2e+00 -5.7e+05# 6.3e-02 0.0e+00 0.0000 f
0.995 3.2e+00 2.2e-12 6.3e-02 0.0e+00 0.0000 p
0.994 3.2e+00 1.6e-12 6.3e-02 0.0e+00 0.0000 p
0.994 3.2e+00 1.0e-12 6.3e-02 0.0e+00 0.0000 p
Gamma value achieved: 0.9936
#+end_example
** Obtained Complementary Filters
The obtained filters are:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H2 = tf(H(1));
H3 = tf(H(2));
H1 = 1 - H2 - H3;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca,'ColorOrderIndex',3)
plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
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$');
set(gca,'ColorOrderIndex',3)
plot(freqs, abs(squeeze(freqresp(H3, freqs, 'Hz'))), '-', 'DisplayName', '$H_3$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([1e-4, 20]);
2021-05-03 17:47:19 +02:00
leg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 18;
2021-04-29 15:03:59 +02:00
% Phase
ax2 = nexttile;
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'))));
set(gca,'ColorOrderIndex',3)
plot(freqs, 180/pi*phase(squeeze(freqresp(H3, freqs, 'Hz'))));
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-360:90:360]); ylim([-270, 270]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/three_complementary_filters_results.pdf', 'width', 'wide', 'height', 'tall');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:three_complementary_filters_results
#+caption : The three complementary filters obtained after $\mathcal{H}_\infty$ synthesis
#+RESULTS :
[[file:figs/three_complementary_filters_results.png ]]
* Implement complementary filters for LIGO
:PROPERTIES:
:header-args:matlab+: :tangle matlab/comp_filters_ligo.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:comp_filters_ligo >>
** Introduction :ignore:
#+begin_note
The Matlab file corresponding to this section is accessible [[file:matlab/comp_filters_ligo.m ][here ]].
#+end_note
Let's try to design complementary filters that are corresponding to the complementary filters design for the LIGO and described in cite:hua05_low_ligo.
The FIR complementary filters designed in cite:hua05_low_ligo are of order 512.
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
2021-04-29 15:03:59 +02:00
<<matlab-dir >>
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
2021-04-29 15:03:59 +02:00
<<matlab-init >>
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab
2021-04-29 15:03:59 +02:00
freqs = logspace(-3, 0, 1000);
#+end_src
#+begin_src matlab :tangle no
addpath('./matlab');
addpath('./matlab/src');
#+end_src
#+begin_src matlab :exec no
addpath('./src');
2020-10-26 21:35:47 +01:00
#+end_src
** Specifications
The specifications for the filters are:
1. From $0$ to $0.008\text{ Hz}$,the magnitude of the filter’
2. From $0.008\text{ Hz}$ to $0.04\text{ Hz}$, it attenuates the input signal proportional to frequency cubed
3. Between $0.04\text{ Hz}$ and $0.1\text{ Hz}$, the magnitude of the transfer function should be less than 3
4. Above $0.1\text{ Hz}$, the maximum of the magnitude of the complement filter should be as close to zero as possible. In our system, we would like to have the magnitude of the complementary filter to be less than $0.1$. As the filters obtained in cite:hua05_low_ligo have a magnitude of $0.045$, we will set that as our requirement
The specifications are translated in upper bounds of the complementary filters are shown on figure [[fig:ligo_specifications ]].
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$');
set(gca,'ColorOrderIndex',1)
plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',1)
plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',2)
plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-4, 10]);
2021-05-03 17:47:19 +02:00
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/ligo_specifications.pdf', 'width', 'wide', 'height', 'normal');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:ligo_specifications
#+caption : Specification for the LIGO complementary filters
#+RESULTS :
[[file:figs/ligo_specifications.png ]]
2021-04-30 11:18:08 +02:00
** FIR Filter
2020-10-26 21:35:47 +01:00
We here try to implement the FIR complementary filter synthesis as explained in cite:hua05_low_ligo.
For that, we use the [[http://cvxr.com/cvx/ ][CVX matlab Toolbox ]].
We setup the CVX toolbox and use the =SeDuMi= solver.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
cvx_startup;
cvx_solver sedumi;
2020-10-26 21:35:47 +01:00
#+end_src
We define the frequency vectors on which we will constrain the norm of the FIR filter.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
w1 = 0:4.06e-4:0.008;
w2 = 0.008:4.06e-4:0.04;
w3 = 0.04:8.12e-4:0.1;
w4 = 0.1:8.12e-4:0.83;
2020-10-26 21:35:47 +01:00
#+end_src
We then define the order of the FIR filter.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
n = 512;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab
2021-04-29 15:03:59 +02:00
A1 = [ones(length(w1),1), cos(kron(w1'.*(2*pi),[1:n-1]))];
A2 = [ones(length(w2),1), cos(kron(w2'.*(2*pi),[1:n-1]))];
A3 = [ones(length(w3),1), cos(kron(w3'.*(2*pi),[1:n-1]))];
A4 = [ones(length(w4),1), cos(kron(w4'.*(2*pi),[1:n-1]))];
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
B1 = [zeros(length(w1),1), sin(kron(w1'.*(2*pi),[1:n-1]))];
B2 = [zeros(length(w2),1), sin(kron(w2'.*(2*pi),[1:n-1]))];
B3 = [zeros(length(w3),1), sin(kron(w3'.*(2*pi),[1:n-1]))];
B4 = [zeros(length(w4),1), sin(kron(w4'.*(2*pi),[1:n-1]))];
2020-10-26 21:35:47 +01:00
#+end_src
We run the convex optimization.
#+begin_src matlab :results output replace :wrap example
2021-04-29 15:03:59 +02:00
cvx_begin
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
variable y(n+1,1)
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
% t
maximize(-y(1))
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
for i = 1:length(w1)
norm([0 A1(i,:); 0 B1(i,:)]*y) <= 8e-3;
end
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
for i = 1:length(w2)
norm([0 A2(i,:); 0 B2(i,:)]*y) <= 8e-3* (2*pi*w2(i)/(0.008*2*pi))^3;
end
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
for i = 1:length(w3)
norm([0 A3(i,:); 0 B3(i,:)]*y) <= 3;
end
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
for i = 1:length(w4)
norm([[1 0]'- [0 A4(i,:); 0 B4(i,:)]*y]) <= y(1);
end
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
cvx_end
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
h = y(2:end);
2020-10-26 21:35:47 +01:00
#+end_src
#+RESULTS :
#+begin_example
cvx_begin
variable y(n+1,1)
% t
maximize(-y(1))
for i = 1:length(w1)
norm([0 A1(i,:); 0 B1(i,:)]*y) <= 8e-3;
end
for i = 1:length(w2)
norm([0 A2(i,:); 0 B2(i,:)]*y) <= 8e-3* (2*pi*w2(i)/(0.008*2*pi))^3;
end
for i = 1:length(w3)
norm([0 A3(i,:); 0 B3(i,:)]*y) <= 3;
end
for i = 1:length(w4)
norm([[1 0]'- [0 A4(i,:); 0 B4(i,:)]*y]) <= y(1);
end
cvx_end
Calling SeDuMi 1.34: 4291 variables, 1586 equality constraints
For improved efficiency, SeDuMi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.34 (beta) by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 1586, order n = 3220, dim = 4292, blocks = 1073
nnz(A) = 1100727 + 0, nnz(ADA) = 1364794, nnz(L) = 683190
it : b*y gap delta rate t/tP* t/tD* feas cg cg prec
0 : 4.11E+02 0.000
1 : -2.58E+00 1.25E+02 0.000 0.3049 0.9000 0.9000 4.87 1 1 3.0E+02
2 : -2.36E+00 3.90E+01 0.000 0.3118 0.9000 0.9000 1.83 1 1 6.6E+01
3 : -1.69E+00 1.31E+01 0.000 0.3354 0.9000 0.9000 1.76 1 1 1.5E+01
4 : -8.60E-01 7.10E+00 0.000 0.5424 0.9000 0.9000 2.48 1 1 4.8E+00
5 : -4.91E-01 5.44E+00 0.000 0.7661 0.9000 0.9000 3.12 1 1 2.5E+00
6 : -2.96E-01 3.88E+00 0.000 0.7140 0.9000 0.9000 2.62 1 1 1.4E+00
7 : -1.98E-01 2.82E+00 0.000 0.7271 0.9000 0.9000 2.14 1 1 8.5E-01
8 : -1.39E-01 2.00E+00 0.000 0.7092 0.9000 0.9000 1.78 1 1 5.4E-01
9 : -9.99E-02 1.30E+00 0.000 0.6494 0.9000 0.9000 1.51 1 1 3.3E-01
10 : -7.57E-02 8.03E-01 0.000 0.6175 0.9000 0.9000 1.31 1 1 2.0E-01
11 : -5.99E-02 4.22E-01 0.000 0.5257 0.9000 0.9000 1.17 1 1 1.0E-01
12 : -5.28E-02 2.45E-01 0.000 0.5808 0.9000 0.9000 1.08 1 1 5.9E-02
13 : -4.82E-02 1.28E-01 0.000 0.5218 0.9000 0.9000 1.05 1 1 3.1E-02
14 : -4.56E-02 5.65E-02 0.000 0.4417 0.9045 0.9000 1.02 1 1 1.4E-02
15 : -4.43E-02 2.41E-02 0.000 0.4265 0.9004 0.9000 1.01 1 1 6.0E-03
16 : -4.37E-02 8.90E-03 0.000 0.3690 0.9070 0.9000 1.00 1 1 2.3E-03
17 : -4.35E-02 3.24E-03 0.000 0.3641 0.9164 0.9000 1.00 1 1 9.5E-04
18 : -4.34E-02 1.55E-03 0.000 0.4788 0.9086 0.9000 1.00 1 1 4.7E-04
19 : -4.34E-02 8.77E-04 0.000 0.5653 0.9169 0.9000 1.00 1 1 2.8E-04
20 : -4.34E-02 5.05E-04 0.000 0.5754 0.9034 0.9000 1.00 1 1 1.6E-04
21 : -4.34E-02 2.94E-04 0.000 0.5829 0.9136 0.9000 1.00 1 1 9.9E-05
22 : -4.34E-02 1.63E-04 0.015 0.5548 0.9000 0.0000 1.00 1 1 6.6E-05
23 : -4.33E-02 9.42E-05 0.000 0.5774 0.9053 0.9000 1.00 1 1 3.9E-05
24 : -4.33E-02 6.27E-05 0.000 0.6658 0.9148 0.9000 1.00 1 1 2.6E-05
25 : -4.33E-02 3.75E-05 0.000 0.5972 0.9187 0.9000 1.00 1 1 1.6E-05
26 : -4.33E-02 1.89E-05 0.000 0.5041 0.9117 0.9000 1.00 1 1 8.6E-06
27 : -4.33E-02 9.72E-06 0.000 0.5149 0.9050 0.9000 1.00 1 1 4.5E-06
28 : -4.33E-02 2.94E-06 0.000 0.3021 0.9194 0.9000 1.00 1 1 1.5E-06
29 : -4.33E-02 9.73E-07 0.000 0.3312 0.9189 0.9000 1.00 2 2 5.3E-07
30 : -4.33E-02 2.82E-07 0.000 0.2895 0.9063 0.9000 1.00 2 2 1.6E-07
31 : -4.33E-02 8.05E-08 0.000 0.2859 0.9049 0.9000 1.00 2 2 4.7E-08
32 : -4.33E-02 1.43E-08 0.000 0.1772 0.9059 0.9000 1.00 2 2 8.8E-09
iter seconds digits c*x b*y
32 49.4 6.8 -4.3334083581e-02 -4.3334090214e-02
|Ax-b| = 3.7e-09, [Ay-c]_+ = 1.1E-10, |x|= 1.0e+00, |y|= 2.6e+00
Detailed timing (sec)
Pre IPM Post
3.902E+00 4.576E+01 1.035E-02
Max-norms: ||b||=1, ||c|| = 3,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 4.26267.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -0.0433341
h = y(2:end);
#+end_example
Finally, we compute the filter response over the frequency vector defined and the result is shown on figure [[fig:fir_filter_ligo ]] which is very close to the filters obtain in cite:hua05_low_ligo.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
w = [w1 w2 w3 w4];
H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
plot(w, abs(H), 'k-');
plot(w, abs(1-H), 'k--');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 5]);
% Phase
ax2 = nexttile;
hold on;
2021-05-03 17:47:19 +02:00
plot(w, 180/pi*unwrap(angle(H)), 'k-');
plot(w, 180/pi*unwrap(angle(1-H)), 'k--');
2021-04-29 15:03:59 +02:00
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
2021-05-03 17:47:19 +02:00
yticks([-450:90:180]); ylim([-450, 200]);
2021-04-29 15:03:59 +02:00
linkaxes([ax1,ax2],'x');
xlim([1e-3, 1]);
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-05-03 17:47:19 +02:00
exportFig('figs/fir_filter_ligo.pdf', 'width', 'wide', 'height', 'tall');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:fir_filter_ligo
#+caption : FIR Complementary filters obtain after convex optimization
#+RESULTS :
[[file:figs/fir_filter_ligo.png ]]
** Weights
We design weights that will be used for the $\mathcal{H}_\infty$ synthesis of the complementary filters.
These weights will determine the order of the obtained filters.
Here are the requirements on the filters:
- reasonable order
- to be as close as possible to the specified upper bounds
- stable minimum phase
The bode plot of the weights is shown on figure [[fig:ligo_weights ]].
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
w1 = 2*pi*0.008; x1 = 0.35;
w2 = 2*pi*0.04; x2 = 0.5;
w3 = 2*pi*0.05; x3 = 0.5;
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
% Slope of +3 from w1
wH = 0.008*(s^2/w1^2 + 2*x1/w1*s + 1)*(s/w1 + 1);
% Little bump from w2 to w3
wH = wH*(s^2/w2^2 + 2*x2/w2*s + 1)/(s^2/w3^2 + 2*x3/w3*s + 1);
% No Slope at high frequencies
wH = wH/(s^2/w3^2 + 2*x3/w3*s + 1)/(s/w3 + 1);
% Little bump between w2 and w3
w0 = 2*pi*0.045; xi = 0.1; A = 2; n = 1;
wH = wH*((s^2 + 2*w0*xi*A^(1/n)*s + w0^2)/(s^2 + 2*w0*xi*s + w0^2))^n;
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
wH = 1/wH;
wH = minreal(ss(wH));
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
n = 20; Rp = 1; Wp = 2*pi*0.102;
[b,a] = cheby1(n, Rp, Wp, 'high', 's');
wL = 0.04*tf(a, b);
2020-10-26 21:35:47 +01:00
2021-04-29 15:03:59 +02:00
wL = 1/wL;
wL = minreal(ss(wL));
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(inv(wH), freqs, 'Hz'))), '-', 'DisplayName', '$|w_H|^{-1}$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(inv(wL), freqs, 'Hz'))), '-', 'DisplayName', '$|w_L|^{-1}$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec.');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off');
plot([0.1, 10], [0.045, 0.045], 'k--', 'HandleVisibility', 'off');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-3, 10]);
2021-05-03 17:47:19 +02:00
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/ligo_weights.pdf', 'width', 'wide', 'height', 'normal');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:ligo_weights
#+caption : Weights for the $\mathcal{H}_\infty$ synthesis
#+RESULTS :
[[file:figs/ligo_weights.png ]]
** H-Infinity Synthesis
We define the generalized plant as shown on figure [[fig:h_infinity_robst_fusion ]].
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P = [0 wL;
wH -wH;
1 0];
2020-10-26 21:35:47 +01:00
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
#+begin_src matlab :results output replace :exports both :wrap example
2021-04-29 15:03:59 +02:00
[Hl, ~, gamma, ~ ] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
2020-10-26 21:35:47 +01:00
#+end_src
#+RESULTS :
#+begin_example
[Hl, ~, 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.3276 < gamma <= 1.8063
gamma hamx_eig xinf_eig hamy_eig yinf_eig nrho_xy p/f
1.806 1.4e-02 -1.7e-16 3.6e-03 -4.8e-12 0.0000 p
1.067 1.3e-02 -4.2e-14 3.6e-03 -1.9e-12 0.0000 p
0.697 1.3e-02 -3.0e-01# 3.6e-03 -3.5e-11 0.0000 f
0.882 1.3e-02 -9.5e-01# 3.6e-03 -1.2e-34 0.0000 f
0.975 1.3e-02 -2.7e+00# 3.6e-03 -1.6e-12 0.0000 f
1.021 1.3e-02 -8.7e+00# 3.6e-03 -4.5e-16 0.0000 f
1.044 1.3e-02 -6.5e-14 3.6e-03 -3.0e-15 0.0000 p
1.032 1.3e-02 -1.8e+01# 3.6e-03 0.0e+00 0.0000 f
1.038 1.3e-02 -3.8e+01# 3.6e-03 0.0e+00 0.0000 f
1.041 1.3e-02 -8.3e+01# 3.6e-03 -2.9e-33 0.0000 f
1.042 1.3e-02 -1.9e+02# 3.6e-03 -3.4e-11 0.0000 f
1.043 1.3e-02 -5.3e+02# 3.6e-03 -7.5e-13 0.0000 f
Gamma value achieved: 1.0439
#+end_example
The high pass filter is defined as $H_H = 1 - H_L$.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
Hh = 1 - Hl;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
Hh = minreal(Hh);
Hl = minreal(Hl);
2020-10-26 21:35:47 +01:00
#+end_src
The size of the filters is shown below.
#+begin_src matlab :exports results :results output replace :wrap example
2021-04-29 15:03:59 +02:00
size(Hh), size(Hl)
2020-10-26 21:35:47 +01:00
#+end_src
#+RESULTS :
#+begin_example
size(Hh), size(Hl)
State-space model with 1 outputs, 1 inputs, and 27 states.
State-space model with 1 outputs, 1 inputs, and 27 states.
#+end_example
The bode plot of the obtained filters as shown on figure [[fig:hinf_synthesis_ligo_results ]].
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
hold on;
set(gca,'ColorOrderIndex',1);
plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$');
set(gca,'ColorOrderIndex',1);
plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',1);
plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',2);
plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$');
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$H_H$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'DisplayName', '$H_L$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-3, 10]);
2021-05-03 17:47:19 +02:00
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/hinf_synthesis_ligo_results.pdf', 'width', 'wide', 'height', 'normal');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:hinf_synthesis_ligo_results
#+caption : Obtained complementary filters using the $\mathcal{H}_\infty$ synthesis
#+RESULTS :
[[file:figs/hinf_synthesis_ligo_results.png ]]
2021-04-30 11:18:08 +02:00
** Compare FIR and H-Infinity Filters
2020-10-26 21:35:47 +01:00
Let's now compare the FIR filters designed in cite:hua05_low_ligo and the one obtained with the $\mathcal{H}_ \infty$ synthesis on figure [[fig:comp_fir_ligo_hinf ]].
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$');
set(gca,'ColorOrderIndex',1);
plot(w, abs(H), '--', ...
'DisplayName', '$H_H(s)$ - FIR');
set(gca,'ColorOrderIndex',2);
plot(w, abs(1-H), '--', ...
'DisplayName', '$H_L(s)$ - FIR');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 10]);
2021-04-30 11:18:08 +02:00
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
2021-04-29 15:03:59 +02:00
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1);
2021-04-30 11:18:08 +02:00
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hh, freqs, 'Hz')))), '-');
2021-04-29 15:03:59 +02:00
set(gca,'ColorOrderIndex',2);
2021-04-30 11:18:08 +02:00
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hl, freqs, 'Hz')))), '-');
2021-04-29 15:03:59 +02:00
set(gca,'ColorOrderIndex',1);
2021-04-30 11:18:08 +02:00
plot(w, 180/pi*unwrap(angle(H)), '--');
set(gca,'ColorOrderIndex',2);
plot(w, 180/pi*unwrap(angle(1-H)), '--');
2021-04-29 15:03:59 +02:00
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
2021-04-30 11:18:08 +02:00
yticks([-450:90:180]); ylim([-450, 200]);
2021-04-29 15:03:59 +02:00
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
2020-10-26 21:35:47 +01:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/comp_fir_ligo_hinf.pdf', 'width', 'wide', 'height', 'tall');
2020-10-26 21:35:47 +01:00
#+end_src
#+name : fig:comp_fir_ligo_hinf
#+caption : Comparison between the FIR filters developped for LIGO and the $\mathcal{H}_\infty$ complementary filters
#+RESULTS :
[[file:figs/comp_fir_ligo_hinf.png ]]
2021-04-28 15:59:35 +02:00
* Alternative Synthesis
2021-04-29 15:03:59 +02:00
** Introduction :ignore:
2021-04-28 15:59:35 +02:00
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
2021-04-29 15:03:59 +02:00
<<matlab-dir >>
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
2021-04-29 15:03:59 +02:00
<<matlab-init >>
#+end_src
#+begin_src matlab :tangle no
addpath('./matlab');
addpath('./matlab/src');
#+end_src
2021-05-03 17:47:19 +02:00
#+begin_src matlab :eval no
2021-04-29 15:03:59 +02:00
addpath('./src');
2021-04-28 15:59:35 +02:00
#+end_src
** Two generalized plants
In order to synthesize the complementary filter using the proposed method, we can use two alternative generalized plant as shown in Figures [[fig:h_infinity_arch_1 ]] and [[fig:h_infinity_arch_2 ]].
\begin{equation}
P_1 = \begin{bmatrix} W_1 & -W_1 \\ 0 & W_2 \\ 1 & 0 \end{bmatrix}
\end{equation}
#+begin_src latex :file h_infinity_arch_1.pdf
\begin{tikzpicture}
\node[block={4.5cm}{3.0cm}, fill=black!20!white, dashed] (P) {};
\node[above] at (P.north) {$P_1(s)$};
\coordinate[] (inputw) at ($(P.south west)!0.75!(P.north west) + (-0.7, 0)$);
\coordinate[] (inputu) at ($(P.south west)!0.35!(P.north west) + (-0.7, 0)$);
\coordinate[] (output1) at ($(P.south east)!0.75!(P.north east) + ( 0.7, 0)$);
\coordinate[] (output2) at ($(P.south east)!0.35!(P.north east) + ( 0.7, 0)$);
\coordinate[] (outputv) at ($(P.south east)!0.1!(P.north east) + ( 0.7, 0)$);
\node[block, left=1.4 of output1] (W1){$W_1(s)$};
\node[block, left=1.4 of output2] (W2){$W_2(s)$};
\node[addb={+}{}{}{}{-}, left=of W1] (sub) {};
\node[block, below=0.3 of P] (H2) {$H_2(s)$};
\draw[->] (inputw) node[above right]{$w$} -- (sub.west);
\draw[->] (H2.west) -| ($(inputu)+(0.35, 0)$) node[above]{$u$} -- (W2.west);
\draw[->] (inputu-|sub) node[branch]{} -- (sub.south);
\draw[->] (sub.east) -- (W1.west);
\draw[->] ($(sub.west)+(-0.6, 0)$) node[branch]{} |- ($(outputv)+ (-0.35, 0)$) node[above]{$v$} |- (H2.east);
\draw[->] (W1.east) -- (output1)node[above left]{$z_1$};
\draw[->] (W2.east) -- (output2)node[above left]{$z_2$};
\end{tikzpicture}
#+end_src
#+name : fig:h_infinity_arch_1
#+caption : Complementary Filter Synthesis - Conf 1
#+RESULTS :
[[file:figs/h_infinity_arch_1.png ]]
\begin{equation}
P_2 = \begin{bmatrix} 0 & W_1 & 1 \\ W_2 & -W_1 & 0 \end{bmatrix}
\end{equation}
#+begin_src latex :file h_infinity_arch_2.pdf
2021-04-29 15:03:59 +02:00
\begin{tikzpicture}
\node[block={4.5cm}{4.5cm}, fill=black!20!white, dashed] (P) {};
\node[above] at (P.north) {$P_2(s)$};
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
\coordinate[] (input2) at ($(P.south west)!0.85!(P.north west) + (-0.7, 0)$);
\coordinate[] (input1) at ($(P.south west)!0.55!(P.north west) + (-0.7, 0)$);
\coordinate[] (inputu) at ($(P.south west)!0.3!( P.north west) + (-0.7, 0)$);
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
\coordinate[] (outputz) at ($(P.south east)!0.3!(P.north east) + (0.7, 0)$);
\coordinate[] (outputv) at ($(P.south east)!0.1!(P.north east) + (0.7, 0)$);
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
\node[block, right=1.4 of input2] (W2){$W_2(s)$};
\node[block, right=1.4 of input1] (W1){$W_1(s)$};
\node[addb={+}{-}{}{}{}, right=of W1] (sub) {};
\node[addb, left=2.5 of outputz] (add) {};
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
\node[block, below=0.3 of P] (H2) {$H_2(s)$};
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
\draw[->] (input2) node[above right]{$w_2$} -- (W2.west);
\draw[->] (input1) node[above right]{$w_1$} -- (W1.west);
\draw[->] (W2.east) -| (sub.north);
\draw[->] (W1.east) -- (sub.west);
\draw[->] (W1-|add)node[branch]{} -- (add.north);
\draw[->] (sub.south) |- (outputv) node[above left]{$v$} |- (H2.east);
\draw[->] (H2.west) -| (inputu) node[above right]{$u$} -- (add.west);
\draw[->] (add.east) -- (outputz) node[above left]{$z$};
\end{tikzpicture}
2021-04-28 15:59:35 +02:00
#+end_src
#+name : fig:h_infinity_arch_2
#+caption : Complementary Filter Synthesis - Conf 2
#+RESULTS :
[[file:figs/h_infinity_arch_2.png ]]
Let's run the $\mathcal{H}_\infty$ synthesis for both generalized plant using the same weights and see if the obtained filters are the same:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = (((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;
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P1 = [W1 -W1;
0 W2;
1 0];
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :results output replace :exports both
2021-04-29 15:03:59 +02:00
[H2, ~, gamma, ~ ] = hinfsyn(P1, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
2021-04-28 15:59:35 +02:00
#+end_src
#+RESULTS :
#+begin_example
[H2, ~, gamma, ~ ] = hinfsyn(P1, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
Test bounds: 0.3263 <= gamma <= 1000
gamma X>=0 Y>=0 rho(XY)<1 p/f
1.807e+01 1.4e-07 0.0e+00 1.185e-18 p
2.428e+00 1.5e-07 0.0e+00 1.285e-18 p
8.902e-01 -2.9e+02 # -7.1e-17 5.168e-19 f
1.470e+00 1.5e-07 0.0e+00 1.462e-14 p
1.144e+00 1.5e-07 0.0e+00 1.260e-14 p
1.009e+00 1.5e-07 0.0e+00 4.120e-13 p
9.478e-01 -6.8e+02 # -2.4e-17 1.449e-14 f
9.780e-01 -1.6e+03 # -7.3e-17 6.791e-14 f
9.934e-01 -4.2e+03 # -1.2e-16 3.524e-14 f
1.001e+00 -2.0e+04 # -2.3e-17 5.717e-20 f
1.005e+00 1.5e-07 0.0e+00 8.953e-18 p
1.003e+00 -2.2e+05 # -1.8e-17 3.225e-12 f
1.004e+00 1.5e-07 0.0e+00 2.445e-12 p
Limiting gains...
1.004e+00 1.6e-07 0.0e+00 5.811e-18 p
Best performance (actual): 1.004
#+end_example
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P2 = [0 W1 1;
W2 -W1 0];
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :results output replace :exports both
2021-04-29 15:03:59 +02:00
[H2b, ~, gamma, ~ ] = hinfsyn(P2, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
2021-04-28 15:59:35 +02:00
#+end_src
#+RESULTS :
#+begin_example
[H2b, ~, gamma, ~ ] = hinfsyn(P2, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
Test bounds: 0.3263 <= gamma <= 1000
gamma X>=0 Y>=0 rho(XY)<1 p/f
1.807e+01 0.0e+00 1.4e-07 2.055e-16 p
2.428e+00 0.0e+00 1.4e-07 1.894e-18 p
8.902e-01 -2.1e-16 -2.7e+02 # 1.466e-16 f
1.470e+00 0.0e+00 1.4e-07 4.118e-16 p
1.144e+00 0.0e+00 1.5e-07 2.105e-18 p
1.009e+00 0.0e+00 1.5e-07 2.590e-13 p
9.478e-01 -9.5e-17 -6.3e+02 # 1.663e-19 f
9.780e-01 -1.1e-16 -1.5e+03 # 1.546e-14 f
9.934e-01 -2.8e-17 -4.0e+03 # 3.934e-14 f
1.001e+00 -3.1e-17 -1.9e+04 # 1.191e-19 f
1.005e+00 0.0e+00 1.5e-07 1.443e-12 p
1.003e+00 -8.3e-17 -2.1e+05 # 8.807e-13 f
1.004e+00 0.0e+00 1.5e-07 1.459e-15 p
Limiting gains...
1.004e+00 0.0e+00 1.5e-07 9.086e-19 p
Best performance (actual): 1.004
#+end_example
And indeed, we can see that the exact same filters are obtained (Figure [[fig:hinf_comp_P1_P2_syn ]]).
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
freqs = logspace(-2, 4, 1000);
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
plot(freqs, abs(squeeze(freqresp(H2b, freqs, 'Hz'))), '--', 'DisplayName', '$H_{2b}$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
legend('location', 'southeast', 'FontSize', 8);
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/hinf_comp_P1_P2_syn.pdf', 'width', 'wide', 'height', 'normal');
2021-04-28 15:59:35 +02:00
#+end_src
#+name : fig:hinf_comp_P1_P2_syn
#+caption : Comparison of $H_2(s)$ when using $P_1(s)$ or $P_2(s)$
#+RESULTS :
[[file:figs/hinf_comp_P1_P2_syn.png ]]
** Shaping the Low pass filter or the high pass filter?
Let's see if there is a difference by explicitly shaping $H_1(s)$ or $H_2(s)$.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = (((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;
2021-04-28 15:59:35 +02:00
#+end_src
Let's first synthesize $H_1(s)$:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P1 = [W2 -W2;
0 W1;
1 0];
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :results output replace :exports both
2021-04-29 15:03:59 +02:00
[H1, ~, gamma, ~ ] = hinfsyn(P1, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
2021-04-28 15:59:35 +02:00
#+end_src
#+RESULTS :
#+begin_example
[H1, ~, gamma, ~ ] = hinfsyn(P1, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
Test bounds: 0.3263 <= gamma <= 1.712
gamma X>=0 Y>=0 rho(XY)<1 p/f
7.476e-01 -2.5e+01 # -8.3e-18 4.938e-20 f
1.131e+00 1.9e-07 0.0e+00 1.566e-16 p
9.197e-01 -1.4e+02 # -7.9e-17 4.241e-17 f
1.020e+00 1.9e-07 0.0e+00 2.095e-16 p
9.686e-01 -3.8e+02 # -7.0e-17 1.463e-23 f
9.940e-01 -1.5e+03 # -1.3e-17 3.168e-19 f
1.007e+00 1.9e-07 0.0e+00 1.696e-15 p
1.000e+00 -4.8e+03 # -7.1e-18 7.203e-20 f
1.004e+00 1.9e-07 0.0e+00 1.491e-14 p
1.002e+00 -1.1e+04 # -2.6e-16 2.579e-14 f
1.003e+00 -2.8e+04 # -6.0e-18 8.558e-20 f
Limiting gains...
1.004e+00 2.0e-07 0.0e+00 5.647e-18 p
1.004e+00 1.0e-06 0.0e+00 5.648e-18 p
Best performance (actual): 1.004
#+end_example
And now $H_2(s)$:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P2 = [W1 -W1;
0 W2;
1 0];
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :results output replace :exports both
2021-04-29 15:03:59 +02:00
[H2, ~, gamma, ~ ] = hinfsyn(P2, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
2021-04-28 15:59:35 +02:00
#+end_src
#+RESULTS :
#+begin_example
[H2b, ~, gamma, ~ ] = hinfsyn(P2, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
Test bounds: 0.3263 <= gamma <= 1000
gamma X>=0 Y>=0 rho(XY)<1 p/f
1.807e+01 1.4e-07 0.0e+00 1.185e-18 p
2.428e+00 1.5e-07 0.0e+00 1.285e-18 p
8.902e-01 -2.9e+02 # -7.1e-17 5.168e-19 f
1.470e+00 1.5e-07 0.0e+00 1.462e-14 p
1.144e+00 1.5e-07 0.0e+00 1.260e-14 p
1.009e+00 1.5e-07 0.0e+00 4.120e-13 p
9.478e-01 -6.8e+02 # -2.4e-17 1.449e-14 f
9.780e-01 -1.6e+03 # -7.3e-17 6.791e-14 f
9.934e-01 -4.2e+03 # -1.2e-16 3.524e-14 f
1.001e+00 -2.0e+04 # -2.3e-17 5.717e-20 f
1.005e+00 1.5e-07 0.0e+00 8.953e-18 p
1.003e+00 -2.2e+05 # -1.8e-17 3.225e-12 f
1.004e+00 1.5e-07 0.0e+00 2.445e-12 p
Limiting gains...
1.004e+00 1.6e-07 0.0e+00 5.811e-18 p
Best performance (actual): 1.004
#+end_example
And compare $H_1(s)$ with $1 - H_2(s)$ and $H_2(s)$ with $1 - H_1(s)$ in Figure [[fig:hinf_comp_H1_H2_syn ]].
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
freqs = logspace(-2, 4, 1000);
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
plot(freqs, abs(squeeze(freqresp(1-H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
plot(freqs, abs(squeeze(freqresp(1-H2, freqs, 'Hz'))), '--', 'DisplayName', '$H_2$');
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '--', 'DisplayName', '$H_1$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/hinf_comp_H1_H2_syn.pdf', 'width', 'wide', 'height', 'normal');
2021-04-28 15:59:35 +02:00
#+end_src
#+name : fig:hinf_comp_H1_H2_syn
#+caption : Comparison of $H_1(s)$ with $1-H_2(s)$, and $H_2(s)$ with $1-H_1(s)$
#+RESULTS :
[[file:figs/hinf_comp_H1_H2_syn.png ]]
2021-04-29 15:03:59 +02:00
** Using Feedback architecture
#+begin_src matlab
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = (((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;
#+end_src
Let's first synthesize $H_1(s)$:
#+begin_src matlab
P = [W1 -W1;
0 W2;
1 -1];
#+end_src
#+begin_src matlab :results output replace :exports both
[K, ~, gamma, ~ ] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
#+end_src
#+begin_src matlab
H1 = inv(1 + K);
H2 = 1 - H1;
#+end_src
#+begin_src matlab :exports none
freqs = logspace(-2, 4, 1000);
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
set(gca,'ColorOrderIndex',1);
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
#+end_src
** Adding feature in the filters
#+begin_src matlab
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = (((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;
#+end_src
#+begin_src matlab
Wf = (1 + s/2/pi/1)/s;
Wf = s/(1 + s/2/pi/1e2);
% W2 = W2/Wf/ (1 + s/2/pi/1e3);
#+end_src
#+begin_src matlab
P = [W1 -Wf*W1;
0 Wf*W2;
1 -Wf];
#+end_src
#+begin_src matlab :results output replace :exports both
[Ka, ~, gamma, ~ ] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
#+end_src
#+begin_src matlab
K = Ka*Wf;
#+end_src
#+begin_src matlab
H1 = inv(1 + K);
H2 = 1 - H1;
#+end_src
#+begin_src matlab :exports none
freqs = logspace(-2, 4, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
set(gca,'ColorOrderIndex',1);
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
% Phase
ax2 = nexttile;
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([-180:90:180]);
ylim([-180, 180]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
#+end_src
2021-05-03 17:47:19 +02:00
** Compare "open-loop" shaping with "close-loop" shaping
*** Simple weights
#+begin_src matlab
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = ((((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);
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = ((((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);
#+end_src
#+begin_src matlab
Pol = [0 W1 1;
W2 -W1 0];
#+end_src
#+begin_src matlab :results value replace :exports results :tangle no
sprintf('The number of states of Pol is %i', length(Pol.StateName))
#+end_src
#+RESULTS :
: The number of states of Pol is 7
#+begin_src matlab :results output replace :exports both
tic;
[Hol, ~, gamma, ~ ] = hinfsyn(Pol, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
toc;
#+end_src
#+begin_src matlab
Hol_2 = Hol;
Hol_1 = 1 - Hol;
#+end_src
#+begin_src matlab
Pcl = [0 W2 1
W1 -W2 -1];
#+end_src
#+begin_src matlab :results value replace :exports results :tangle no
sprintf('The number of states of Pcl is %i', length(Pcl.StateName))
#+end_src
#+RESULTS :
: The number of states of Pcl is 8
#+begin_src matlab :results output replace :exports both
tic;
[Hcl, ~, gamma, ~ ] = hinfsyn(Pcl, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
toc;
#+end_src
#+begin_src matlab
Hcl_1 = 1 - 1/(1 + Hcl);
Hcl_2 = 1/(1 + Hcl);
#+end_src
#+begin_src matlab :results output replace :exports results
size(Hol_1)
size(Hol_2)
size(Hcl_1)
size(Hcl_2)
#+end_src
#+RESULTS :
#+begin_example
size(Hol_1)
State-space model with 1 outputs, 1 inputs, and 5 states.
size(Hol_2)
State-space model with 1 outputs, 1 inputs, and 5 states.
size(Hcl_1)
State-space model with 1 outputs, 1 inputs, and 5 states.
size(Hcl_2)
State-space model with 1 outputs, 1 inputs, and 5 states.
'org_babel_eoe'
ans =
'org_babel_eoe'
#+end_example
#+begin_src matlab :exports none
freqs = logspace(-1, 3, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hol_1, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$ OL');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hol_2, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$ OL');
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hcl_1, freqs, 'Hz'))), '--', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$ CL');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hcl_2, freqs, 'Hz'))), '--', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$ CL');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hol_1, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hol_2, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hcl_1, freqs, 'Hz')))), '--');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hcl_2, freqs, 'Hz')))), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-450:90:180]); ylim([-450, 200]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
#+end_src
*** Simple weights with LMI
#+begin_src matlab
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = ((((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);
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = ((((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);
#+end_src
#+begin_src matlab
Pol = [0 W1 1;
W2 -W1 0];
#+end_src
#+begin_src matlab :results value replace :exports results :tangle no
sprintf('The number of states of Pol is %i', length(Pol.StateName))
#+end_src
#+RESULTS :
: The number of states of Pol is 7
#+begin_src matlab :results output replace :exports both
tic;
[Hol, ~, gamma, ~ ] = hinfsyn(Pol, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
toc;
#+end_src
#+begin_src matlab
Hol_2 = Hol;
Hol_1 = 1 - Hol;
#+end_src
#+begin_src matlab
Pcl = [0 W2 1
W1 -W2 -1];
#+end_src
#+begin_src matlab :results value replace :exports results :tangle no
sprintf('The number of states of Pcl is %i', length(Pcl.StateName))
#+end_src
#+RESULTS :
: The number of states of Pcl is 8
#+begin_src matlab :results output replace :exports both
tic;
[Hcl, ~, gamma, ~ ] = hinfsyn(Pcl, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
toc;
#+end_src
#+begin_src matlab
Hcl_1 = 1 - 1/(1 + Hcl);
Hcl_2 = 1/(1 + Hcl);
#+end_src
#+begin_src matlab :results output replace :exports results
size(Hol_1)
size(Hol_2)
size(Hcl_1)
size(Hcl_2)
#+end_src
#+RESULTS :
#+begin_example
size(Hol_1)
State-space model with 1 outputs, 1 inputs, and 5 states.
size(Hol_2)
State-space model with 1 outputs, 1 inputs, and 5 states.
size(Hcl_1)
State-space model with 1 outputs, 1 inputs, and 5 states.
size(Hcl_2)
State-space model with 1 outputs, 1 inputs, and 5 states.
'org_babel_eoe'
ans =
'org_babel_eoe'
#+end_example
#+begin_src matlab :exports none
freqs = logspace(-1, 3, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hol_1, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$ OL');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hol_2, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$ OL');
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hcl_1, freqs, 'Hz'))), '--', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$ CL');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hcl_2, freqs, 'Hz'))), '--', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$ CL');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hol_1, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hol_2, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hcl_1, freqs, 'Hz')))), '--');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hcl_2, freqs, 'Hz')))), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-450:90:180]); ylim([-450, 200]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
#+end_src
*** Complex weights
#+begin_src matlab :exports none
w1 = 2*pi*0.008; x1 = 0.35;
w2 = 2*pi*0.04; x2 = 0.5;
w3 = 2*pi*0.05; x3 = 0.5;
% Slope of +3 from w1
wH = 0.008*(s^2/w1^2 + 2*x1/w1*s + 1)*(s/w1 + 1);
% Little bump from w2 to w3
wH = wH*(s^2/w2^2 + 2*x2/w2*s + 1)/(s^2/w3^2 + 2*x3/w3*s + 1);
% No Slope at high frequencies
wH = wH/(s^2/w3^2 + 2*x3/w3*s + 1)/(s/w3 + 1);
% Little bump between w2 and w3
w0 = 2*pi*0.045; xi = 0.1; A = 2; n = 1;
wH = wH*((s^2 + 2*w0*xi*A^(1/n)*s + w0^2)/(s^2 + 2*w0*xi*s + w0^2))^n;
wH = 1/wH;
W1 = wH;
#+end_src
#+begin_src matlab :exports none
n = 20; Rp = 1; Wp = 2*pi*0.102;
[b,a] = cheby1(n, Rp, Wp, 'high', 's');
wL = 0.04*tf(a, b);
wL = 1/wL;
W2 = wL;
#+end_src
#+begin_src matlab
Pol = ss([0 W1 1;
W2 -W1 0]);
#+end_src
#+begin_src matlab
Pcl = ss([0 W2 1
W1 -W2 -1]);
#+end_src
#+begin_src matlab :results output replace :exports results
size(Pol)
size(Pcl)
#+end_src
#+RESULTS :
: size(Pol)
: State-space model with 2 outputs, 3 inputs, and 27 states.
: size(Pcl)
: State-space model with 2 outputs, 3 inputs, and 27 states.
#+begin_src matlab :results output replace :exports both
tic;
for iter = 1:10
[Hol, ~, gamma, ~ ] = hinfsyn(Pol, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'off');
end;
toc;
#+end_src
#+begin_src matlab
Hol_1 = 1 - Hol;
Hol_2 = Hol;
#+end_src
#+begin_src matlab :results output replace :exports both
tic;
for iter = 1:10
[Hcl, ~, gamma, ~ ] = hinfsyn(Pcl, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'off');
end
toc;
#+end_src
#+begin_src matlab
Hcl_1 = 1 - 1/(1 + Hcl);
Hcl_2 = 1/(1 + Hcl);
#+end_src
#+begin_src matlab :results output replace :exports results
size(Hol_1)
size(Hol_2)
size(Hcl_1)
size(Hcl_2)
#+end_src
#+RESULTS :
#+begin_example
size(Hol_1)
State-space model with 1 outputs, 1 inputs, and 27 states.
size(Hol_2)
State-space model with 1 outputs, 1 inputs, and 27 states.
size(Hcl_1)
State-space model with 1 outputs, 1 inputs, and 27 states.
size(Hcl_2)
State-space model with 1 outputs, 1 inputs, and 27 states.
'org_babel_eoe'
ans =
'org_babel_eoe'
#+end_example
#+begin_src matlab :exports none
freqs = logspace(-3, 1, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hol_1, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$ OL');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hol_2, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$ OL');
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hcl_1, freqs, 'Hz'))), '--', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$ CL');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hcl_2, freqs, 'Hz'))), '--', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$ CL');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hol_1, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hol_2, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hcl_1, freqs, 'Hz')))), '--');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hcl_2, freqs, 'Hz')))), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-450:90:180]); ylim([-450, 200]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
#+end_src
*** Complex weights with alternative conf
#+begin_src matlab
Pcl = ss([W1 -W1;
0 W2;
1 -1]);
#+end_src
#+begin_src matlab :results output replace :exports both
tic;
[Hcl, ~, gamma, ~ ] = hinfsyn(Pcl, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
toc;
#+end_src
#+begin_src matlab
Hcl_1 = 1 - 1/(1 + Hcl);
Hcl_2 = 1/(1 + Hcl);
#+end_src
*** Complex weights with LMI
#+begin_src matlab :exports none
w1 = 2*pi*0.008; x1 = 0.35;
w2 = 2*pi*0.04; x2 = 0.5;
w3 = 2*pi*0.05; x3 = 0.5;
% Slope of +3 from w1
wH = 0.008*(s^2/w1^2 + 2*x1/w1*s + 1)*(s/w1 + 1);
% Little bump from w2 to w3
wH = wH*(s^2/w2^2 + 2*x2/w2*s + 1)/(s^2/w3^2 + 2*x3/w3*s + 1);
% No Slope at high frequencies
wH = wH/(s^2/w3^2 + 2*x3/w3*s + 1)/(s/w3 + 1);
% Little bump between w2 and w3
w0 = 2*pi*0.045; xi = 0.1; A = 2; n = 1;
wH = wH*((s^2 + 2*w0*xi*A^(1/n)*s + w0^2)/(s^2 + 2*w0*xi*s + w0^2))^n;
wH = 1/wH;
W1 = ss(minreal(wH));
#+end_src
#+begin_src matlab :exports none
n = 20; Rp = 1; Wp = 2*pi*0.102;
[b,a] = cheby1(n, Rp, Wp, 'high', 's');
wL = 0.04*tf(a, b);
wL = 1/wL;
W2 = ss(minreal(wL));
#+end_src
#+begin_src matlab
Pol = [0 W1 1;
W2 -W1 0];
#+end_src
#+begin_src matlab
Pcl = [0 W2 1
W1 -W2 -1];
#+end_src
#+begin_src matlab :results output replace :exports both
tic;
[Hol, ~, gamma, ~ ] = hinfsyn(Pol, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
toc;
#+end_src
#+begin_src matlab
Hol_1 = 1 - Hol;
Hol_2 = Hol;
#+end_src
#+begin_src matlab :results output replace :exports both
tic;
[Hcl, ~, gamma, ~ ] = hinfsyn(Pcl, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
toc;
#+end_src
#+begin_src matlab
Hcl_1 = 1 - 1/(1 + Hcl);
Hcl_2 = 1/(1 + Hcl);
#+end_src
#+begin_src matlab :results output replace :exports results
size(Hol_1)
size(Hol_2)
size(Hcl_1)
size(Hcl_2)
#+end_src
#+RESULTS :
#+begin_example
size(Hol_1)
State-space model with 1 outputs, 1 inputs, and 34 states.
size(Hol_2)
State-space model with 1 outputs, 1 inputs, and 34 states.
size(Hcl_1)
State-space model with 1 outputs, 1 inputs, and 47 states.
size(Hcl_2)
State-space model with 1 outputs, 1 inputs, and 47 states.
'org_babel_eoe'
ans =
'org_babel_eoe'
#+end_example
#+begin_src matlab :exports none
freqs = logspace(-3, 1, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hol_1, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$ OL');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hol_2, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$ OL');
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hcl_1, freqs, 'Hz'))), '--', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_ \infty$ CL');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hcl_2, freqs, 'Hz'))), '--', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_ \infty$ CL');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hol_1, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hol_2, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hcl_1, freqs, 'Hz')))), '--');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hcl_2, freqs, 'Hz')))), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-450:90:180]); ylim([-450, 200]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
#+end_src
*** Conclusion
#+begin_important
There is no difference between " open-loop" shaping and "close-loop" shaping:
- same " solving" time
- same obtained filter orders
#+end_important
2021-04-28 15:59:35 +02:00
* Impose a positive slope at DC or a negative slope at infinite frequency
2021-04-29 15:03:59 +02:00
** Introduction :ignore:
2021-04-28 15:59:35 +02:00
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
2021-04-29 15:03:59 +02:00
<<matlab-dir >>
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
2021-04-29 15:03:59 +02:00
<<matlab-init >>
#+end_src
#+begin_src matlab :tangle no
addpath('./matlab');
addpath('./matlab/src');
#+end_src
#+begin_src matlab :exec no
addpath('./src');
2021-04-28 15:59:35 +02:00
#+end_src
** Manually shift zeros to the origin after synthesis
Suppose we want $H_2(s)$ to be an high pass filter with a slope of +2 at low frequency (from 0Hz).
We cannot impose that using the weight $W_2(s)$ as it would be improper.
However, we may manually shift 2 of the low frequency zeros to the origin.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
n = 3; w0 = 2*pi*10; G0 = 1e4; G1 = 0.1; Gc = 1/2;
W2 = (((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;
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P = [W1 -W1;
0 W2;
1 0];
2021-04-28 15:59:35 +02:00
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
[H2, ~, gamma, ~ ] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab
2021-04-29 15:03:59 +02:00
[z,p,k] = zpkdata(H2)
2021-04-28 15:59:35 +02:00
#+end_src
Looking at the zeros, we see two low frequency complex conjugate zeros.
#+begin_src matlab :results output replace :exports results
2021-04-29 15:03:59 +02:00
z{1}
2021-04-28 15:59:35 +02:00
#+end_src
#+RESULTS :
#+begin_example
z{1}
ans =
-4690930.24283199 + 0i
-163.420524657426 + 0i
-0.853192261081498 + 0.713416012479897i
-0.853192261081498 - 0.713416012479897i
-3.15812268762265 + 0i
#+end_example
We manually put these zeros at the origin:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
z{1}([3,4]) = 0;
2021-04-28 15:59:35 +02:00
#+end_src
And we create a modified filter $H_{2z}(s)$:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H2z = zpk(z,p,k);
2021-04-28 15:59:35 +02:00
#+end_src
And as usual, $H_{1z}(s)$ is defined as the complementary of $H_ {2z}(s)$:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H1z = 1 - H2z;
2021-04-28 15:59:35 +02:00
#+end_src
The bode plots of $H_1(s)$, $H_2(s)$, $H_ {1z}(s)$ and $H_{2z}(s)$ are shown in Figure [[fig:comp_filters_shift_zero ]].
And we see that $H_{1z}(s)$ is slightly modified when setting the zeros at the origin for $H_ {2z}(s)$.
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
freqs = logspace(-2, 4, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, abs(squeeze(freqresp(1-H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
set(gca,'ColorOrderIndex',1)
plot(freqs, abs(squeeze(freqresp(1-H2z, freqs, 'Hz'))), '--', 'DisplayName', '$H_{1z}$');
plot(freqs, abs(squeeze(freqresp(H2z, freqs, 'Hz'))), '--', 'DisplayName', '$H_{2z}$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(1-H2, freqs, 'Hz'))), '-');
plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-');
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(H1z, freqs, 'Hz'))), '--');
plot(freqs, 360+180/pi*phase(squeeze(freqresp(H2z, freqs, 'Hz'))), '--');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-180:90:180]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/comp_filters_shift_zero.pdf', 'width', 'wide', 'height', 'tall');
2021-04-28 15:59:35 +02:00
#+end_src
#+name : fig:comp_filters_shift_zero
#+caption : Bode plots of $H_1(s)$, $H_2(s)$, $H_{1z}(s)$ and $H_{2z}(s)$
#+RESULTS :
[[file:figs/comp_filters_shift_zero.png ]]
** Imposing a positive slope at DC during the synthesis phase
Suppose we want to synthesize $H_2(s)$ such that it has a slope of +2 from DC.
We can include this "feature" in the generalized plant as shown in Figure [[fig:h_infinity_arch_H2_feature ]].
#+begin_src latex :file h_infinity_arch_H2_feature.pdf
\begin{tikzpicture}
\node[block={4.5cm}{4.0cm}, fill=black!20!white, dashed] (P) {};
\node[above] at (P.north) {$P(s)$};
\coordinate[] (inputw) at ($(P.south west)!0.8!(P.north west) + (-0.7, 0)$);
\coordinate[] (inputu) at ($(P.south west)!0.5!(P.north west) + (-0.7, 0)$);
\coordinate[] (output1) at ($(P.south east)!0.8!(P.north east) + ( 0.7, 0)$);
\coordinate[] (output2) at ($(P.south east)!0.5!(P.north east) + ( 0.7, 0)$);
\coordinate[] (outputv) at ($(P.south east)!0.2!(P.north east) + ( 0.7, 0)$);
\node[block, left=1.4 of output1] (W1){$W_1(s)$};
\node[block, left=1.4 of output2] (W2){$W_2(s)$};
\node[block, left=1.4 of outputv] (Hw){$H_{2w}(s)$};
\node[addb={+}{}{}{}{-}, left=of W1] (sub) {};
\node[block, below=0.3 of P] (H2) {$H_2^\prime(s)$};
\draw[->] (inputw) node[above right]{$w$} -- (sub.west);
\draw[->] (H2.west) -| ($(inputu)+(0.35, 0)$) node[above]{$u$} -- (W2.west);
\draw[->] (inputu-|sub) node[branch]{} -- (sub.south);
\draw[->] (sub.east) -- (W1.west);
\draw[->] ($(sub.west)+(-0.6, 0)$) node[branch]{} |- (Hw.west);
\draw[->] (Hw.east) -- ($(outputv)+(-0.35, 0)$) node[above]{$v$} |- (H2.east);
\draw[->] (W1.east) -- (output1)node[above left]{$z_1$};
\draw[->] (W2.east) -- (output2)node[above left]{$z_2$};
\end{tikzpicture}
#+end_src
#+name : fig:h_infinity_arch_H2_feature
#+caption : Generalized plant with included wanted feature represented by $H_{2w}(s)$
#+RESULTS :
[[file:figs/h_infinity_arch_H2_feature.png ]]
After synthesis, the obtained filter will be:
\begin{equation}
H_2(s) = H_2^\prime(s) H_ {2w}(s)
\end{equation}
and therefore the "feature" will be included in the filter.
For $H_1(s)$ nothing is changed: $H_1(s) = 1 - H_2(s)$.
The weighting functions are defined as usual:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
n = 3; w0 = 2*pi*10; G0 = 1e4; G1 = 0.1; Gc = 1/2;
W2 = (((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;
2021-04-28 15:59:35 +02:00
#+end_src
The wanted feature here is a +2 slope at low frequency.
For that, we use an high pass filter with a slope of +2 at low frequency.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
w0 = 2*pi*50;
H2w = (s/w0/ (s/w0+1))^2;
2021-04-28 15:59:35 +02:00
#+end_src
We define the generalized plant as shown in Figure [[fig:h_infinity_arch_H2_feature ]].
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P = [W1 -W1;
0 W2;
H2w 0];
2021-04-28 15:59:35 +02:00
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
#+begin_src matlab :results output replace :exports both
2021-04-29 15:03:59 +02:00
[H2p, ~, gamma, ~ ] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
2021-04-28 15:59:35 +02:00
#+end_src
Finally, we define $H_2(s)$ as the product of the synthesized filter and the wanted "feature":
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H2 = H2p*H2w;
2021-04-28 15:59:35 +02:00
#+end_src
And we define $H_1(s)$ to be the complementary of $H_2(s)$:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H1 = 1 - H2;
2021-04-28 15:59:35 +02:00
#+end_src
The obtained complementary filters are shown in Figure [[fig:comp_filters_H2_feature ]].
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
freqs = logspace(-3, 3, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 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$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-');
set(gca,'ColorOrderIndex',2)
plot(freqs, 360+180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-180:90:180]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/comp_filters_H2_feature.pdf', 'width', 'wide', 'height', 'tall');
2021-04-28 15:59:35 +02:00
#+end_src
#+name : fig:comp_filters_H2_feature
#+caption : Obtained complementary fitlers
#+RESULTS :
[[file:figs/comp_filters_H2_feature.png ]]
** Imposing a negative slope at infinity frequency during the synthesis phase
Let's suppose we now want to shape a low pass filter that as a negative slope until infinite frequency.
The used technique is the same as in the previous section, and the generalized plant is shown in Figure [[fig:h_infinity_arch_H2_feature ]].
The weights are defined as usual.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W1 = (((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;
2021-04-28 15:59:35 +02:00
2021-04-29 15:03:59 +02:00
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W2 = (((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;
2021-04-28 15:59:35 +02:00
#+end_src
This time, the feature is a low pass filter with a slope of -2 at high frequency.
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H2w = 1/(s/ (2*pi*10) + 1)^2;
2021-04-28 15:59:35 +02:00
#+end_src
The generalized plant is defined:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
P = [W1 -W1;
0 W2;
H2w 0];
2021-04-28 15:59:35 +02:00
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
#+begin_src matlab :results output replace :exports both
2021-04-29 15:03:59 +02:00
[H2p, ~, gamma, ~ ] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on');
2021-04-28 15:59:35 +02:00
#+end_src
The feature is added to the synthesized filter:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H2 = H2p*H2w;
2021-04-28 15:59:35 +02:00
#+end_src
And $H_1(s)$ is defined as follows:
#+begin_src matlab
2021-04-29 15:03:59 +02:00
H1 = 1 - H2;
2021-04-28 15:59:35 +02:00
#+end_src
The obtained complementary filters are shown in Figure [[fig:comp_filters_H2_feature_neg_slope ]].
#+begin_src matlab :exports none
2021-04-29 15:03:59 +02:00
freqs = logspace(-1, 4, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 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$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
% Phase
ax2 = nexttile;
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([-180:90:180]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
2021-04-28 15:59:35 +02:00
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
2021-04-29 15:03:59 +02:00
exportFig('figs/comp_filters_H2_feature_neg_slope.pdf', 'width', 'wide', 'height', 'tall');
2021-04-28 15:59:35 +02:00
#+end_src
#+name : fig:comp_filters_H2_feature_neg_slope
#+caption : Obtained complementary fitlers
#+RESULTS :
[[file:figs/comp_filters_H2_feature_neg_slope.png ]]
2020-10-26 21:35:47 +01:00
* Bibliography :ignore:
bibliographystyle:unsrt
bibliography:ref.bib
2021-04-29 15:03:59 +02:00
* Functions
** =generateWF=: Generate Weighting Functions
:PROPERTIES:
:header-args:matlab+: :tangle matlab/src/generateWF.m
:header-args:matlab+: :comments none :mkdirp yes :eval no
:END:
<<sec:generateWF >>
This Matlab function is accessible [[file:matlab/src/generateWF.m ][here ]].
*** Function description
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
function [W] = generateWF(args)
% createWeight -
%
% Syntax: [W] = generateWeight(args)
%
% Inputs:
% - n - Weight Order (integer)
% - G0 - Low frequency Gain
% - G1 - High frequency Gain
% - Gc - Gain of the weight at frequency w0
% - w0 - Frequency at which |W(j w0)| = Gc [rad/s]
%
% Outputs:
% - W - Generated Weight
#+end_src
*** Optional Parameters
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
arguments
args.n (1,1) double {mustBeInteger, mustBePositive} = 1
args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1
args.Ginf (1,1) double {mustBeNumeric, mustBePositive} = 10
args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1
args.w0 (1,1) double {mustBeNumeric, mustBePositive} = 1
end
mustBeBetween(args.G0, args.Gc, args.Ginf);
#+end_src
*** Generate the Weighting function
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
s = zpk('s');
#+end_src
#+begin_src matlab
W = (((1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/ (1-(args.Gc/args.Ginf)^(2/args.n)))*s + ...
(args.G0/args.Gc)^(1/args.n))/ ...
((1/args.Ginf)^(1/args.n)*(1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/ (1-(args.Gc/args.Ginf)^(2/args.n)))*s + ...
(1/args.Gc)^(1/args.n)))^args.n;
#+end_src
*** Verification of the $G_0$, $G_c$ and $G_\infty$ gains
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
% Custom validation function
function mustBeBetween(a,b,c)
if ~((a > b && b > c) || (c > b && b > a))
eid = 'createWeight:inputError';
msg = 'Gc should be between G0 and Ginf.';
throwAsCaller(MException(eid,msg))
end
#+end_src
** =generateCF=: Generate Complementary Filters
:PROPERTIES:
:header-args:matlab+: :tangle matlab/src/generateCF.m
:header-args:matlab+: :comments none :mkdirp yes :eval no
:END:
<<sec:generateCF >>
This Matlab function is accessible [[file:matlab/src/generateCF.m ][here ]].
*** Function description
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
function [H1, H2] = generateCF(W1, W2, args)
% createWeight -
%
% Syntax: [W] = generateCF(args)
%
% Inputs:
% - W1 - Weighting Function for H1
% - W2 - Weighting Function for H2
2021-04-30 11:18:08 +02:00
% - args:
% - method - H-Infinity solver ('lmi' or 'ric')
% - display - Display synthesis results ('on' or 'off')
2021-04-29 15:03:59 +02:00
%
% Outputs:
% - H1 - Generated H1 Filter
% - H2 - Generated H2 Filter
#+end_src
*** Optional Parameters
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
arguments
W1
W2
args.method char {mustBeMember(args.method,{'lmi', 'ric'})} = 'ric'
args.display char {mustBeMember(args.display,{'on', 'off'})} = 'on'
end
#+end_src
*** H-Infinity Synthesis
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
P = [W1 -W1;
0 W2;
1 0];
#+end_src
#+begin_src matlab :results output replace :exports both
[H2, ~, gamma, ~ ] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', args.method, 'DISPLAY', args.display);
#+end_src
#+begin_src matlab
H1 = 1 - H2;
#+end_src
