Completely rework matlab files

This commit is contained in:
Thomas Dehaeze 2021-09-01 15:00:04 +02:00
parent c073d945d1
commit 5592f7a37d
8 changed files with 2975 additions and 530 deletions

98
matlab/.latexmkrc Normal file
View File

@ -0,0 +1,98 @@
#!/bin/env perl
# Shebang is only to get syntax highlighting right across GitLab, GitHub and IDEs.
# This file is not meant to be run, but read by `latexmk`.
# ======================================================================================
# Perl `latexmk` configuration file
# ======================================================================================
# ======================================================================================
# PDF Generation/Building/Compilation
# ======================================================================================
@default_files=('dehaeze21_desig_compl_filte_matlab.tex');
# PDF-generating modes are:
# 1: pdflatex, as specified by $pdflatex variable (still largely in use)
# 2: postscript conversion, as specified by the $ps2pdf variable (useless)
# 3: dvi conversion, as specified by the $dvipdf variable (useless)
# 4: lualatex, as specified by the $lualatex variable (best)
# 5: xelatex, as specified by the $xelatex variable (second best)
$pdf_mode = 1;
# Treat undefined references and citations as well as multiply defined references as
# ERRORS instead of WARNINGS.
# This is only checked in the *last* run, since naturally, there are undefined references
# in initial runs.
# This setting is potentially annoying when debugging/editing, but highly desirable
# in the CI pipeline, where such a warning should result in a failed pipeline, since the
# final document is incomplete/corrupted.
#
# However, I could not eradicate all warnings, so that `latexmk` currently fails with
# this option enabled.
# Specifically, `microtype` fails together with `fontawesome`/`fontawesome5`, see:
# https://tex.stackexchange.com/a/547514/120853
# The fix in that answer did not help.
# Setting `verbose=silent` to mute `microtype` warnings did not work.
# Switching between `fontawesome` and `fontawesome5` did not help.
$warnings_as_errors = 0;
# Show used CPU time. Looks like: https://tex.stackexchange.com/a/312224/120853
$show_time = 1;
# Default is 5; we seem to need more owed to the complexity of the document.
# Actual documents probably don't need this many since they won't use all features,
# plus won't be compiling from cold each time.
$max_repeat=7;
# --shell-escape option (execution of code outside of latex) is required for the
#'svg' package.
# It converts raw SVG files to the PDF+PDF_TEX combo using InkScape.
#
# SyncTeX allows to jump between source (code) and output (PDF) in IDEs with support
# (many have it). A value of `1` is enabled (gzipped), `-1` is enabled but uncompressed,
# `0` is off.
# Testing in VSCode w/ LaTeX Workshop only worked for the compressed version.
# Adjust this as needed. Of course, only relevant for local use, no effect on a remote
# CI pipeline (except for slower compilation, probably).
#
# %O and %S will forward Options and the Source file, respectively, given to latexmk.
#
# `set_tex_cmds` applies to all *latex commands (latex, xelatex, lualatex, ...), so
# no need to specify these each. This allows to simply change `$pdf_mode` to get a
# different engine. Check if this works with `latexmk --commands`.
set_tex_cmds("--shell-escape -interaction=nonstopmode --synctex=1 %O %S");
# Use default pdf viewer
$pdf_previewer = 'zathura';
# option 2 is same as 1 (run biber when necessary), but also deletes the
# regeneratable bbl-file in a clenaup (`latexmk -c`). Do not use if original
# bib file is not available!
$bibtex_use = 2; # default: 1
# Change default `biber` call, help catch errors faster/clearer. See
# https://web.archive.org/web/20200526101657/https://www.semipol.de/2018/06/12/latex-best-practices.html#database-entries
$biber = "biber --validate-datamodel %O %S";
# ======================================================================================
# Auxiliary Files
# ======================================================================================
# Let latexmk know about generated files, so they can be used to detect if a
# rerun is required, or be deleted in a cleanup.
# loe: List of Examples (KOMAScript)
# lol: List of Listings (`listings` and `minted` packages)
# run.xml: biber runs
# glg: glossaries log
# glstex: generated from glossaries-extra
push @generated_exts, 'loe', 'lol', 'run.xml', 'glg', 'glstex';
# Also delete the *.glstex files from package glossaries-extra. Problem is,
# that that package generates files of the form "basename-digit.glstex" if
# multiple glossaries are present. Latexmk looks for "basename.glstex" and so
# does not find those. For that purpose, use wildcard.
# Also delete files generated by gnuplot/pgfplots contour plots
# (.dat, .script, .table).
$clean_ext = "%R-*.glstex %R_contourtmp*.*";

File diff suppressed because it is too large Load Diff

View File

@ -1,13 +1,22 @@
#+TITLE: A new method of designing complementary filters for sensor fusion using the $\mathcal{H}_\infty$ synthesis - Matlab Computation #+TITLE: A new method of designing complementary filters for sensor fusion using the $\mathcal{H}_\infty$ synthesis - Matlab Computation
:DRAWER: :DRAWER:
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas
#+HTML_LINK_HOME: ../index.html #+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html #+HTML_LINK_UP: ../index.html
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="https://research.tdehaeze.xyz/css/style.css"/> #+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> #+HTML_HEAD: <script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.js"></script>
#+LaTeX_CLASS: scrreprt #+BIND: org-latex-image-default-option "scale=1"
#+LaTeX_HEADER_EXTRA: \input{/home/thomas/Cloud/org-theme/preamble.tex} #+BIND: org-latex-image-default-width ""
#+LATEX_CLASS: scrreprt
#+LATEX_CLASS_OPTIONS: [a4paper, 10pt, DIV=12, parskip=full]
#+LATEX_HEADER_EXTRA: \input{preamble.tex}
#+LATEX_HEADER_EXTRA: \addbibresource{ref.bib}
#+PROPERTY: header-args:matlab :session *MATLAB* #+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments none #+PROPERTY: header-args:matlab+ :comments none
@ -19,15 +28,27 @@
#+PROPERTY: header-args:matlab+ :output-dir figs #+PROPERTY: header-args:matlab+ :output-dir figs
:END: :END:
#+begin_export html
<hr>
<p>This report is also available as a <a href="./dehaeze21_desig_compl_filte_matlab.pdf">pdf</a>.</p>
<hr>
#+end_export
#+begin_export latex
\clearpage
#+end_export
* Introduction :ignore: * Introduction :ignore:
This file is the Matlab file for the paper cite:dehaeze21_new_method_desig_compl_filter. The present document is a companion file for the journal paper cite:dehaeze21_new_method_desig_compl_filter.
All the Matlab cite:matlab20 scripts used in the paper are here shared and explained.
This document is divided into several sections: This document is divided into the following sections also corresponding to the paper sections:
- in section [[#sec:h_inf_synthesis_complementary_filters]], the $\mathcal{H}_\infty$ synthesis is used for generating two complementary filters - Section [[#sec:h_inf_synthesis_complementary_filters]]: the shaping of complementary filters is written as an $\mathcal{H}_\infty$ optimization problem using weighting functions. The weighting function design is discussed and the method is applied for the design of a set of simple complementary filters.
- in section [[#sec:three_comp_filters]], a method using the $\mathcal{H}_\infty$ synthesis is proposed to shape three of more complementary filters - Section [[#sec:comp_filters_ligo]]: the effectiveness of the proposed complementary filter synthesis strategy is demonstrated by designing complex complementary filters used in the first isolation stage at the LIGO
- in section [[#sec:comp_filters_ligo]], the $\mathcal{H}_\infty$ synthesis is used and compared with FIR complementary filters used for LIGO - Section [[#sec:closed_loop_complementary_filters]]: complementary filters are designed using the classical feedback loop
- in section [[#sec:closed_loop_complementary_filters]] - Section [[#sec:three_comp_filters]]: the proposed design method is generalized for the design of a set of three complementary filters
- Section [[#sec:matlab_scripts]]: complete Matlab scripts and functions developed are listed
* H-Infinity synthesis of complementary filters * H-Infinity synthesis of complementary filters
:PROPERTIES: :PROPERTIES:
@ -35,10 +56,46 @@ This document is divided into several sections:
:CUSTOM_ID: sec:h_inf_synthesis_complementary_filters :CUSTOM_ID: sec:h_inf_synthesis_complementary_filters
:END: :END:
** Introduction :ignore: ** Synthesis Architecture
#+begin_note In order to generate two complementary filters with a wanted shape, the generalized plant of Figure [[fig:h_infinity_robust_fusion_plant]] can be used.
The Matlab file corresponding to this section is accessible [[file:matlab/1_synthesis_complementary_filters.m][here]]. The included weights $W_1(s)$ and $W_2(s)$ are used to specify the upper bounds of the complementary filters being generated.
#+end_note
#+name: fig:h_infinity_robust_fusion_plant
#+caption: Generalized plant used for the $\mathcal{H}_\infty$ synthesis of a set of two complementary fiters
[[file:figs-journal/h_infinity_robust_fusion_plant.png]]
Applied the standard $\mathcal{H}_\infty$ synthesis on this generalized plant will give a transfer function $H_2(s)$ (see Figure [[fig:h_infinity_robust_fusion_fb]]) such that the $\mathcal{H}_\infty$ norm of the transfer function from $w$ to $[z_1,\ z_2]$ is less than one eqref:eq:h_inf_objective.
#+name: fig:h_infinity_robust_fusion_fb
#+caption: Generalized plant with the synthesized filter obtained after the $\mathcal{H}_\infty$ synthesis
[[file:figs-journal/h_infinity_robust_fusion_fb.png]]
\begin{equation}
\left\| \begin{array}{c} (1 - H_2(s)) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_\infty < 1 \label{eq:h_inf_objective}
\end{equation}
Thus, if the synthesis is successful and the above condition is verified, we can define $H_1(s)$ to be the complementary of $H_2(s)$ eqref:eq:H1_complementary_of_H2 and we have condition eqref:eq:shaping_comp_filters verified.
\begin{equation}
H_1(s) = 1 - H_2(s) \label{eq:H1_complementary_of_H2}
\end{equation}
\begin{equation}
\left\| \begin{array}{c} H_1(s) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_\infty < 1
\quad \Longrightarrow \quad
\left\{
\begin{array}{c}
|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{array}
\right. \label{eq:shaping_comp_filters}
\end{equation}
We then see that $W_1(s)$ and $W_2(s)$ can be used to set the wanted upper bounds of the magnitudes of both $H_1(s)$ and $H_2(s)$.
The presented synthesis method therefore allows to shape two filters $H_1(s)$ and $H_2(s)$ eqref:eq:shaping_comp_filters while ensuring their complementary property eqref:eq:H1_complementary_of_H2.
The complete Matlab script for this part is given in Section [[#sec:1_synthesis_complementary_filters]].
** Matlab Init :noexport:ignore: ** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
@ -64,37 +121,14 @@ addpath('./matlab/src');
addpath('./src'); addpath('./src');
#+end_src #+end_src
** Synthesis Architecture ** Design of Weighting Function - Proposed formula
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
[[file:figs-journal/h_infinity_robust_fusion_plant.png]]
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: A formula is proposed to help the design of the weighting functions:
\begin{equation} \begin{equation}
W(s) = \left( \frac{ 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}} \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}} \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 }\right)^n \label{eq:weighting_function_formula}
\end{equation} \end{equation}
The parameters permits to specify: The parameters permits to specify:
@ -152,11 +186,16 @@ exportFig('figs/weight_formula.pdf', 'width', 'wide', 'height', 'normal');
#+end_src #+end_src
#+name: fig:weight_formula #+name: fig:weight_formula
#+caption: Gain of the Weighting Function formula #+caption: Magnitude of the weighting function generated using formula eqref:eq:weighting_function_formula
#+RESULTS: #+RESULTS:
[[file:figs/weight_formula.png]] [[file:figs/weight_formula.png]]
** Example ** Weighting functions for the design of two complementary filters
<<sec:weighting_functions_example>>
The weighting function formula eqref:eq:weighting_function_formula is used to generate the upper bounds of two complementary filters that we wish to design.
The matlab function =generateWF= is described in Section [[#sec:generateWF]].
#+begin_src matlab #+begin_src matlab
%% Design of the Weighting Functions %% Design of the Weighting Functions
@ -164,6 +203,8 @@ W1 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.45);
W2 = generateWF('n', 2, 'w0', 2*pi*10, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.45); W2 = generateWF('n', 2, 'w0', 2*pi*10, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.45);
#+end_src #+end_src
The inverse magnitude of these two weighting functions are shown in Figure [[fig:weights_W1_W2]].
#+begin_src matlab :exports none #+begin_src matlab :exports none
%% Plot of the Weighting function magnitude %% Plot of the Weighting function magnitude
figure; figure;
@ -196,12 +237,12 @@ exportFig('figs/weights_W1_W2.pdf', 'width', 'wide', 'height', 'normal');
#+end_src #+end_src
#+name: fig:weights_W1_W2 #+name: fig:weights_W1_W2
#+caption: Weights on the complementary filters $W_1$ and $W_2$ and the associated performance weights #+caption: Inverse magnitude of the design weighting functions
#+RESULTS: #+RESULTS:
[[file:figs/weights_W1_W2.png]] [[file:figs/weights_W1_W2.png]]
** H-Infinity Synthesis ** Synthesis of the complementary filters
We define the generalized plant $P$ on matlab. The generalized plant of Figure [[fig:h_infinity_robust_fusion_plant]] is defined as follows:
#+begin_src matlab #+begin_src matlab
%% Generalized Plant %% Generalized Plant
P = [W1 -W1; P = [W1 -W1;
@ -209,7 +250,7 @@ P = [W1 -W1;
1 0]; 1 0];
#+end_src #+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command. And the $\mathcal{H}_\infty$ synthesis is performed using the =hinfsyn= command.
#+begin_src matlab :results output replace :exports both #+begin_src matlab :results output replace :exports both
%% H-Infinity Synthesis %% H-Infinity Synthesis
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); [H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
@ -242,20 +283,23 @@ And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
Best performance (actual): 0.9897 Best performance (actual): 0.9897
#+end_example #+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]]. As shown above, the obtained $\mathcal{H}_\infty$ norm of the transfer function from $w$ to $[z_1,\ z_2]$ is found to be less than one meaning the synthesis is successful.
We then define the filter $H_1(s)$ to be the complementary of $H_2(s)$ eqref:eq:H1_complementary_of_H2.
#+begin_src matlab #+begin_src matlab
%% Define H1 to be the complementary of H2 %% Define H1 to be the complementary of H2
H1 = 1 - H2; H1 = 1 - H2;
#+end_src #+end_src
Or one can just used to =generateCF= Matlab function: The function =generateCF= can also be used to synthesize the complementary filters.
This function is described in Section [[#sec:generateCF]].
#+begin_src matlab :tangle no :eval no #+begin_src matlab :tangle no :eval no
[H1, H2] = generateCF(W1, W2); [H1, H2] = generateCF(W1, W2);
#+end_src #+end_src
** Obtained Complementary Filters ** Obtained Complementary Filters
The obtained complementary filters are shown on figure [[fig:hinf_filters_results]]. The obtained complementary filters are shown below and are found to be of order 5.
Their bode plots are shown in figure [[fig:hinf_filters_results]] and compare with the defined upper bounds.
#+begin_src matlab :results output replace :exports results :tangle no #+begin_src matlab :results output replace :exports results :tangle no
zpk(H1) zpk(H1)
@ -343,13 +387,9 @@ exportFig('figs/hinf_filters_results.pdf', 'width', 'wide', 'height', 'tall');
:END: :END:
** Introduction :ignore: ** Introduction :ignore:
#+begin_note In this section, the proposed method for the design of complementary filters is validated for the design of a set of two complex complementary filters used for the first isolation stage at the LIGO cite:hua04_low_ligo.
The Matlab file corresponding to this section is accessible [[file:matlab/2_ligo_complementary_filters.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 complete Matlab script for this part is given in Section [[#sec:2_ligo_complementary_filters]].
The FIR complementary filters designed in cite:hua05_low_ligo are of order 512.
** Matlab Init :noexport:ignore: ** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
@ -376,13 +416,13 @@ addpath('./src');
#+end_src #+end_src
** Specifications ** Specifications
The specifications for the filters are: The specifications for the set of complementary filters ($L_1,H_1$) used at the LIGO are summarized below (for further details, refer to cite:hua04_polyp_fir_compl_filter_contr_system):
1. From $0$ to $0.008\text{ Hz}$,the magnitude of the filters transfer function should be less than or equal to $8 \times 10^{-3}$ - From 0 to 0.008 Hz, the magnitude $|L_1(j\omega)|$ should be less or equal to $8 \times 10^{-4}$
2. From $0.008\text{ Hz}$ to $0.04\text{ Hz}$, it attenuates the input signal proportional to frequency cubed - Between 0.008 Hz to 0.04 Hz, the filter $L_1(s)$ should attenuate 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 - Between 0.04 Hz to 0.1 Hz, the magnitude $|L_1(j\omega)|$ 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 - Above 0.1 Hz, the magnitude $|H_1(j\omega)|$ should be less than $0.045$
The specifications are translated in upper bounds of the complementary filters are shown on figure [[fig:ligo_specifications]]. The specifications are translated into upper bounds of the complementary filters and are shown in Figure [[fig:ligo_specifications]].
#+begin_src matlab :exports none #+begin_src matlab :exports none
%% Upper bounds for the complementary filters %% Upper bounds for the complementary filters
@ -415,10 +455,9 @@ exportFig('figs/ligo_specifications.pdf', 'width', 'wide', 'height', 'normal');
[[file:figs/ligo_specifications.png]] [[file:figs/ligo_specifications.png]]
** FIR Filter ** FIR Filter
We here try to implement the FIR complementary filter synthesis as explained in cite:hua05_low_ligo. To replicated the complementary filters developed in cite:hua04_low_ligo, the CVX Matlab toolbox cite:grant14_cvx is used.
For that, we use the [[http://cvxr.com/cvx/][CVX matlab Toolbox]].
We setup the CVX toolbox and use the =SeDuMi= solver. The CVX toolbox is initialized and the =SeDuMi= solver cite:sturm99_using_sedum is used.
#+begin_src matlab #+begin_src matlab
%% Initialized CVX %% Initialized CVX
cvx_startup; cvx_startup;
@ -434,7 +473,7 @@ w3 = 0.04:8.12e-4:0.1;
w4 = 0.1:8.12e-4:0.83; w4 = 0.1:8.12e-4:0.83;
#+end_src #+end_src
We then define the order of the FIR filter. The order $n$ of the FIR filter is defined.
#+begin_src matlab #+begin_src matlab
%% Filter order %% Filter order
n = 512; n = 512;
@ -453,7 +492,7 @@ 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]))]; B4 = [zeros(length(w4),1), sin(kron(w4'.*(2*pi),[1:n-1]))];
#+end_src #+end_src
We run the convex optimization. And the convex optimization is run.
#+begin_src matlab :results output replace :wrap example #+begin_src matlab :results output replace :wrap example
%% Convex optimization %% Convex optimization
cvx_begin cvx_begin
@ -561,7 +600,7 @@ Optimal value (cvx_optval): -0.0433341
h = y(2:end); h = y(2:end);
#+end_example #+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. Finally, the filter response is computed 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:hua04_low_ligo.
#+begin_src matlab #+begin_src matlab
%% Combine the frequency vectors to form the obtained filter %% Combine the frequency vectors to form the obtained filter
@ -621,17 +660,18 @@ exportFig('figs/fir_filter_ligo.pdf', 'width', 'wide', 'height', 'tall');
#+RESULTS: #+RESULTS:
[[file:figs/fir_filter_ligo.png]] [[file:figs/fir_filter_ligo.png]]
** Weights ** Weighting function design
We design weights that will be used for the $\mathcal{H}_\infty$ synthesis of the complementary filters. The weightings function that will be used for the $\mathcal{H}_\infty$ synthesis of the complementary filters are now designed.
These weights will determine the order of the obtained filters. These weights will determine the order of the obtained filters.
Here are the requirements on the filters: Here are the requirements on the filters:
- reasonable order - reasonable order
- to be as close as possible to the specified upper bounds - to be as close as possible to the specified upper bounds
- stable minimum phase - stable and minimum phase
The bode plot of the weights is shown on figure [[fig:ligo_weights]]. The weighting function for the High Pass filter is defined as follows:
#+begin_src matlab
#+begin_src matlab :exports none
%% Design of the weight for the high pass filter %% Design of the weight for the high pass filter
w1 = 2*pi*0.008; x1 = 0.35; w1 = 2*pi*0.008; x1 = 0.35;
w2 = 2*pi*0.04; x2 = 0.5; w2 = 2*pi*0.04; x2 = 0.5;
@ -651,7 +691,8 @@ wH = 1/wH;
wH = minreal(ss(wH)); wH = minreal(ss(wH));
#+end_src #+end_src
#+begin_src matlab :exports none And the weighting function for the Low pass filter is taken as a Chebyshev Type I filter.
#+begin_src matlab
%% Design of the weight for the low pass filter %% Design of the weight for the low pass filter
n = 20; % Filter order n = 20; % Filter order
Rp = 1; % Peak to peak passband ripple Rp = 1; % Peak to peak passband ripple
@ -665,6 +706,8 @@ wL = 1/wL;
wL = minreal(ss(wL)); wL = minreal(ss(wL));
#+end_src #+end_src
The inverse magnitude of the weighting functions are shown in Figure [[fig:ligo_weights]].
#+begin_src matlab :exports none #+begin_src matlab :exports none
%% Magnitude of the designed Weights and initial specifications %% Magnitude of the designed Weights and initial specifications
figure; figure;
@ -701,8 +744,8 @@ exportFig('figs/ligo_weights.pdf', 'width', 'wide', 'height', 'normal');
#+RESULTS: #+RESULTS:
[[file:figs/ligo_weights.png]] [[file:figs/ligo_weights.png]]
** H-Infinity Synthesis ** Synthesis of the complementary filters
We define the generalized plant as shown on figure [[fig:h_infinity_robst_fusion]]. The generalized plant of figure [[fig:h_infinity_robust_fusion_plant]] is defined.
#+begin_src matlab #+begin_src matlab
%% Generalized plant for the H-infinity Synthesis %% Generalized plant for the H-infinity Synthesis
P = [0 wL; P = [0 wL;
@ -710,7 +753,7 @@ P = [0 wL;
1 0]; 1 0];
#+end_src #+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command. And the standard $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command is performed.
#+begin_src matlab :results output replace :exports both :wrap example #+begin_src matlab :results output replace :exports both :wrap example
%% Standard H-Infinity synthesis %% Standard H-Infinity synthesis
[Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); [Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
@ -740,7 +783,13 @@ Test bounds: 0.3276 < gamma <= 1.8063
Gamma value achieved: 1.0439 Gamma value achieved: 1.0439
#+end_example #+end_example
The high pass filter is defined as $H_H = 1 - H_L$. The obtained $\mathcal{H}_\infty$ norm is found to be close than one meaning the synthesis is successful.
The high pass filter $H_H(s)$ is defined to be the complementary of the synthesized low pass filter $H_L(s)$:
\begin{equation}
H_H(s) = 1 - H_L(s)
\end{equation}
#+begin_src matlab #+begin_src matlab
%% High pass filter as the complementary of the low pass filter %% High pass filter as the complementary of the low pass filter
Hh = 1 - Hl; Hh = 1 - Hl;
@ -752,7 +801,7 @@ Hh = minreal(Hh);
Hl = minreal(Hl); Hl = minreal(Hl);
#+end_src #+end_src
The size of the filters is shown below. The size of the filters is shown to be equal to the sum of the weighting functions orders.
#+begin_src matlab :exports results :results output replace :wrap example :tangle no #+begin_src matlab :exports results :results output replace :wrap example :tangle no
size(Hh), size(Hl) size(Hh), size(Hl)
@ -765,7 +814,7 @@ State-space model with 1 outputs, 1 inputs, and 27 states.
State-space model with 1 outputs, 1 inputs, and 27 states. State-space model with 1 outputs, 1 inputs, and 27 states.
#+end_example #+end_example
The bode plot of the obtained filters as shown on figure [[fig:hinf_synthesis_ligo_results]]. The magnitude of the obtained filters as well as the requirements are shown in Figure [[fig:hinf_synthesis_ligo_results]].
#+begin_src matlab :exports none #+begin_src matlab :exports none
%% Bode plot of the obtained filters and comparison with the upper bounds %% Bode plot of the obtained filters and comparison with the upper bounds
@ -799,8 +848,11 @@ exportFig('figs/hinf_synthesis_ligo_results.pdf', 'width', 'wide', 'height', 'no
#+RESULTS: #+RESULTS:
[[file:figs/hinf_synthesis_ligo_results.png]] [[file:figs/hinf_synthesis_ligo_results.png]]
** Compare FIR and H-Infinity Filters ** Comparison of the FIR filters and synthesized filters
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]].
Let's now compare the FIR filters designed in cite:hua04_low_ligo with the with complementary filters obtained with the $\mathcal{H}_\infty$ synthesis.
This is done in Figure [[fig:comp_fir_ligo_hinf]], and both set of filters are found to be very close to each other.
#+begin_src matlab :exports none #+begin_src matlab :exports none
%% Comparison of the complementary filters obtained with H-infinity and with CVX %% Comparison of the complementary filters obtained with H-infinity and with CVX
@ -872,9 +924,14 @@ exportFig('figs/comp_fir_ligo_hinf.pdf', 'width', 'wide', 'height', 'tall');
:END: :END:
** Introduction :ignore: ** Introduction :ignore:
#+begin_note
The Matlab file corresponding to this section is accessible [[file:matlab/3_closed_loop_complementary_filters.m][here]]. In this section, the classical feedback architecture shown in Figure [[fig:feedback_sensor_fusion]] is used for the design of complementary filters.
#+end_note
#+name: fig:feedback_sensor_fusion
#+caption: "Closed-Loop" complementary filters
[[file:figs-journal/feedback_sensor_fusion.png]]
The complete Matlab script for this part is given in Section [[#sec:3_closed_loop_complementary_filters]].
** Matlab Init :noexport:ignore: ** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
@ -900,26 +957,62 @@ addpath('./matlab/src');
addpath('./src'); addpath('./src');
#+end_src #+end_src
** Using Feedback architecture ** Weighting Function design
Weighting functions using the =generateWF= Matlab function are designed to specify the upper bounds of the complementary filters to be designed.
These weighting functions are the same as the ones used in Section [[sec:weighting_functions_example]].
#+begin_src matlab #+begin_src matlab
%% Design of the Weighting Functions %% Design of the Weighting Functions
W1 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.45); W1 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.45);
W2 = generateWF('n', 2, 'w0', 2*pi*10, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.45); W2 = generateWF('n', 2, 'w0', 2*pi*10, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.45);
#+end_src #+end_src
Let's first synthesize $H_1(s)$: ** Generalized plant
The generalized plant of Figure [[fig:feedback_synthesis_architecture_generalized_plant]] is defined below:
#+begin_src matlab #+begin_src matlab
%% Generalized plant for "closed-loop" complementary filter synthesis %% Generalized plant for "closed-loop" complementary filter synthesis
P = [ W1 0 1; P = [ W1 0 1;
-W1 W2 -1]; -W1 W2 -1];
#+end_src #+end_src
#+name: fig:feedback_synthesis_architecture_generalized_plant
#+caption: Generalized plant used for the $\mathcal{H}_\infty$ synthesis of "closed-loop" complementary filters
[[file:figs-journal/feedback_synthesis_architecture_generalized_plant.png]]
** Synthesis of the closed-loop complementary filters
And the standard $\mathcal{H}_\infty$ synthesis is performed.
#+begin_src matlab :results output replace :exports both #+begin_src matlab :results output replace :exports both
%% Standard H-Infinity Synthesis %% Standard H-Infinity Synthesis
[L, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); [L, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
#+end_src #+end_src
#+RESULTS:
#+begin_example
Test bounds: 0.3191 <= gamma <= 1.669
gamma X>=0 Y>=0 rho(XY)<1 p/f
7.299e-01 -1.5e-19 -2.4e+01 # 1.555e-18 f
1.104e+00 0.0e+00 1.6e-07 2.037e-19 p
8.976e-01 -3.2e-16 -1.4e+02 # 5.561e-16 f
9.954e-01 0.0e+00 1.6e-07 1.041e-15 p
9.452e-01 -1.1e-15 -3.8e+02 # 4.267e-15 f
9.700e-01 -6.5e-16 -1.6e+03 # 9.876e-15 f
9.826e-01 0.0e+00 1.6e-07 8.775e-39 p
9.763e-01 -5.0e-16 -6.2e+03 # 3.519e-14 f
9.795e-01 0.0e+00 1.6e-07 6.971e-20 p
9.779e-01 -1.9e-31 -2.2e+04 # 5.600e-18 f
9.787e-01 0.0e+00 1.6e-07 5.546e-19 p
Limiting gains...
9.789e-01 0.0e+00 1.6e-07 1.084e-13 p
9.789e-01 0.0e+00 9.7e-07 1.137e-13 p
Best performance (actual): 0.9789
#+end_example
** Synthesized filters
The obtained filter $L(s)$ can then be included in the feedback architecture shown in Figure [[fig:hinf_filters_results_mixed_sensitivity]].
The closed-loop transfer functions from $\hat{x}_1$ to $\hat{x}$ and from $\hat{x}_2$ to $\hat{x}$ corresponding respectively to the sensitivity and complementary sensitivity transfer functions are defined below:
#+begin_src matlab #+begin_src matlab
%% Complementary filters %% Complementary filters
H1 = inv(1 + L); H1 = inv(1 + L);
@ -944,6 +1037,8 @@ zpk(H2) =
(s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272) (s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)
#+end_example #+end_example
The bode plots of the synthesized complementary filters are compared with the upper bounds in Figure [[fig:hinf_filters_results_mixed_sensitivity]].
#+begin_src matlab :exports none #+begin_src matlab :exports none
%% Bode plot of the obtained Complementary filters with upper-bounds %% Bode plot of the obtained Complementary filters with upper-bounds
freqs = logspace(-1, 3, 1000); freqs = logspace(-1, 3, 1000);
@ -1000,21 +1095,21 @@ exportFig('figs/hinf_filters_results_mixed_sensitivity.pdf', 'width', 'wide', 'h
#+end_src #+end_src
#+name: fig:hinf_filters_results_mixed_sensitivity #+name: fig:hinf_filters_results_mixed_sensitivity
#+caption: #+caption: Bode plot of the obtained complementary filters
#+RESULTS: #+RESULTS:
[[file:figs/hinf_filters_results_mixed_sensitivity.png]] [[file:figs/hinf_filters_results_mixed_sensitivity.png]]
* Synthesis of three complementary filters * Synthesis of three complementary filters
:PROPERTIES: :PROPERTIES:
:header-args:matlab+: :tangle matlab/4_three_complementary_filters.m :header-args:matlab+: :tangle matlab/4_three_complementary_filters.m
:header-args:matlab+: :comments org :mkdirp yes
:CUSTOM_ID: sec:three_comp_filters :CUSTOM_ID: sec:three_comp_filters
:END: :END:
** Introduction :ignore: ** Introduction :ignore:
#+begin_note
The Matlab file corresponding to this section is accessible [[file:matlab/4_three_complementary_filters.m][here]]. In this section, the proposed synthesis method of complementary filters is generalized for the synthesis of a set of three complementary filters.
#+end_note
The complete Matlab script for this part is given in Section [[#sec:4_three_complementary_filters]].
** Matlab Init :noexport:ignore: ** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
@ -1038,32 +1133,35 @@ addpath('./matlab/src');
addpath('./src'); addpath('./src');
#+end_src #+end_src
** Theory ** Synthesis Architecture
We want: The synthesis objective is to shape three filters that are complementary.
\begin{align*} This corresponds to the conditions eqref:eq:obj_three_cf where $W_1(s)$, $W_2(s)$ and $W_3(s)$ are weighting functions used to specify the maximum wanted magnitude of the three complementary filters.
& |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]]. \begin{equation}
\begin{aligned}
& |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 \\
& |H_3(j\omega)| < \frac{1}{|W_3(j\omega)|}, \quad \forall\omega \\
& H_1(s) + H_2(s) + H_3(s) = 1
\end{aligned} \label{eq:obj_three_cf}
\end{equation}
This synthesis can be done by performing the standard $\mathcal{H}_\infty$ synthesis with on the generalized plant in Figure [[fig:comp_filter_three_hinf]].
#+name: fig:comp_filter_three_hinf #+name: fig:comp_filter_three_hinf
#+caption: Generalized architecture for generating 3 complementary filters #+caption: Generalized architecture for generating 3 complementary filters
[[file:figs-journal/comp_filter_three_hinf_fb.png]] [[file:figs-journal/comp_filter_three_hinf_fb.png]]
The $\mathcal{H}_\infty$ objective is: After synthesis, filter $H_2(s)$ and $H_3(s)$ are obtained as shown in Figure [[fig:comp_filter_three_hinf]].
\begin{align*} The last filter $H_1(s)$ is defined as the complementary of the two others as in eqref:eq:H1_complementary_of_H2_H3.
& |(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. \begin{equation}
H_1(s) = 1 - H_2(s) - H_3(s) \label{eq:H1_complementary_of_H2_H3}
\end{equation}
** Weights ** Weights
First we define the weights. The three weighting functions are defined as shown below.
#+begin_src matlab #+begin_src matlab
%% Design of the Weighting Functions %% Design of the Weighting Functions
W1 = generateWF('n', 2, 'w0', 2*pi*1, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.5); W1 = generateWF('n', 2, 'w0', 2*pi*1, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.5);
@ -1071,6 +1169,8 @@ 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/
W3 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.5); W3 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.5);
#+end_src #+end_src
Their inverse magnitudes are displayed in Figure [[fig:three_weighting_functions]].
#+begin_src matlab :exports none #+begin_src matlab :exports none
%% Inverse magnitude of the weighting functions %% Inverse magnitude of the weighting functions
figure; figure;
@ -1099,7 +1199,9 @@ exportFig('figs/three_weighting_functions.pdf', 'width', 'wide', 'height', 'norm
[[file:figs/three_weighting_functions.png]] [[file:figs/three_weighting_functions.png]]
** H-Infinity Synthesis ** H-Infinity Synthesis
Then we create the generalized plant =P=.
The generalized plant in Figure [[fig:comp_filter_three_hinf]] containing the weighting functions is defined below.
#+begin_src matlab #+begin_src matlab
%% Generalized plant for the synthesis of 3 complementary filters %% Generalized plant for the synthesis of 3 complementary filters
P = [W1 -W1 -W1; P = [W1 -W1 -W1;
@ -1108,7 +1210,7 @@ P = [W1 -W1 -W1;
1 0 0]; 1 0 0];
#+end_src #+end_src
And we do the $\mathcal{H}_\infty$ synthesis. And the standard $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command is performed.
#+begin_src matlab :results output replace :exports both #+begin_src matlab :results output replace :exports both
%% Standard H-Infinity Synthesis %% Standard H-Infinity Synthesis
[H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); [H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
@ -1116,7 +1218,6 @@ And we do the $\mathcal{H}_\infty$ synthesis.
#+RESULTS: #+RESULTS:
#+begin_example #+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 Resetting value of Gamma min based on D_11, D_12, D_21 terms
Test bounds: 0.1000 < gamma <= 1050.0000 Test bounds: 0.1000 < gamma <= 1050.0000
@ -1148,15 +1249,21 @@ Test bounds: 0.1000 < gamma <= 1050.0000
Gamma value achieved: 0.9936 Gamma value achieved: 0.9936
#+end_example #+end_example
** Obtained Complementary Filters The two synthesized filters $H_2(s)$ and $H_3(s)$ are defined below:
The obtained filters are: #+begin_src matlab :exports none
#+begin_src matlab %% Synthesized H2 and H3 filters
%%
H2 = tf(H(1)); H2 = tf(H(1));
H3 = tf(H(2)); H3 = tf(H(2));
#+end_src
And the third filter $H_1(s)$ is defined using eqref:eq:H1_complementary_of_H2_H3.
#+begin_src matlab
%% H1 is defined as the complementary filter of H2 and H3
H1 = 1 - H2 - H3; H1 = 1 - H2 - H3;
#+end_src #+end_src
The bode plots of the three obtained complementary filters are shown in Figure [[fig:three_complementary_filters_results]].
#+begin_src matlab :exports none #+begin_src matlab :exports none
%% Bode plot of the obtained complementary filters %% Bode plot of the obtained complementary filters
figure; figure;
@ -1217,31 +1324,53 @@ exportFig('figs/three_complementary_filters_results.pdf', 'width', 'wide', 'heig
#+RESULTS: #+RESULTS:
[[file:figs/three_complementary_filters_results.png]] [[file:figs/three_complementary_filters_results.png]]
* Bibliography :ignore: * Matlab Scripts
:PROPERTIES: :PROPERTIES:
:CUSTOM_ID: sec:bibliography :CUSTOM_ID: sec:matlab_scripts
:END: :END:
bibliographystyle:unsrt ** =1_synthesis_complementary_filters.m=
bibliography:ref.bib :PROPERTIES:
:CUSTOM_ID: sec:1_synthesis_complementary_filters
:END:
This scripts corresponds to section 3 of cite:dehaeze21_new_method_desig_compl_filter.
* Functions #+INCLUDE: "matlab/1_synthesis_complementary_filters.m" src matlab -n
** =2_ligo_complementary_filters.m=
:PROPERTIES: :PROPERTIES:
:header-args:matlab+: :comments none :mkdirp yes :eval no :CUSTOM_ID: sec:2_ligo_complementary_filters
:CUSTOM_ID: sec:functions
:END: :END:
This scripts corresponds to section 4 of cite:dehaeze21_new_method_desig_compl_filter.
#+INCLUDE: "matlab/2_ligo_complementary_filters.m" src matlab -n
** =3_closed_loop_complementary_filters.m=
:PROPERTIES:
:CUSTOM_ID: sec:3_closed_loop_complementary_filters
:END:
This scripts corresponds to section 5.1 of cite:dehaeze21_new_method_desig_compl_filter.
#+INCLUDE: "matlab/3_closed_loop_complementary_filters.m" src matlab -n
** =4_three_complementary_filters.m=
:PROPERTIES:
:CUSTOM_ID: sec:4_three_complementary_filters
:END:
This scripts corresponds to section 5.2 of cite:dehaeze21_new_method_desig_compl_filter.
#+INCLUDE: "matlab/4_three_complementary_filters.m" src matlab -n
** =generateWF=: Generate Weighting Functions ** =generateWF=: Generate Weighting Functions
:PROPERTIES: :PROPERTIES:
:header-args:matlab+: :tangle matlab/src/generateWF.m :header-args:matlab+: :tangle matlab/src/generateWF.m
:header-args:matlab+: :comments none :mkdirp yes :eval no
:CUSTOM_ID: sec:generateWF
:END: :END:
<<sec:generateWF>>
This Matlab function is accessible [[file:matlab/src/generateWF.m][here]]. This function is used to easily generate weighting functions from classical requirements.
*** Function description
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab #+begin_src matlab
function [W] = generateWF(args) function [W] = generateWF(args)
@ -1258,13 +1387,7 @@ function [W] = generateWF(args)
% %
% Outputs: % Outputs:
% - W - Generated Weighting Function % - W - Generated Weighting Function
#+end_src
*** Optional Parameters
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
%% Argument validation %% Argument validation
arguments arguments
args.n (1,1) double {mustBeInteger, mustBePositive} = 1 args.n (1,1) double {mustBeInteger, mustBePositive} = 1
@ -1273,54 +1396,19 @@ arguments
args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1 args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1
args.w0 (1,1) double {mustBeNumeric, mustBePositive} = 1 args.w0 (1,1) double {mustBeNumeric, mustBePositive} = 1
end end
#+end_src
Verification that the parameters $G_0$, $G_c$ and $G_\infty$ are satisfy condition eqref:eq:cond_formula_1 or eqref:eq:cond_formula_2.
\begin{equation}
G_0 < 1 < G_\infty \text{ and } G_0 < G_c < G_\infty \label{eq:cond_formula_1}
\end{equation}
\begin{equation}
G_\infty < 1 < G_0 \text{ and } G_\infty < G_c < G_0 \label{eq:cond_formula_2}
\end{equation}
#+begin_src matlab
% Verification of correct relation between G0, Gc and Ginf % Verification of correct relation between G0, Gc and Ginf
mustBeBetween(args.G0, args.Gc, args.Ginf); mustBeBetween(args.G0, args.Gc, args.Ginf);
#+end_src
*** Generate the Weighting function
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
%% Initialize the Laplace variable %% Initialize the Laplace variable
s = zpk('s'); s = zpk('s');
#+end_src
The weighting function formula use is:
#+name: eq:weight_formula
\begin{equation}
W(s) = \left( \frac{
\frac{1}{\omega_c} \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_c} \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}
#+begin_src matlab
%% Create the weighting function according to formula %% Create the weighting function according to formula
W = (((1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.Ginf)^(2/args.n)))*s + ... 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))/... (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.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; (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 %% Custom validation function
function mustBeBetween(a,b,c) function mustBeBetween(a,b,c)
if ~((a > b && b > c) || (c > b && b > a)) if ~((a > b && b > c) || (c > b && b > a))
@ -1333,15 +1421,11 @@ function mustBeBetween(a,b,c)
** =generateCF=: Generate Complementary Filters ** =generateCF=: Generate Complementary Filters
:PROPERTIES: :PROPERTIES:
:header-args:matlab+: :tangle matlab/src/generateCF.m :header-args:matlab+: :tangle matlab/src/generateCF.m
:header-args:matlab+: :comments none :mkdirp yes :eval no
:CUSTOM_ID: sec:generateCF
:END: :END:
<<sec:generateCF>>
This Matlab function is accessible [[file:matlab/src/generateCF.m][here]]. This function is used to easily synthesize a set of two complementary filters using the $\mathcal{H}_\infty$ synthesis.
*** Function description
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab #+begin_src matlab
function [H1, H2] = generateCF(W1, W2, args) function [H1, H2] = generateCF(W1, W2, args)
@ -1359,13 +1443,7 @@ function [H1, H2] = generateCF(W1, W2, args)
% Outputs: % Outputs:
% - H1 - Generated H1 Filter % - H1 - Generated H1 Filter
% - H2 - Generated H2 Filter % - H2 - Generated H2 Filter
#+end_src
*** Optional Parameters
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
%% Argument validation %% Argument validation
arguments arguments
W1 W1
@ -1373,25 +1451,17 @@ arguments
args.method char {mustBeMember(args.method,{'lmi', 'ric'})} = 'ric' args.method char {mustBeMember(args.method,{'lmi', 'ric'})} = 'ric'
args.display char {mustBeMember(args.display,{'on', 'off'})} = 'on' args.display char {mustBeMember(args.display,{'on', 'off'})} = 'on'
end end
#+end_src
*** H-Infinity Synthesis
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
%% The generalized plant is defined %% The generalized plant is defined
P = [W1 -W1; P = [W1 -W1;
0 W2; 0 W2;
1 0]; 1 0];
#+end_src
#+begin_src matlab :results output replace :exports both
%% The standard H-infinity synthesis is performed %% The standard H-infinity synthesis is performed
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', args.method, 'DISPLAY', args.display); [H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', args.method, 'DISPLAY', args.display);
#+end_src
#+begin_src matlab
%% H1 is defined as the complementary of H2 %% H1 is defined as the complementary of H2
H1 = 1 - H2; H1 = 1 - H2;
#+end_src #+end_src
* Bibliography :ignore:
#+latex: \printbibliography

Binary file not shown.

File diff suppressed because it is too large Load Diff

Binary file not shown.

114
matlab/preamble.tex Normal file
View File

@ -0,0 +1,114 @@
\usepackage{float}
\usepackage{caption,tabularx,booktabs}
\usepackage{biblatex}
\usepackage{fontawesome}
\usepackage{caption}
\usepackage{subcaption}
\captionsetup[figure]{labelfont=bf}
\captionsetup[subfigure]{labelfont=bf}
\captionsetup[listing]{labelfont=bf}
\captionsetup[table]{labelfont=bf}
\usepackage{xcolor}
\definecolor{my-blue}{HTML}{6b7adb}
\definecolor{my-pale-blue}{HTML}{e6e9f9}
\definecolor{my-red}{HTML}{db6b6b}
\definecolor{my-pale-red}{HTML}{f9e6e6}
\definecolor{my-green}{HTML}{6bdbb6}
\definecolor{my-pale-green}{HTML}{e6f9f3}
\definecolor{my-yellow}{HTML}{dbd26b}
\definecolor{my-pale-yellow}{HTML}{f9f7e6}
\definecolor{my-orange}{HTML}{dba76b}
\definecolor{my-pale-orange}{HTML}{f9f0e6}
\definecolor{my-grey}{HTML}{a3a3a3}
\definecolor{my-pale-grey}{HTML}{f0f0f0}
\definecolor{my-turq}{HTML}{6bc7db}
\definecolor{my-pale-turq}{HTML}{e6f6f9}
\usepackage{inconsolata}
\usepackage[newfloat=true, chapter]{minted}
\usemintedstyle{autumn}
\setminted{frame=lines,breaklines=true,tabsize=4,fontsize=\scriptsize,autogobble=true,labelposition=topline,bgcolor=my-pale-grey}
\setminted[matlab]{label=Matlab}
\setminted[latex]{label=LaTeX}
\setminted[bash]{label=Bash}
\setminted[python]{label=Python}
\setminted[text]{label=Results}
\setminted[md]{label=Org Mode}
\setmintedinline{fontsize=\normalsize,bgcolor=my-pale-grey}
\usepackage[most]{tcolorbox}
\tcbuselibrary{minted}
\newtcolorbox{seealso}{ enhanced,breakable,colback=my-pale-grey,colframe=my-grey,fonttitle=\bfseries,title=See Also}
\newtcolorbox{hint}{ enhanced,breakable,colback=my-pale-grey,colframe=my-grey,fonttitle=\bfseries,title=Hint}
\newtcolorbox{definition}{enhanced,breakable,colback=my-pale-red, colframe=my-red, fonttitle=\bfseries,title=Definition}
\newtcolorbox{important}{ enhanced,breakable,colback=my-pale-red, colframe=my-red, fonttitle=\bfseries,title=Important}
\newtcolorbox{exampl}[1][]{ enhanced,breakable,colback=my-pale-green,colframe=my-green,fonttitle=\bfseries,title=Example,#1}
\newtcolorbox{exercice}{ enhanced,breakable,colback=my-pale-yellow,colframe=my-yellow,fonttitle=\bfseries,title=Exercice}
\newtcolorbox{question}{ enhanced,breakable,colback=my-pale-yellow,colframe=my-yellow,fonttitle=\bfseries,title=Question}
\newtcolorbox{answer}{ enhanced,breakable,colback=my-pale-turq,colframe=my-turq,fonttitle=\bfseries,title=Answer}
\newtcolorbox{summary}{ enhanced,breakable,colback=my-pale-blue,colframe=my-blue,fonttitle=\bfseries,title=Summary}
\newtcolorbox{note}{ enhanced,breakable,colback=my-pale-blue,colframe=my-blue,fonttitle=\bfseries,title=Note}
\newtcolorbox{caution}{ enhanced,breakable,colback=my-pale-orange,colframe=my-orange,fonttitle=\bfseries,title=Caution}
\newtcolorbox{warning}{ enhanced,breakable,colback=my-pale-orange,colframe=my-orange,fonttitle=\bfseries,title=Warning}
\newtcolorbox{my-quote}[1]{%
colback=my-pale-grey,
grow to right by=-10mm,
grow to left by=-10mm,
boxrule=0pt,
boxsep=0pt,
breakable,
enhanced jigsaw,
borderline west={4pt}{0pt}{my-grey}}
\renewenvironment{quote}{\begin{my-quote}}{\end{my-quote}}
\newtcolorbox{my-verse}[1]{%
colback=my-pale-grey,
grow to right by=-10mm,
grow to left by=-10mm,
boxrule=0pt,
boxsep=0pt,
breakable,
enhanced jigsaw,
borderline west={4pt}{0pt}{my-grey}}
\renewenvironment{verse}{\begin{my-verse}}{\end{my-verse}}
\usepackage{environ}% http://ctan.org/pkg/environ
\NewEnviron{aside}{%
\marginpar{\BODY}
}
\renewenvironment{verbatim}{\VerbatimEnvironment\begin{minted}[]{text}}{\end{minted}}
\usepackage{soul}
\sethlcolor{my-pale-grey}
\let\OldTexttt\texttt
\renewcommand{\texttt}[1]{{\ttfamily\hl{\mbox{\,#1\,}}}}
\makeatletter
\preto\Gin@extensions{png,}
\DeclareGraphicsRule{.png}{pdf}{.pdf}{\noexpand\Gin@base.pdf}
\makeatother
\usepackage{hyperref}
\hypersetup{
colorlinks = true,
allcolors = my-blue
}
\usepackage{hypcap}

View File

@ -1,16 +1,64 @@
@phdthesis{hua05_low_ligo, @inproceedings{hua04_low_ligo,
author = {Hua, Wensheng}, author = {Hua, Wensheng and Adhikari, R and DeBra, Daniel B and
school = {stanford university}, Giaime, Joseph A and Hammond, Giles Dominic and Hardham, C and
title = {Low frequency vibration isolation and alignment system for Hennessy, Mike and How, Jonathan P and Lantz, Brian T and
advanced LIGO}, Macinnis, M and others},
year = 2005, title = {Low-frequency active vibration isolation for advanced
{LIGO}},
booktitle = {Gravitational Wave and Particle Astrophysics Detectors},
year = 2004,
volume = 5500,
pages = {194--205},
organization = {International Society for Optics and Photonics},
} }
@article{dehaeze21_new_method_desig_compl_filter, @inproceedings{dehaeze21_new_method_desig_compl_filter,
author = {Dehaeze, Thomas and Vermat, Mohit and Collette, Christophe}, author = {Dehaeze, Thomas and Vermat, Mohit and Collette, Christophe},
title = {A New Method of Designing Complementary Filters for Sensor title = {A New Method of Designing Complementary Fil},
Fusion Using the $\mathcal{H}_\infty$ Synthesis}, booktitle = {Gravitational Wave and Particle Astrophysics Detectors},
journal = {Mechanical Systems and Signal Processing}, year = 2004,
year = 2021, volume = 5500,
month = {Nov}, pages = {194--205},
organization = {International Society for Optics and Photonics},
}
@book{matlab20,
author = {MATLAB},
title = {version 9.9.0 (R2020b)},
year = 2020,
publisher = {The MathWorks Inc.},
address = {Natick, Massachusetts},
}
@article{sturm99_using_sedum,
author = {Sturm, Jos F},
title = {Using {SeDuMi} 1.02, a Matlab Toolbox for Optimization Over
Symmetric Cones},
journal = {Optimization methods and software},
volume = 11,
number = {1-4},
pages = {625--653},
year = 1999,
doi = {10.1080/10556789908805766},
url = {https://doi.org/10.1080/10556789908805766},
publisher = {Taylor \& Francis},
}
@inproceedings{hua04_polyp_fir_compl_filter_contr_system,
author = {Hua, Wensheng and Debra, B. and Hardham, T. and Lantz, T.
and Giaime, A.},
title = {Polyphase {FIR} Complementary Filters for Control Systems},
booktitle = {Proceedings of ASPE Spring Topical Meeting on Control of
Precision Systems},
year = 2004,
pages = {109--114},
}
@misc{grant14_cvx,
author = {Michael Grant and Stephen Boyd},
howpublished = {\url{http://cvxr.com/cvx}},
month = 3,
title = {{CVX}: Matlab Software for Disciplined Convex Programming,
version 2.1},
year = 2014,
} }