\title{A new method of designing complementary filters for sensor fusion using \(\mathcal{H}_\infty\) synthesis}
\begin{document}
\hypersetup{allcolors=teal}
\begin{frontmatter}
\begin{abstract}
Sensor have limited bandwidth and are accurate only in a certain frequency band.
In many applications, the signals of different sensor are fused together in order to either enhance the stability or improve the operational bandwidth of the system.
The sensor signals can be fused using complementary filters.
The tuning of complementary filters is a complex task and is the subject of this paper.
The filters needs to meet design specifications while satisfying the complementary property.
This paper presents a framework to shape the norm of complementary filters using the \(\mathcal{H}_\infty\) norm minimization.
The design specifications are imposed as constraints in the optimization problem by appropriate selection of weighting functions.
The proposed method is quite general and easily extendable to cases where more than two sensors are fused.
Finally, the proposed method is applied to the design of complementary filter design for active vibration isolation of the Laser Interferometer Gravitation-wave Observatory (LIGO).
\item\cite{pascoal99_navig_system_desig_using_time} use LMI to generate complementary filters
\item\cite{hua05_low_ligo,hua04_polyp_fir_compl_filter_contr_system}: FIR + convex optimization
\item Similar to feedback system:
\begin{itemize}
\item\cite{plummer06_optim_compl_filter_their_applic_motion_measur} use H-Infinity to optimize complementary filters (flatten the super sensor noise spectral density)
\item Robustness problems: \cite{zimmermann92_high_bandw_orien_measur_contr} change of phase near the merging frequency
\item Trial and error
\item Although many design methods of complementary filters have been proposed in the literature, no simple method that allows to shape the norm of the complementary filters is available.
A general sensor fusion architecture is shown in Figure \ref{fig:sensor_fusion_overview} where several sensors (here two) are measuring the same physical quantity \(x\).
The two sensors output signals are estimates \(\hat{x}_1\) and \(\hat{x}_2\) of \(x\).
Each of these estimates are then filtered out by complementary filters and combined to form a new estimate \(\hat{x}\).
We further call the overall system from \(x\) to \(\hat{x}\) the ``super sensor''.
The filters \(H_1(s)\) and \(H_2(s)\) are complementary which implies that:
\begin{equation}
\label{eq:comp_filter}
H_1(s) + H_2(s) = 1
\end{equation}
It will soon become clear why the complementary property is important.
\subsection{Sensor Models and Sensor Normalization}
\label{sec:orgfc7a65c}
\label{sec:sensor_models}
In order to study such sensor fusion architecture, a model of the sensor is required.
The sensor model is shown in Figure \ref{fig:sensor_model}.
It consists of a Linear Time Invariant system (LTI) \(G_i(s)\) representing the dynamics of the sensor and an additive noise input \(n_i\) representing its noise.
The model input \(x\) is the measured quantity and its output \(\tilde{x}_i\) is the ``raw'' output of the sensor.
Before filtering the sensor outputs \(\tilde{x}_i\) by the complementary filters, the sensors are usually normalized.
This normalization consists of obtaining an estimate \(\hat{G}_i(s)\) of the sensor dynamics \(G_i(s)\).
The raw output of the sensor \(\tilde{x}_i\) is then passed through the inverse of the sensor dynamics estimate as shown in Figure \ref{fig:sensor_model_calibrated}.
This way, the units of the estimates \(\hat{x}_i\) are equal to the units of the physical quantity \(x\).
The sensor dynamics estimate \(\hat{G}_1(s)\) can be a simple gain or more complex transfer functions.
Let's now combine the two calibrated sensors models (Figure \ref{fig:sensor_model_calibrated}) with the sensor fusion architecture of figure \ref{fig:sensor_fusion_overview}.
The result is shown in Figure \ref{fig:fusion_super_sensor}.
The two sensors are measuring the same physical quantity \(x\) with dynamics \(G_1(s)\) and \(G_2(s)\), and with \emph{uncorrelated} noises \(n_1\) and \(n_2\).
The signals from both calibrated sensors are fed into two complementary filters \(H_1(s)\) and \(H_2(s)\) and then combined to yield an estimate \(\hat{x}\) of \(x\) as shown in Fig. \ref{fig:fusion_super_sensor}.
Thus, this sensor fusion architecture permits to filter the noise of both sensors without introducing any distortion in the physical quantity to be measured.
Let's define the estimation error \(\delta x\) by \eqref{eq:estimate_error}.
As shown in \eqref{eq:noise_filtering_psd}, the Power Spectral Density (PSD) of the estimation error \(\Phi_{\delta x}\) depends both on the norm of the two complementary filters and on the PSD of the noise sources \(\Phi_{n_1}\) and \(\Phi_{n_2}\).
Usually, the two sensors have high noise levels over distinct frequency regions.
In order to lower the noise of the super sensor, the value of the norm \(|H_1|\) has to be lowered when \(\Phi_{n_1}\) is larger than \(\Phi_{n_2}\) and that of \(|H_2|\) lowered when \(\Phi_{n_2}\) is larger than \(\Phi_{n_1}\).
In order to study such imperfection, the sensor dynamical uncertainty is modeled using multiplicative input uncertainty (Figure \ref{fig:sensor_model_uncertainty}), where the nominal model is taken as the estimated model for the normalization \(\hat{G}_i(s)\), \(\Delta_i\) is any stable transfer function satisfying \(|\Delta_i(j\omega)| \le1,\ \forall\omega\), and \(w_i(s)\) is a weight representing the magnitude of the uncertainty.
As the nominal sensor dynamics is taken as the normalized filter, the normalized sensor can be further simplified as shown in Figure \ref{fig:sensor_model_uncertainty_simplified}.
\caption{\label{fig:sensor_fusion_dynamic_uncertainty}Sensor fusion architecture with sensor dynamics uncertainty}
\end{figure}
The super sensor dynamics \eqref{eq:super_sensor_dyn_uncertainty} is no longer equal to \(1\) and now depends on the sensor dynamics uncertainty weights \(w_i(s)\) as well as on the complementary filters \(H_i(s)\).
The uncertainty region of the super sensor can be represented in the complex plane by a circle centered on \(1\) with a radius equal to \(|w_1(j\omega) H_1(j\omega)| + |w_2(j\omega) H_2(j\omega)|\) as shown in Figure \ref{fig:uncertainty_set_super_sensor}.
\caption{\label{fig:uncertainty_set_super_sensor}Uncertainty region of the super sensor dynamics in the complex plane (solid circle). The contribution of both sensors 1 and 2 to the uncertainty are represented respectively by a blue circle and a red circle}
As it is generally desired to limit the maximum phase added by the super sensor, \(H_1(s)\) and \(H_2(s)\) should be designed such that \eqref{eq:max_uncertainty_super_sensor} is satisfied.
The synthesis objective is to shape the norm of two filters \(H_1(s)\) and \(H_2(s)\) while ensuring their complementary property \eqref{eq:comp_filter}.
This is equivalent as to finding stable transfer functions \(H_1(s)\) and \(H_2(s)\) such that conditions \eqref{eq:comp_filter_problem_form} are satisfied.
In order to express this optimization problem as a standard \(\mathcal{H}_\infty\) problem, the architecture shown in Fig. \ref{fig:h_infinity_robust_fusion} is used where the generalized plant \(P\) is described by \eqref{eq:generalized_plant}.
\caption{\label{fig:h_infinity_robust_fusion}Architecture used for \(\mathcal{H}_\infty\) synthesis of complementary filters}
\end{figure}
The \(\mathcal{H}_\infty\) filter design problem is then to find a stable filter \(H_2(s)\) which based on \(v\), generates a signal \(u\) such that the \(\mathcal{H}_\infty\) norm from \(w\) to \([z_1, \ z_2]\) is less than one \eqref{eq:hinf_syn_obj}.
The complementary condition \eqref{eq:hinf_cond_complementarity} is ensured by \eqref{eq:definition_H1}.
The conditions \eqref{eq:hinf_cond_h1} and \eqref{eq:hinf_cond_h2} on the filters shapes are satisfied by \eqref{eq:hinf_problem}.
Therefore, all the conditions \eqref{eq:comp_filter_problem_form} are satisfied using this synthesis method based on \(\mathcal{H}_\infty\) synthesis, and thus it permits to shape complementary filters as desired.
The proper design of the weighting functions is of primary importance for the success of the presented complementary filters \(\mathcal{H}_\infty\) synthesis.
First, only proper, stable and minimum phase transfer functions should be used.
Second, the order of the weights should stay reasonably small in order to reduce the computational costs associated with the solving of the optimization problem and for the physical implementation of the filters (the order of the synthesized filters being equal to the sum of the weighting functions order).
Third, one should not forget the fundamental limitations imposed by the complementary property \eqref{eq:comp_filter}.
This implies for instance that \(|H_1(j\omega)|\) and \(|H_2(j\omega)|\) cannot be made small at the same time.
When designing complementary filters, it is usually desired to specify the slope of the filter, its crossover frequency and its gain at low and high frequency.
To help with the design of the weighting functions such that the above specification can be easily expressed, the following formula is proposed.
\begin{equation}
\label{eq:weight_formula}
W(s) = \left( \frac{
\hfill{}\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 permit to specify:
\begin{itemize}
\item the low frequency gain: \(G_0= lim_{\omega\to0} |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 parameters \(G_0\), \(G_c\) and \(G_\infty\) should either satisfy condition \eqref{eq:cond_formula_1} or \eqref{eq:cond_formula_2}.
\caption{\label{fig:weight_formula}Magnitude of a weighting function generated using the proposed formula \eqref{eq:weight_formula}, \(G_0=1e^{-3}\), \(G_\infty=10\), \(\omega_c =\SI{10}{Hz}\), \(G_c =2\), \(n =3\)}
\end{figure}
\subsection{Validation of the proposed synthesis method}
Let's validate the proposed design method of complementary filters with a simple example where two complementary filters \(H_1(s)\) and \(H_2(s)\) have to be designed such that:
\begin{itemize}
\item the merging frequency is around \(\SI{10}{Hz}\)
\item the slope of \(|H_1(j\omega)|\) is \(-2\) above \(\SI{10}{Hz}\)
\item the slope of \(|H_2(j\omega)|\) is \(+3\) below \(\SI{10}{Hz}\)
\item the gain of both filters is equal to \(10^{-3}\) away from the merging frequency
\end{itemize}
The weighting functions \(W_1(s)\) and \(W_2(s)\) are designed using \eqref{eq:weight_formula}.
The parameters used are summarized in table \ref{tab:weights_params} and the magnitude of the weighting functions is shown in Fig. \ref{fig:hinf_filters_results}.
The bode plots of the obtained complementary filters are shown in Fig. \ref{fig:hinf_filters_results} and their transfer functions in the Laplace domain are given below.
\caption{\label{fig:hinf_filters_results}Frequency response of the weighting functions and complementary filters obtained using \(\mathcal{H}_\infty\) synthesis}
Several complementary filters are used in the active isolation system at the LIGO \cite{hua05_low_ligo,hua04_polyp_fir_compl_filter_contr_system}.
The requirements on those filters are very tight and thus their design is complex.
The approach used in \cite{hua05_low_ligo} for their design is to write the synthesis of complementary FIR filters as a convex optimization problem.
The obtained FIR filters are compliant with the requirements. However they are of very high order so their implementation is quite complex.
The effectiveness of the proposed method is demonstrated by designing complementary filters with the same requirements as the one described in \cite{hua05_low_ligo}.
The specifications for one pair of complementary filters used at the LIGO are summarized below (for further details, refer to \cite{hua04_polyp_fir_compl_filter_contr_system}) and shown in Fig. \ref{fig:ligo_weights}:
\begin{itemize}
\item From \(0\) to \(\SI{0.008}{Hz}\), the magnitude of the filter's transfer function should be less or equal to \(8\times10^{-4}\)
\item Between \(\SI{0.008}{Hz}\) to \(\SI{0.04}{Hz}\), the filter should attenuate the input signal proportional to frequency cubed
\item Between \(\SI{0.04}{Hz}\) to \(\SI{0.1}{Hz}\), the magnitude of the transfer function should be less than \(3\)
\item Above \(\SI{0.1}{Hz}\), the magnitude of the complementary filter should be less than \(0.045\)
The weighting functions should be designed such that their inverse magnitude is as close as possible to the specifications in order to not over-constrain the synthesis problem.
However, the order of each weight should stay reasonably small in order to reduce the computational costs of the optimization problem as well as for the physical implementation of the filters.
A Type I Chebyshev filter of order \(20\) is used as the weighting transfer function \(w_L(s)\) corresponding to the low pass filter.
For the one corresponding to the high pass filter \(w_H(s)\), a \(7^{\text{th}}\) order transfer function is designed.
The magnitudes of the weighting functions are shown in Fig. \ref{fig:ligo_weights}.
\caption{\label{fig:comp_fir_ligo_hinf}Comparison of the FIR filters (solid) designed in \cite{hua05_low_ligo} with the filters obtained with \(\mathcal{H}_\infty\) synthesis (dashed)}
&\left| H_i(j\omega) \right| < \frac{1}{\left| W_i(j\omega) \right|}, \quad\forall\omega,\ i = 1 \dots n \label{eq:hinf_cond_perf_gen}
\end{align}
\end{subequations}
The synthesis method is generalized here for the synthesis of three complementary filters using the architecture shown in Fig. \ref{fig:comp_filter_three_hinf}.
The \(\mathcal{H}_\infty\) synthesis objective applied on \(P(s)\) is to design two stable filters \(H_2(s)\) and \(H_3(s)\) such that the \(\mathcal{H}_\infty\) norm of the transfer function from \(w\) to \([z_1,\ z_2, \ z_3]\) is less than one \eqref{eq:hinf_syn_obj_three}.
\caption{\label{fig:comp_filter_three_hinf}Architecture for \(\mathcal{H}_\infty\) synthesis of three complementary filters}
\end{figure}
By choosing \(H_1(s)\triangleq1- H_2(s)- H_3(s)\), the proposed \(\mathcal{H}_\infty\) synthesis solves the design problem \eqref{eq:hinf_problem_gen}. \par
An example is given to validate the method where three sensors are used in different frequency bands (up to \(\SI{1}{Hz}\), from \(1\) to \(\SI{10}{Hz}\) and above \(\SI{10}{Hz}\) respectively).
Three weighting functions are designed using \eqref{eq:weight_formula} and shown by dashed curves in Fig. \ref{fig:three_complementary_filters_results}.
The bode plots of the obtained complementary filters are shown in Fig. \ref{fig:three_complementary_filters_results}.
\caption{\label{fig:three_complementary_filters_results}Frequency response of the weighting functions and three complementary filters obtained using \(\mathcal{H}_\infty\) synthesis}
This paper has shown how complementary filters can be used to combine multiple sensors in order to obtain a super sensor.
Typical specification on the super sensor noise and on the robustness of the sensor fusion has been shown to be linked to the norm of the complementary filters.
Therefore, a synthesis method that permits the shaping of the complementary filters norms has been proposed and has been successfully applied for the design of complex filters.
Future work will aim at further developing this synthesis method for the robust and optimal synthesis of complementary filters used in sensor fusion.
