A new method of designing complementary filters for sensor fusion using the \(\mathcal{H}_\infty\) synthesis - Matlab Computation
+Table of Contents
+-
+
- 1. H-Infinity synthesis of complementary filters + + +
- 2. Design of complementary filters used in the Active Vibration Isolation System at the LIGO + + +
- 3. “Closed-Loop” complementary filters + + +
- 4. Synthesis of three complementary filters + + +
- 5. Functions + + +
+This file is the Matlab file for the paper (Dehaeze, Vermat, and Collette 2021). +
+ ++This document is divided into several sections: +
+-
+
- in section 1, the \(\mathcal{H}_\infty\) synthesis is used for generating two complementary filters +
- in section 4, a method using the \(\mathcal{H}_\infty\) synthesis is proposed to shape three of more complementary filters +
- in section 2, the \(\mathcal{H}_\infty\) synthesis is used and compared with FIR complementary filters used for LIGO +
- in section 3 +
1. H-Infinity synthesis of complementary filters
++The Matlab file corresponding to this section is accessible here. +
+ +1.1. Synthesis Architecture
++We here synthesize two complementary filters using the \(\mathcal{H}_\infty\) synthesis. +The goal is to specify upper bounds on the norms of the two complementary filters \(H_1(s)\) and \(H_2(s)\) while ensuring their complementary property (\(H_1(s) + H_2(s) = 1\)). +
+ ++In order to do so, we use the generalized plant shown on figure 1 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. +
+ + +
+
Figure 1: \(\mathcal{H}_\infty\) synthesis of the complementary filters
++The \(\mathcal{H}_\infty\) synthesis applied on this generalized plant will give a transfer function \(H_2\) (figure 1) 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)\). +
+1.2. Design of Weighting Function
++A formula is proposed to help the design of the weighting functions: +
+\begin{equation} + W(s) = \left( \frac{ + \frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_\infty}\right)^{\frac{2}{n}}}} s + \left(\frac{G_0}{G_c}\right)^{\frac{1}{n}} + }{ + \left(\frac{1}{G_\infty}\right)^{\frac{1}{n}} \frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_\infty}\right)^{\frac{2}{n}}}} s + \left(\frac{1}{G_c}\right)^{\frac{1}{n}} + }\right)^n +\end{equation} + ++The parameters permits to specify: +
+-
+
- the low frequency gain: \(G_0 = lim_{\omega \to 0} |W(j\omega)|\) +
- the high frequency gain: \(G_\infty = lim_{\omega \to \infty} |W(j\omega)|\) +
- the absolute gain at \(\omega_0\): \(G_c = |W(j\omega_0)|\) +
- the absolute slope between high and low frequency: \(n\) +
+The general shape of a weighting function generated using the formula is shown in figure 2. +
+ + +
+
Figure 2: Gain of the Weighting Function formula
+1.3. Example
+%% 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); ++
+
Figure 3: Weights on the complementary filters \(W_1\) and \(W_2\) and the associated performance weights
+1.4. H-Infinity Synthesis
++We define the generalized plant \(P\) on matlab. +
+%% Generalized Plant +P = [W1 -W1; + 0 W2; + 1 0]; ++
+And we do the \(\mathcal{H}_\infty\) synthesis using the hinfsyn command.
+
%% 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'); + + 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 ++ +
+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 4. +
+ +%% Define H1 to be the complementary of H2 +H1 = 1 - H2; ++
+Or one can just used to generateCF Matlab function:
+
[H1, H2] = generateCF(W1, W2); ++
1.5. Obtained Complementary Filters
++The obtained complementary filters are shown on figure 4. +
+ ++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) ++ + +
+
Figure 4: Obtained complementary filters using \(\mathcal{H}_\infty\) synthesis
+2. Design of complementary filters used in the Active Vibration Isolation System at the LIGO
++The Matlab file corresponding to this section is accessible here. +
+ ++Let’s try to design complementary filters that are corresponding to the complementary filters design for the LIGO and described in (Hua 2005). +
+ ++The FIR complementary filters designed in (Hua 2005) are of order 512. +
+2.1. Specifications
++The specifications for the filters are: +
+-
+
- From \(0\) to \(0.008\text{ Hz}\),the magnitude of the filter’s transfer function should be less than or equal to \(8 \times 10^{-3}\) +
- From \(0.008\text{ Hz}\) to \(0.04\text{ Hz}\), it attenuates the input signal proportional to frequency cubed +
- Between \(0.04\text{ Hz}\) and \(0.1\text{ Hz}\), the magnitude of the transfer function should be less than 3 +
- 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 (Hua 2005) have a magnitude of \(0.045\), we will set that as our requirement +
+The specifications are translated in upper bounds of the complementary filters are shown on figure 5. +
+ + +
+
Figure 5: Specification for the LIGO complementary filters
+2.2. FIR Filter
++We here try to implement the FIR complementary filter synthesis as explained in (Hua 2005). +For that, we use the CVX matlab Toolbox. +
+ +
+We setup the CVX toolbox and use the SeDuMi solver.
+
%% Initialized CVX
+cvx_startup;
+cvx_solver sedumi;
+
++We define the frequency vectors on which we will constrain the norm of the FIR filter. +
+%% 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; ++
+We then define the order of the FIR filter. +
+%% 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]))]; ++
+We run the convex optimization. +
+%% 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); ++
+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); ++ +
+Finally, we compute the filter response over the frequency vector defined and the result is shown on figure 6 which is very close to the filters obtain in (Hua 2005). +
+ +%% 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; ++
+
Figure 6: FIR Complementary filters obtain after convex optimization
+2.3. Weights
++We design weights that will be used for the \(\mathcal{H}_\infty\) synthesis of the complementary filters. +These weights will determine the order of the obtained filters. +Here are the requirements on the filters: +
+-
+
- reasonable order +
- to be as close as possible to the specified upper bounds +
- stable minimum phase +
+The bode plot of the weights is shown on figure 7. +
+ + +
+
Figure 7: Weights for the \(\mathcal{H}_\infty\) synthesis
+2.4. H-Infinity Synthesis
++We define the generalized plant as shown on figure 1. +
+%% Generalized plant for the H-infinity Synthesis +P = [0 wL; + wH -wH; + 1 0]; ++
+And we do the \(\mathcal{H}_\infty\) synthesis using the hinfsyn command.
+
%% 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'); +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 ++ +
+The high pass filter is defined as \(H_H = 1 - H_L\). +
+%% High pass filter as the complementary of the low pass filter +Hh = 1 - Hl; ++
+The size of the filters is shown below. +
+ ++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. ++ +
+The bode plot of the obtained filters as shown on figure 8. +
+ + +
+
Figure 8: Obtained complementary filters using the \(\mathcal{H}_\infty\) synthesis
+2.5. Compare FIR and H-Infinity Filters
+ +
+3. “Closed-Loop” complementary filters
++The Matlab file corresponding to this section is accessible here. +
+ +3.1. Using Feedback architecture
+%% 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); ++
+Let’s first synthesize \(H_1(s)\): +
+%% 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; ++
+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) ++ + +
+
4. Synthesis of three complementary filters
++The Matlab file corresponding to this section is accessible here. +
+ +4.1. Theory
++We want: +
+\begin{align*} + & |H_1(j\omega)| < 1/|W_1(j\omega)|, \quad \forall\omega\\ + & |H_2(j\omega)| < 1/|W_2(j\omega)|, \quad \forall\omega\\ + & |H_3(j\omega)| < 1/|W_3(j\omega)|, \quad \forall\omega\\ + & H_1(s) + H_2(s) + H_3(s) = 1 +\end{align*} + ++For that, we use the \(\mathcal{H}_\infty\) synthesis with the architecture shown on figure 11. +
+ + +
+
Figure 11: Generalized architecture for generating 3 complementary filters
++The \(\mathcal{H}_\infty\) objective is: +
+\begin{align*} + & |(1 - H_2(j\omega) - H_3(j\omega)) W_1(j\omega)| < 1, \quad \forall\omega\\ + & |H_2(j\omega) W_2(j\omega)| < 1, \quad \forall\omega\\ + & |H_3(j\omega) W_3(j\omega)| < 1, \quad \forall\omega\\ +\end{align*} + ++And thus if we choose \(H_1 = 1 - H_2 - H_3\) we have solved the problem. +
+4.2. Weights
++First we define the weights. +
+%% 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); ++
+
Figure 12: Three weighting functions used for the \(\mathcal{H}_\infty\) synthesis of the complementary filters
+4.3. H-Infinity Synthesis
+
+Then we create the generalized plant P.
+
%% Generalized plant for the synthesis of 3 complementary filters +P = [W1 -W1 -W1; + 0 W2 0 ; + 0 0 W3; + 1 0 0]; ++
+And we do the \(\mathcal{H}_\infty\) 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'); +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 ++
4.4. Obtained Complementary Filters
++The obtained filters are: +
+%% +H2 = tf(H(1)); +H3 = tf(H(2)); +H1 = 1 - H2 - H3; ++
+
Figure 13: The three complementary filters obtained after \(\mathcal{H}_\infty\) synthesis
++ +
+ +Bibliography
+5. Functions
+5.1. generateWF: Generate Weighting Functions
++This Matlab function is accessible here. +
+Function description
+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 ++
Optional Parameters
+%% 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 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} + +% Verification of correct relation between G0, Gc and Ginf +mustBeBetween(args.G0, args.Gc, args.Ginf); ++
Generate the Weighting function
+%% Initialize the Laplace variable +s = zpk('s'); ++
+The weighting function formula use is: +
+\begin{equation} +\label{orge02c446} + 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} + +%% 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; ++
Verification of the \(G_0\), \(G_c\) and \(G_\infty\) gains
+%% 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 ++
5.2. generateCF: Generate Complementary Filters
++This Matlab function is accessible here. +
+Function description
+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 ++
Optional Parameters
+%% Argument validation +arguments + W1 + W2 + args.method char {mustBeMember(args.method,{'lmi', 'ric'})} = 'ric' + args.display char {mustBeMember(args.display,{'on', 'off'})} = 'on' +end ++
H-Infinity Synthesis
+%% 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; ++
-
-
-
-
-
-
-
-
-
-