tdehaeze/dehaeze21_desig_compl_filte
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
master
dehaeze21_desig_compl_filte/matlab/dehaeze21_desig_compl_filte_matlab.tex
Thomas Dehaeze 572e1d5abd Update title
2021-09-02 10:04:16 +02:00

1405 lines
52 KiB
TeX
Raw Permalink Blame History

% Created 2021-09-02 jeu. 10:03
% Intended LaTeX compiler: pdflatex
\documentclass[a4paper, 10pt, DIV=12, parskip=full]{scrreprt}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{grffile}
\usepackage{longtable}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage[normalem]{ulem}
\usepackage{amsmath}
\usepackage{textcomp}
\usepackage{amssymb}
\usepackage{capt-of}
\usepackage{hyperref}
\usepackage[most]{tcolorbox}
\usepackage{bm}
\usepackage{booktabs}
\usepackage{tabularx}
\usepackage{array}
\usepackage{siunitx}
\input{preamble.tex}
\addbibresource{ref.bib}
\author{Dehaeze Thomas}
\date{\today}
\title{Designing Complementary filters for sensor fusion using \(\mathcal{H}_\infty\) synthesis - Matlab Computation}
\begin{document}
\maketitle
\tableofcontents
\clearpage
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 the following sections also corresponding to the paper sections:
\begin{itemize}
\item Section \ref{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.
\item Section \ref{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
\item Section \ref{sec:closed_loop_complementary_filters}: complementary filters are designed using the classical feedback loop
\item Section \ref{sec:three_comp_filters}: the proposed design method is generalized for the design of a set of three complementary filters
\item Section \ref{sec:matlab_scripts}: complete Matlab scripts and functions developed are listed
\end{itemize}
\chapter{H-Infinity synthesis of complementary filters}
\label{sec:h_inf_synthesis_complementary_filters}
\section{Synthesis Architecture}
In order to generate two complementary filters with a wanted shape, the generalized plant of Figure \ref{fig:h_infinity_robust_fusion_plant} can be used.
The included weights \(W_1(s)\) and \(W_2(s)\) are used to specify the upper bounds of the complementary filters being generated.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs-journal/h_infinity_robust_fusion_plant.png}
\caption{\label{fig:h_infinity_robust_fusion_plant}Generalized plant used for the \(\mathcal{H}_\infty\) synthesis of a set of two complementary fiters}
\end{figure}
Applied the standard \(\mathcal{H}_\infty\) synthesis on this generalized plant will give a transfer function \(H_2(s)\) (see Figure \ref{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}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs-journal/h_infinity_robust_fusion_fb.png}
\caption{\label{fig:h_infinity_robust_fusion_fb}Generalized plant with the synthesized filter obtained after the \(\mathcal{H}_\infty\) synthesis}
\end{figure}
\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 \ref{sec:1_synthesis_complementary_filters}.
\section{Design of Weighting Function - Proposed formula}
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 \label{eq:weighting_function_formula}
\end{equation}
The parameters permits to specify:
\begin{itemize}
\item the low frequency gain: \(G_0 = lim_{\omega \to 0} |W(j\omega)|\)
\item the high frequency gain: \(G_\infty = lim_{\omega \to \infty} |W(j\omega)|\)
\item the absolute gain at \(\omega_0\): \(G_c = |W(j\omega_0)|\)
\item the absolute slope between high and low frequency: \(n\)
\end{itemize}
The general shape of a weighting function generated using the formula is shown in figure \ref{fig:weight_formula}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/weight_formula.png}
\caption{\label{fig:weight_formula}Magnitude of the weighting function generated using formula \eqref{eq:weighting_function_formula}}
\end{figure}
\section{Weighting functions for the design of two complementary filters}
\label{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 \texttt{generateWF} is described in Section \ref{sec:generateWF}.
\begin{minted}[]{matlab}
%% Design of the Weighting Functions
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);
\end{minted}
The inverse magnitude of these two weighting functions are shown in Figure \ref{fig:weights_W1_W2}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/weights_W1_W2.png}
\caption{\label{fig:weights_W1_W2}Inverse magnitude of the design weighting functions}
\end{figure}
\section{Synthesis of the complementary filters}
The generalized plant of Figure \ref{fig:h_infinity_robust_fusion_plant} is defined as follows:
\begin{minted}[]{matlab}
%% Generalized Plant
P = [W1 -W1;
0 W2;
1 0];
\end{minted}
And the \(\mathcal{H}_\infty\) synthesis is performed using the \texttt{hinfsyn} command.
\begin{minted}[]{matlab}
%% H-Infinity Synthesis
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
\end{minted}
\begin{verbatim}
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
Test bounds: 0.3223 <= gamma <= 1000
gamma X>=0 Y>=0 rho(XY)<1 p/f
1.795e+01 1.4e-07 0.0e+00 1.481e-16 p
2.406e+00 1.4e-07 0.0e+00 3.604e-15 p
8.806e-01 -3.1e+02 # -1.4e-16 7.370e-19 f
1.456e+00 1.4e-07 0.0e+00 1.499e-18 p
1.132e+00 1.4e-07 0.0e+00 8.587e-15 p
9.985e-01 1.4e-07 0.0e+00 2.331e-13 p
9.377e-01 -7.7e+02 # -6.6e-17 3.744e-14 f
9.676e-01 -2.0e+03 # -5.7e-17 1.046e-13 f
9.829e-01 -6.6e+03 # -1.1e-16 2.949e-13 f
9.907e-01 1.4e-07 0.0e+00 2.374e-19 p
9.868e-01 -1.6e+04 # -6.4e-17 5.331e-14 f
9.887e-01 -5.1e+04 # -1.5e-17 2.703e-19 f
9.897e-01 1.4e-07 0.0e+00 1.583e-11 p
Limiting gains...
9.897e-01 1.5e-07 0.0e+00 1.183e-12 p
9.897e-01 6.9e-07 0.0e+00 1.365e-12 p
Best performance (actual): 0.9897
\end{verbatim}
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{minted}[]{matlab}
%% Define H1 to be the complementary of H2
H1 = 1 - H2;
\end{minted}
The function \texttt{generateCF} can also be used to synthesize the complementary filters.
This function is described in Section \ref{sec:generateCF}.
\begin{minted}[]{matlab}
[H1, H2] = generateCF(W1, W2);
\end{minted}
\section{Obtained Complementary Filters}
The obtained complementary filters are shown below and are found to be of order 5.
Their bode plots are shown in figure \ref{fig:hinf_filters_results} and compare with the defined upper bounds.
\begin{verbatim}
zpk(H1)
ans =
(s+1.289e05) (s+153.6) (s+3.842)^3
-------------------------------------------------------
(s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)
zpk(H2)
ans =
125.61 (s+3358)^2 (s^2 + 46.61s + 813.8)
-------------------------------------------------------
(s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)
\end{verbatim}
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/hinf_filters_results.png}
\caption{\label{fig:hinf_filters_results}Obtained complementary filters using \(\mathcal{H}_\infty\) synthesis}
\end{figure}
\chapter{Design of complementary filters used in the Active Vibration Isolation System at the LIGO}
\label{sec:comp_filters_ligo}
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 complete Matlab script for this part is given in Section \ref{sec:2_ligo_complementary_filters}.
\section{Specifications}
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}):
\begin{itemize}
\item From 0 to 0.008 Hz, the magnitude \(|L_1(j\omega)|\) should be less or equal to \(8 \times 10^{-4}\)
\item Between 0.008 Hz to 0.04 Hz, the filter \(L_1(s)\) should attenuate the input signal proportional to frequency cubed
\item Between 0.04 Hz to 0.1 Hz, the magnitude \(|L_1(j\omega)|\) should be less than \(3\)
\item Above 0.1 Hz, the magnitude \(|H_1(j\omega)|\) should be less than \(0.045\)
\end{itemize}
The specifications are translated into upper bounds of the complementary filters and are shown in Figure \ref{fig:ligo_specifications}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/ligo_specifications.png}
\caption{\label{fig:ligo_specifications}Specification for the LIGO complementary filters}
\end{figure}
\section{FIR Filter}
To replicated the complementary filters developed in \cite{hua04_low_ligo}, the CVX Matlab toolbox \cite{grant14_cvx} is used.
The CVX toolbox is initialized and the \texttt{SeDuMi} solver \cite{sturm99_using_sedum} is used.
\begin{minted}[]{matlab}
%% Initialized CVX
cvx_startup;
cvx_solver sedumi;
\end{minted}
We define the frequency vectors on which we will constrain the norm of the FIR filter.
\begin{minted}[]{matlab}
%% Frequency vectors
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;
\end{minted}
The order \(n\) of the FIR filter is defined.
\begin{minted}[]{matlab}
%% Filter order
n = 512;
\end{minted}
\begin{minted}[]{matlab}
%% Initialization of filter responses
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]))];
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]))];
\end{minted}
And the convex optimization is run.
\begin{minted}[]{matlab}
%% Convex optimization
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
h = y(2:end);
\end{minted}
\begin{verbatim}
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{verbatim}
Finally, the filter response is computed over the frequency vector defined and the result is shown on figure \ref{fig:fir_filter_ligo} which is very close to the filters obtain in \cite{hua04_low_ligo}.
\begin{minted}[]{matlab}
%% Combine the frequency vectors to form the obtained filter
w = [w1 w2 w3 w4];
H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h;
\end{minted}
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/fir_filter_ligo.png}
\caption{\label{fig:fir_filter_ligo}FIR Complementary filters obtain after convex optimization}
\end{figure}
\section{Weighting function design}
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.
Here are the requirements on the filters:
\begin{itemize}
\item reasonable order
\item to be as close as possible to the specified upper bounds
\item stable and minimum phase
\end{itemize}
The weighting function for the High Pass filter is defined as follows:
\begin{minted}[]{matlab}
%% Design of the weight for the high pass filter
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;
wH = minreal(ss(wH));
\end{minted}
And the weighting function for the Low pass filter is taken as a Chebyshev Type I filter.
\begin{minted}[]{matlab}
%% Design of the weight for the low pass filter
n = 20; % Filter order
Rp = 1; % Peak to peak passband ripple
Wp = 2*pi*0.102; % Edge frequency
% Chebyshev Type I filter design
[b,a] = cheby1(n, Rp, Wp, 'high', 's');
wL = 0.04*tf(a, b);
wL = 1/wL;
wL = minreal(ss(wL));
\end{minted}
The inverse magnitude of the weighting functions are shown in Figure \ref{fig:ligo_weights}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/ligo_weights.png}
\caption{\label{fig:ligo_weights}Weights for the \(\mathcal{H}_\infty\) synthesis}
\end{figure}
\section{Synthesis of the complementary filters}
The generalized plant of figure \ref{fig:h_infinity_robust_fusion_plant} is defined.
\begin{minted}[]{matlab}
%% Generalized plant for the H-infinity Synthesis
P = [0 wL;
wH -wH;
1 0];
\end{minted}
And the standard \(\mathcal{H}_\infty\) synthesis using the \texttt{hinfsyn} command is performed.
\begin{minted}[]{matlab}
%% Standard H-Infinity synthesis
[Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
\end{minted}
\begin{verbatim}
[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{verbatim}
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{minted}[]{matlab}
%% High pass filter as the complementary of the low pass filter
Hh = 1 - Hl;
\end{minted}
The size of the filters is shown to be equal to the sum of the weighting functions orders.
\begin{verbatim}
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{verbatim}
The magnitude of the obtained filters as well as the requirements are shown in Figure \ref{fig:hinf_synthesis_ligo_results}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/hinf_synthesis_ligo_results.png}
\caption{\label{fig:hinf_synthesis_ligo_results}Obtained complementary filters using the \(\mathcal{H}_\infty\) synthesis}
\end{figure}
\section{Comparison of the FIR filters and synthesized filters}
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 \ref{fig:comp_fir_ligo_hinf}, and both set of filters are found to be very close to each other.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/comp_fir_ligo_hinf.png}
\caption{\label{fig:comp_fir_ligo_hinf}Comparison between the FIR filters developped for LIGO and the \(\mathcal{H}_\infty\) complementary filters}
\end{figure}
\chapter{``Closed-Loop'' complementary filters}
\label{sec:closed_loop_complementary_filters}
In this section, the classical feedback architecture shown in Figure \ref{fig:feedback_sensor_fusion} is used for the design of complementary filters.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs-journal/feedback_sensor_fusion.png}
\caption{\label{fig:feedback_sensor_fusion}``Closed-Loop'' complementary filters}
\end{figure}
The complete Matlab script for this part is given in Section \ref{sec:3_closed_loop_complementary_filters}.
\section{Weighting Function design}
Weighting functions using the \texttt{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 \ref{sec:weighting_functions_example}.
\begin{minted}[]{matlab}
%% Design of the Weighting Functions
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);
\end{minted}
\section{Generalized plant}
The generalized plant of Figure \ref{fig:feedback_synthesis_architecture_generalized_plant} is defined below:
\begin{minted}[]{matlab}
%% Generalized plant for "closed-loop" complementary filter synthesis
P = [ W1 0 1;
-W1 W2 -1];
\end{minted}
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs-journal/feedback_synthesis_architecture_generalized_plant.png}
\caption{\label{fig:feedback_synthesis_architecture_generalized_plant}Generalized plant used for the \(\mathcal{H}_\infty\) synthesis of ``closed-loop'' complementary filters}
\end{figure}
\section{Synthesis of the closed-loop complementary filters}
And the standard \(\mathcal{H}_\infty\) synthesis is performed.
\begin{minted}[]{matlab}
%% Standard H-Infinity Synthesis
[L, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
\end{minted}
\begin{verbatim}
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{verbatim}
\section{Synthesized filters}
The obtained filter \(L(s)\) can then be included in the feedback architecture shown in Figure \ref{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{minted}[]{matlab}
%% Complementary filters
H1 = inv(1 + L);
H2 = 1 - H1;
\end{minted}
\begin{verbatim}
zpk(H1) =
(s+3.842)^3 (s+153.6) (s+1.289e05)
-------------------------------------------------------
(s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)
zpk(H2) =
125.61 (s+3358)^2 (s^2 + 46.61s + 813.8)
-------------------------------------------------------
(s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)
\end{verbatim}
The bode plots of the synthesized complementary filters are compared with the upper bounds in Figure \ref{fig:hinf_filters_results_mixed_sensitivity}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/hinf_filters_results_mixed_sensitivity.png}
\caption{\label{fig:hinf_filters_results_mixed_sensitivity}Bode plot of the obtained complementary filters}
\end{figure}
\chapter{Synthesis of three complementary filters}
\label{sec:three_comp_filters}
In this section, the proposed synthesis method of complementary filters is generalized for the synthesis of a set of three complementary filters.
The complete Matlab script for this part is given in Section \ref{sec:4_three_complementary_filters}.
\section{Synthesis Architecture}
The synthesis objective is to shape three filters that are complementary.
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.
\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 \ref{fig:comp_filter_three_hinf}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs-journal/comp_filter_three_hinf_fb.png}
\caption{\label{fig:comp_filter_three_hinf}Generalized architecture for generating 3 complementary filters}
\end{figure}
After synthesis, filter \(H_2(s)\) and \(H_3(s)\) are obtained as shown in Figure \ref{fig:comp_filter_three_hinf}.
The last filter \(H_1(s)\) is defined as the complementary of the two others as in \eqref{eq:H1_complementary_of_H2_H3}.
\begin{equation}
H_1(s) = 1 - H_2(s) - H_3(s) \label{eq:H1_complementary_of_H2_H3}
\end{equation}
\section{Weights}
The three weighting functions are defined as shown below.
\begin{minted}[]{matlab}
%% Design of the Weighting Functions
W1 = generateWF('n', 2, 'w0', 2*pi*1, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.5);
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;
W3 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.5);
\end{minted}
Their inverse magnitudes are displayed in Figure \ref{fig:three_weighting_functions}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/three_weighting_functions.png}
\caption{\label{fig:three_weighting_functions}Three weighting functions used for the \(\mathcal{H}_\infty\) synthesis of the complementary filters}
\end{figure}
\section{H-Infinity Synthesis}
The generalized plant in Figure \ref{fig:comp_filter_three_hinf} containing the weighting functions is defined below.
\begin{minted}[]{matlab}
%% Generalized plant for the synthesis of 3 complementary filters
P = [W1 -W1 -W1;
0 W2 0 ;
0 0 W3;
1 0 0];
\end{minted}
And the standard \(\mathcal{H}_\infty\) synthesis using the \texttt{hinfsyn} command is performed.
\begin{minted}[]{matlab}
%% Standard H-Infinity Synthesis
[H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
\end{minted}
\begin{verbatim}
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{verbatim}
The two synthesized filters \(H_2(s)\) and \(H_3(s)\) are defined below:
And the third filter \(H_1(s)\) is defined using \eqref{eq:H1_complementary_of_H2_H3}.
\begin{minted}[]{matlab}
%% H1 is defined as the complementary filter of H2 and H3
H1 = 1 - H2 - H3;
\end{minted}
The bode plots of the three obtained complementary filters are shown in Figure \ref{fig:three_complementary_filters_results}.
\begin{figure}[htbp]
\centering
\includegraphics[scale=1]{figs/three_complementary_filters_results.png}
\caption{\label{fig:three_complementary_filters_results}The three complementary filters obtained after \(\mathcal{H}_\infty\) synthesis}
\end{figure}
\chapter{Matlab Scripts}
\label{sec:matlab_scripts}
\section{\texttt{1\_synthesis\_complementary\_filters.m}}
\label{sec:1_synthesis_complementary_filters}
This scripts corresponds to section 3 of \cite{dehaeze21_new_method_desig_compl_filter}.
\begin{minted}[linenos,firstnumber=1]{matlab}
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
%% Initialize Frequency Vector
freqs = logspace(-1, 3, 1000);
%% Weighting Function Design
% Parameters
n = 3; w0 = 2*pi*10; G0 = 1e-3; G1 = 1e1; Gc = 2;
% Formulas
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;
%% Magnitude of the weighting function with parameters
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(W, freqs, 'Hz'))), 'k-');
plot([1e-3 1e0], [G0 G0], 'k--', 'LineWidth', 1)
text(1e0, G0, '$\quad G_0$')
plot([1e1 1e3], [G1 G1], 'k--', 'LineWidth', 1)
text(1e1,G1,'$G_{\infty}\quad$','HorizontalAlignment', 'right')
plot([w0/2/pi w0/2/pi], [1 2*Gc], 'k--', 'LineWidth', 1)
text(w0/2/pi,1,'$\omega_c$','VerticalAlignment', 'top', 'HorizontalAlignment', 'center')
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')
text(w0/5/pi/2, abs(evalfr(W, j*w0/5)), 'Slope: $n \quad$', 'HorizontalAlignment', 'right')
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]);
yticks([1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1]);
%% Design of the Weighting Functions
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);
%% Plot of the Weighting function magnitude
figure;
tiledlayout(1, 1, 'TileSpacing', 'None', 'Padding', 'None');
ax1 = nexttile();
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]', 'FontSize', 10); ylabel('Magnitude', 'FontSize', 10);
hold off;
xlim([freqs(1), freqs(end)]);
xticks([0.1, 1, 10, 100, 1000]);
ylim([8e-4, 20]);
yticks([1e-3, 1e-2, 1e-1, 1, 1e1]);
yticklabels({'', '$10^{-2}$', '', '$10^0$', ''});
ax1.FontSize = 9;
leg = legend('location', 'south', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
%% Generalized Plant
P = [W1 -W1;
0 W2;
1 0];
%% H-Infinity Synthesis
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
%% Define H1 to be the complementary of H2
H1 = 1 - H2;
%% Bode plot of the complementary filters
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',1)
plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
set(gca,'ColorOrderIndex',2)
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]); ylabel('Magnitude');
ylim([8e-4, 20]);
yticks([1e-3, 1e-2, 1e-1, 1, 1e1]);
yticklabels({'', '$10^{-2}$', '', '$10^0$', ''})
leg = legend('location', 'south', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 18;
% 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;
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
yticks([-180:90:180]);
ylim([-180, 200])
yticklabels({'-180', '', '0', '', '180'})
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
\end{minted}
\section{\texttt{2\_ligo\_complementary\_filters.m}}
\label{sec:2_ligo_complementary_filters}
This scripts corresponds to section 4 of \cite{dehaeze21_new_method_desig_compl_filter}.
\begin{minted}[linenos,firstnumber=1]{matlab}
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
%% Initialize Frequency Vector
freqs = logspace(-3, 0, 1000);
%% Upper bounds for the complementary filters
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]);
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
%% Initialized CVX
cvx_startup;
cvx_solver sedumi;
%% Frequency vectors
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;
%% Filter order
n = 512;
%% Initialization of filter responses
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]))];
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]))];
%% Convex optimization
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
h = y(2:end);
%% Combine the frequency vectors to form the obtained filter
w = [w1 w2 w3 w4];
H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h;
%% Bode plot of the obtained complementary filters
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1)
plot(w, abs(1-H), '-', 'DisplayName', '$L_1$');
plot([0.1, 10], [0.045, 0.045], 'k:', 'DisplayName', 'Spec. on $L_1$');
set(gca,'ColorOrderIndex',2)
plot(w, abs(H), '-', 'DisplayName', '$H_1$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]); ylabel('Magnitude');
hold off;
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
% Phase
ax2 = nexttile;
hold on;
plot(w, 180/pi*unwrap(angle(1-H)), '-');
plot(w, 180/pi*unwrap(angle(H)), '-');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-360:180:180]); ylim([-380, 200]);
linkaxes([ax1,ax2],'x');
xlim([1e-3, 1]);
%% Design of the weight for the high pass filter
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;
wH = minreal(ss(wH));
%% Design of the weight for the low pass filter
n = 20; % Filter order
Rp = 1; % Peak to peak passband ripple
Wp = 2*pi*0.102; % Edge frequency
% Chebyshev Type I filter design
[b,a] = cheby1(n, Rp, Wp, 'high', 's');
wL = 0.04*tf(a, b);
wL = 1/wL;
wL = minreal(ss(wL));
%% Magnitude of the designed Weights and initial specifications
figure;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(inv(wL), freqs, 'Hz'))), '-', 'DisplayName', '$|W_L|^{-1}$');
plot([0.1, 10], [0.045, 0.045], 'k:', 'DisplayName', 'Spec. on $L_1$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(inv(wH), freqs, 'Hz'))), '-', 'DisplayName', '$|W_H|^{-1}$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
%% Generalized plant for the H-infinity Synthesis
P = [0 wL;
wH -wH;
1 0];
%% Standard H-Infinity synthesis
[Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
%% High pass filter as the complementary of the low pass filter
Hh = 1 - Hl;
%% Minimum realization of the filters
Hh = minreal(Hh);
Hl = minreal(Hl);
%% Bode plot of the obtained filters and comparison with the upper bounds
figure;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'DisplayName', '$L_1^\prime$');
plot([0.1, 10], [0.045, 0.045], 'k:', 'DisplayName', 'Spec. on $L_1$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$H_1^\prime$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
%% Comparison of the complementary filters obtained with H-infinity and with CVX
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', ...
'DisplayName', '$L_1(s)$ - $\mathcal{H}_\infty$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_1(s)$ - $\mathcal{H}_\infty$');
set(gca,'ColorOrderIndex',1);
plot(w, abs(1-H), '--', ...
'DisplayName', '$L_1(s)$ - FIR');
set(gca,'ColorOrderIndex',2);
plot(w, abs(H), '--', ...
'DisplayName', '$H_1(s)$ - FIR');
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(Hl, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hh, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',1);
plot(w, 180/pi*unwrap(angle(1-H)), '--');
set(gca,'ColorOrderIndex',2);
plot(w, 180/pi*unwrap(angle(H)), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-360:180:180]); ylim([-380, 200]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
\end{minted}
\section{\texttt{3\_closed\_loop\_complementary\_filters.m}}
\label{sec:3_closed_loop_complementary_filters}
This scripts corresponds to section 5.1 of \cite{dehaeze21_new_method_desig_compl_filter}.
\begin{minted}[linenos,firstnumber=1]{matlab}
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
%% Initialize Frequency Vector
freqs = logspace(-1, 3, 1000);
%% Design of the Weighting Functions
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);
%% Generalized plant for "closed-loop" complementary filter synthesis
P = [ W1 0 1;
-W1 W2 -1];
%% Standard H-Infinity Synthesis
[L, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
%% Complementary filters
H1 = inv(1 + L);
H2 = 1 - H1;
%% Bode plot of the obtained Complementary filters with upper-bounds
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, 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',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$');
plot(freqs, abs(squeeze(freqresp(L, freqs, 'Hz'))), 'k--', 'DisplayName', '$|L|$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]); ylabel('Magnitude');
ylim([1e-3, 1e3]);
yticks([1e-3, 1e-2, 1e-1, 1, 1e1, 1e2, 1e3]);
yticklabels({'', '$10^{-2}$', '', '$10^0$', '', '$10^2$', ''});
leg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 3);
leg.ItemTokenSize(1) = 18;
% 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;
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
yticks([-180:90:180]);
ylim([-180, 200])
yticklabels({'-180', '', '0', '', '180'})
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
\end{minted}
\section{\texttt{4\_three\_complementary\_filters.m}}
\label{sec:4_three_complementary_filters}
This scripts corresponds to section 5.2 of \cite{dehaeze21_new_method_desig_compl_filter}.
\begin{minted}[linenos,firstnumber=1]{matlab}
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
freqs = logspace(-2, 3, 1000);
%% Design of the Weighting Functions
W1 = generateWF('n', 2, 'w0', 2*pi*1, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.5);
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;
W3 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.5);
%% Inverse magnitude of the weighting functions
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;
xlim([freqs(1), freqs(end)]); ylim([2e-4, 1.3e1])
leg = legend('location', 'northeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
%% Generalized plant for the synthesis of 3 complementary filters
P = [W1 -W1 -W1;
0 W2 0 ;
0 0 W3;
1 0 0];
%% Standard H-Infinity Synthesis
[H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
%% Synthesized H2 and H3 filters
H2 = tf(H(1));
H3 = tf(H(2));
%% H1 is defined as the complementary filter of H2 and H3
H1 = 1 - H2 - H3;
%% Bode plot of the obtained complementary filters
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]);
leg = legend('location', 'northeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
% 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([-180:90:180]); ylim([-220, 220]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
\end{minted}
\section{\texttt{generateWF}: Generate Weighting Functions}
\label{sec:generateWF}
This function is used to easily generate weighting functions from classical requirements.
\begin{minted}[]{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 Weighting Function
%% Argument validation
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
% Verification of correct relation between G0, Gc and Ginf
mustBeBetween(args.G0, args.Gc, args.Ginf);
%% Initialize the Laplace variable
s = zpk('s');
%% 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 + ...
(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;
%% 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{minted}
\section{\texttt{generateCF}: Generate Complementary Filters}
\label{sec:generateCF}
This function is used to easily synthesize a set of two complementary filters using the \(\mathcal{H}_\infty\) synthesis.
\begin{minted}[]{matlab}
function [H1, H2] = generateCF(W1, W2, args)
% createWeight -
%
% Syntax: [H1, H2] = generateCF(W1, W2, args)
%
% Inputs:
% - W1 - Weighting Function for H1
% - W2 - Weighting Function for H2
% - args:
% - method - H-Infinity solver ('lmi' or 'ric')
% - display - Display synthesis results ('on' or 'off')
%
% Outputs:
% - H1 - Generated H1 Filter
% - H2 - Generated H2 Filter
%% Argument validation
arguments
W1
W2
args.method char {mustBeMember(args.method,{'lmi', 'ric'})} = 'ric'
args.display char {mustBeMember(args.display,{'on', 'off'})} = 'on'
end
%% The generalized plant is defined
P = [W1 -W1;
0 W2;
1 0];
%% The standard H-infinity synthesis is performed
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', args.method, 'DISPLAY', args.display);
%% H1 is defined as the complementary of H2
H1 = 1 - H2;
\end{minted}
\printbibliography
\end{document}
Reference in New Issue View Git Blame Copy Permalink