2021-02-02 19:11:20 +01:00
% Created 2021-02-02 mar. 19:10
2021-02-02 18:54:35 +01:00
% Intended LaTeX compiler: pdflatex
2021-02-02 19:11:20 +01:00
\documentclass [a4paper, 10pt, DIV=12, parskip=full] { scrreprt}
2021-02-02 18:54:35 +01:00
\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}
2021-02-02 19:11:20 +01:00
\input { preamble.tex}
2021-02-02 18:54:35 +01:00
\addbibresource { ref.bib}
\author { Dehaeze Thomas}
\date { \today }
\title { Sensor Fusion - Test Bench}
\hypersetup {
pdfauthor={ Dehaeze Thomas} ,
pdftitle={ Sensor Fusion - Test Bench} ,
pdfkeywords={ } ,
pdfsubject={ } ,
pdfcreator={ Emacs 27.1 (Org mode 9.5)} ,
pdflang={ English} }
\begin { document}
\maketitle
\setcounter { tocdepth} { 2}
\tableofcontents
In this document, we wish the experimentally validate sensor fusion of inertial sensors.
This document is divided into the following sections:
\begin { itemize}
\item Section \ref { sec:experimental_ setup} : the experimental setup is described
\item Section \ref { sec:first_ identification} : a first identification of the system dynamics is performed
\item Section \ref { sec:integral_ force_ feedback} : the integral force feedback active damping technique is applied on the system
\item Section \ref { sec:optimal_ excitation} : the optimal excitation signal is determine in order to have the best possible system dynamics estimation
\item Section \ref { sec:inertial_ sensor_ dynamics} : the inertial sensor dynamics are experimentally estimated
\item Section \ref { sec:inertial_ sensor_ noise} : the inertial sensor noises are estimated and the \( \mathcal { H } _ 2 \) synthesis of complementary filters is performed in order to yield a super sensor with minimal noise
\item Section \ref { sec:inertial_ sensor_ uncertainty} : the dynamical uncertainty of the inertial sensors is estimated. Then the \( \mathcal { H } _ \infty \) synthesis of complementary filters is performed in order to minimize the super sensor dynamical uncertainty
\item Section \ref { sec:optimal_ sensor_ fusion} : Optimal sensor fusion is performed using the \( \mathcal { H } _ 2 / \mathcal { H } _ \infty \) synthesis
\end { itemize}
2021-02-02 19:11:20 +01:00
\chapter { Experimental Setup}
\label { sec:org2eab8ee}
2021-02-02 18:54:35 +01:00
\label { sec:experimental_ setup}
The goal of this experimental setup is to experimentally merge inertial sensors.
To merge the sensors, optimal and robust complementary filters are designed.
A schematic of the test-bench used is shown in Figure \ref { fig:exp_ setup_ sensor_ fusion} and a picture of it is shown in Figure \ref { fig:test-bench-sensor-fusion-picture} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/exp_ setup_ sensor_ fusion.png}
\caption { \label { fig:exp_ setup_ sensor_ fusion} Schematic of the test-bench}
\end { figure}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/test-bench-sensor-fusion-picture.png}
\caption { \label { fig:test-bench-sensor-fusion-picture} Picture of the test-bench}
\end { figure}
Two inertial sensors are used:
\begin { itemize}
\item An vertical accelerometer \emph { PCB 393B05} (\href { doc/TM-VIB-Seismic\_ Lowres.pdf} { doc} )
\item A vertical geophone \emph { Mark Product L-22}
\end { itemize}
Basic characteristics of both sensors are shown in Tables \ref { tab:393B05_ spec} and \ref { tab:L22_ spec} .
\begin { table} [htbp]
\caption { \label { tab:393B05_ spec} Accelerometer (393B05) Specifications}
\centering
\begin { tabular} { ll}
\textbf { Specification} & \textbf { Value} \\
\hline
Sensitivity & 1.02 [V/(m/s2)]\\
Resonant Frequency & > 2.5 [kHz]\\
Resolution (1 to 10kHz) & 0.00004 [m/s2 rms]\\
\end { tabular}
\end { table}
\begin { table} [htbp]
\caption { \label { tab:L22_ spec} Geophone (L22) Specifications}
\centering
\begin { tabular} { ll}
\textbf { Specification} & \textbf { Value} \\
\hline
Sensitivity & To be measured [V/(m/s)]\\
Resonant Frequency & 2 [Hz]\\
\end { tabular}
\end { table}
The ADC used are the IO131 Speedgoat module (\href { https://www.speedgoat.com/products/io-connectivity-analog-io131} { link} ) with a 16bit resolution over +/- 10V.
The geophone signals are amplified using a DLPVA-100-B-D voltage amplified from Femto (\href { doc/de-dlpva-100-b.pdf} { doc} ).
The force sensor signal is amplified using a Low Noise Voltage Preamplifier from Ametek (\href { doc/model\_ 5113.pdf} { doc} ).
The excitation signal is amplified by a linear amplified from Cedrat (LA75B) with a gain equals to 20 (\href { doc/LA75B.pdf} { doc} ).
Geophone electronics:
\begin { itemize}
\item gain: 10 (20dB)
\item low pass filter: 1.5Hz
\item hifh pass filter: 100kHz (2nd order)
\end { itemize}
Force Sensor electronics:
\begin { itemize}
\item gain: 10 (20dB)
\item low pass filter: 1st order at 3Hz
\item high pass filter: 1st order at 30kHz
\end { itemize}
2021-02-02 19:11:20 +01:00
\chapter { First identification of the system}
\label { sec:org8c87e73}
2021-02-02 18:54:35 +01:00
\label { sec:first_ identification}
In this section, a first identification of each elements of the system is performed.
This include the dynamics from the actuator to the force sensor, interferometer and inertial sensors.
Each of the dynamics is compared with the dynamics identified form a Simscape model.
2021-02-02 19:11:20 +01:00
\section { Load Data}
\label { sec:orged7fae1}
2021-02-02 18:54:35 +01:00
The data is loaded in the Matlab workspace.
\begin { minted} []{ matlab}
id_ ol = load('identification_ noise_ bis.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
\end { minted}
Then, any offset is removed.
\begin { minted} []{ matlab}
id_ ol.d = detrend(id_ ol.d, 0);
id_ ol.acc_ 1 = detrend(id_ ol.acc_ 1, 0);
id_ ol.acc_ 2 = detrend(id_ ol.acc_ 2, 0);
id_ ol.geo_ 1 = detrend(id_ ol.geo_ 1, 0);
id_ ol.geo_ 2 = detrend(id_ ol.geo_ 2, 0);
id_ ol.f_ meas = detrend(id_ ol.f_ meas, 0);
id_ ol.u = detrend(id_ ol.u, 0);
\end { minted}
2021-02-02 19:11:20 +01:00
\section { Excitation Signal}
\label { sec:org8b54688}
2021-02-02 18:54:35 +01:00
The generated voltage used to excite the system is a white noise and can be seen in Figure \ref { fig:excitation_ signal_ first_ identification} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/excitation_ signal_ first_ identification.png}
\caption { \label { fig:excitation_ signal_ first_ identification} Voltage excitation signal}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Identified Plant}
\label { sec:org3c04eac}
2021-02-02 18:54:35 +01:00
The transfer function from the excitation voltage to the mass displacement and to the force sensor stack voltage are identified using the \texttt { tfestimate} command.
\begin { minted} []{ matlab}
Ts = id_ ol.t(2) - id_ ol.t(1);
win = hann(ceil(10/Ts));
\end { minted}
\begin { minted} []{ matlab}
[tf_ fmeas_ est, f] = tfestimate(id_ ol.u, id_ ol.f_ meas, win, [], [], 1/Ts); % [V/V]
[tf_ G_ ol_ est, ~] = tfestimate(id_ ol.u, id_ ol.d, win, [], [], 1/Ts); % [m/V]
\end { minted}
The bode plots of the obtained dynamics are shown in Figures \ref { fig:force_ sensor_ bode_ plot} and \ref { fig:displacement_ sensor_ bode_ plot} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/force_ sensor_ bode_ plot.png}
\caption { \label { fig:force_ sensor_ bode_ plot} Bode plot of the dynamics from excitation voltage to measured force sensor stack voltage}
\end { figure}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/displacement_ sensor_ bode_ plot.png}
\caption { \label { fig:displacement_ sensor_ bode_ plot} Bode plot of the dynamics from excitation voltage to displacement of the mass as measured by the interferometer}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Simscape Model - Comparison}
\label { sec:orgfbeadc5}
2021-02-02 18:54:35 +01:00
A simscape model representing the test-bench has been developed.
The same transfer functions as the one identified using the test-bench can be obtained thanks to the simscape model.
They are compared in Figure \ref { fig:simscape_ comp_ iff_ plant} and \ref { fig:simscape_ comp_ disp_ plant} .
It is shown that there is a good agreement between the model and the experiment.
\begin { minted} []{ matlab}
load('piezo_ amplified_ 3d.mat', 'int_ xyz', 'int_ i', 'n_ xyz', 'n_ i', 'nodes', 'M', 'K');
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/simscape_ comp_ iff_ plant.png}
\caption { \label { fig:simscape_ comp_ iff_ plant} Comparison of the dynamics from excitation voltage to measured force sensor stack voltage - Identified dynamics and Simscape Model}
\end { figure}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/simscape_ comp_ disp_ plant.png}
\caption { \label { fig:simscape_ comp_ disp_ plant} Comparison of the dynamics from excitation voltage to measured mass displacement - Identified dynamics and Simscape Model}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Integral Force Feedback}
\label { sec:org5839beb}
2021-02-02 18:54:35 +01:00
The force sensor stack can be used to damp the system.
This makes the system easier to excite properly without too much amplification near resonances.
This is done thanks to the integral force feedback control architecture.
The force sensor stack signal is integrated (or rather low pass filtered) and fed back to the force sensor stacks.
The low pass filter used as the controller is defined below:
\begin { minted} []{ matlab}
Kiff = 102/(s + 2*pi*2);
\end { minted}
The integral force feedback control strategy is applied to the simscape model as well as to the real test bench.
The damped system is then identified again using a noise excitation.
The data is loaded into Matlab and any offset is removed.
\begin { minted} []{ matlab}
id_ cl = load('identification_ noise_ iff_ bis.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
id_ cl.d = detrend(id_ cl.d, 0);
id_ cl.acc_ 1 = detrend(id_ cl.acc_ 1, 0);
id_ cl.acc_ 2 = detrend(id_ cl.acc_ 2, 0);
id_ cl.geo_ 1 = detrend(id_ cl.geo_ 1, 0);
id_ cl.geo_ 2 = detrend(id_ cl.geo_ 2, 0);
id_ cl.f_ meas = detrend(id_ cl.f_ meas, 0);
id_ cl.u = detrend(id_ cl.u, 0);
\end { minted}
The transfer functions are estimated using \texttt { tfestimate} .
\begin { minted} []{ matlab}
[tf_ G_ cl_ est, ~] = tfestimate(id_ cl.u, id_ cl.d, win, [], [], 1/Ts);
[co_ G_ cl_ est, ~] = mscohere( id_ cl.u, id_ cl.d, win, [], [], 1/Ts);
\end { minted}
The dynamics from driving voltage to the measured displacement are compared both in the open-loop and IFF case, and for the test-bench experimental identification and for the Simscape model in Figure \ref { fig:iff_ ol_ cl_ identified_ simscape_ comp} .
This shows that the Integral Force Feedback architecture effectively damps the first resonance of the system.
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/iff_ ol_ cl_ identified_ simscape_ comp.png}
\caption { \label { fig:iff_ ol_ cl_ identified_ simscape_ comp} Comparison of the open-loop and closed-loop (IFF) dynamics for both the real identification and the Simscape one}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Inertial Sensors}
\label { sec:org8d9fa37}
2021-02-02 18:54:35 +01:00
In order to estimate the dynamics of the inertial sensor (the transfer function from the ``absolute'' displacement to the measured voltage), the following experiment can be performed:
\begin { itemize}
\item The mass is excited such that is relative displacement as measured by the interferometer is much larger that the ground ``absolute'' motion.
\item The transfer function from the measured displacement by the interferometer to the measured voltage generated by the inertial sensors can be estimated.
\end { itemize}
The first point is quite important in order to have a good coherence between the interferometer measurement and the inertial sensor measurement.
Here, a first identification is performed were the excitation signal is a white noise.
As usual, the data is loaded and any offset is removed.
\begin { minted} []{ matlab}
id = load('identification_ noise_ opt_ iff.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
id.d = detrend(id.d, 0);
id.acc_ 1 = detrend(id.acc_ 1, 0);
id.acc_ 2 = detrend(id.acc_ 2, 0);
id.geo_ 1 = detrend(id.geo_ 1, 0);
id.geo_ 2 = detrend(id.geo_ 2, 0);
id.f_ meas = detrend(id.f_ meas, 0);
\end { minted}
Then the transfer functions from the measured displacement by the interferometer to the generated voltage of the inertial sensors are computed..
\begin { minted} []{ matlab}
Ts = id.t(2) - id.t(1);
win = hann(ceil(10/Ts));
\end { minted}
\begin { minted} []{ matlab}
[tf_ acc1_ est, f] = tfestimate(id.d, id.acc_ 1, win, [], [], 1/Ts);
[co_ acc1_ est, ~] = mscohere( id.d, id.acc_ 1, win, [], [], 1/Ts);
[tf_ acc2_ est, ~] = tfestimate(id.d, id.acc_ 2, win, [], [], 1/Ts);
[co_ acc2_ est, ~] = mscohere( id.d, id.acc_ 2, win, [], [], 1/Ts);
[tf_ geo1_ est, ~] = tfestimate(id.d, id.geo_ 1, win, [], [], 1/Ts);
[co_ geo1_ est, ~] = mscohere( id.d, id.geo_ 1, win, [], [], 1/Ts);
[tf_ geo2_ est, ~] = tfestimate(id.d, id.geo_ 2, win, [], [], 1/Ts);
[co_ geo2_ est, ~] = mscohere( id.d, id.geo_ 2, win, [], [], 1/Ts);
\end { minted}
The same transfer functions are estimated using the Simscape model.
The obtained dynamics of the accelerometer are compared in Figure \ref { fig:comp_ dynamics_ accelerometer} while the one of the geophones are compared in Figure \ref { fig:comp_ dynamics_ geophone} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/comp_ dynamics_ accelerometer.png}
\caption { \label { fig:comp_ dynamics_ accelerometer} Comparison of the measured accelerometer dynamics}
\end { figure}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/comp_ dynamics_ geophone.png}
\caption { \label { fig:comp_ dynamics_ geophone} Comparison of the measured geophone dynamics}
\end { figure}
2021-02-02 19:11:20 +01:00
\chapter { Optimal IFF Development}
\label { sec:org0fdc9f6}
2021-02-02 18:54:35 +01:00
\label { sec:integral_ force_ feedback}
In this section, a proper identification of the transfer function from the force actuator to the force sensor is performed.
Then, an optimal IFF controller is developed and applied experimentally.
The damped system is identified to verified the effectiveness of the added method.
2021-02-02 19:11:20 +01:00
\section { Load Data}
\label { sec:org478ca36}
2021-02-02 18:54:35 +01:00
The experimental data is loaded and any offset is removed.
\begin { minted} []{ matlab}
id_ ol = load('identification_ noise_ bis.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
\end { minted}
\begin { minted} []{ matlab}
id_ ol.d = detrend(id_ ol.d, 0);
id_ ol.acc_ 1 = detrend(id_ ol.acc_ 1, 0);
id_ ol.acc_ 2 = detrend(id_ ol.acc_ 2, 0);
id_ ol.geo_ 1 = detrend(id_ ol.geo_ 1, 0);
id_ ol.geo_ 2 = detrend(id_ ol.geo_ 2, 0);
id_ ol.f_ meas = detrend(id_ ol.f_ meas, 0);
id_ ol.u = detrend(id_ ol.u, 0);
\end { minted}
2021-02-02 19:11:20 +01:00
\section { Experimental Data}
\label { sec:org705ead9}
2021-02-02 18:54:35 +01:00
The transfer function from force actuator to force sensors is estimated.
The coherence shown in Figure \ref { fig:iff_ identification_ coh} shows that the excitation signal is good enough.
\begin { minted} []{ matlab}
Ts = id_ ol.t(2) - id_ ol.t(1);
win = hann(ceil(10/Ts));
\end { minted}
\begin { minted} []{ matlab}
[tf_ fmeas_ est, f] = tfestimate(id_ ol.u, id_ ol.f_ meas, win, [], [], 1/Ts); % [V/m]
[co_ fmeas_ est, ~] = mscohere( id_ ol.u, id_ ol.f_ meas, win, [], [], 1/Ts);
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/iff_ identification_ coh.png}
\caption { \label { fig:iff_ identification_ coh} Coherence for the identification of the IFF plant}
\end { figure}
The obtained dynamics is shown in Figure \ref { fig:iff_ identification_ bode_ plot} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/iff_ identification_ bode_ plot.png}
\caption { \label { fig:iff_ identification_ bode_ plot} Bode plot of the identified IFF plant}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Model of the IFF Plant}
\label { sec:orge2595c0}
2021-02-02 18:54:35 +01:00
In order to plot the root locus for the IFF control strategy, a model of the identified plant is developed.
It consists of several poles and zeros are shown below.
\begin { minted} []{ matlab}
wz = 2*pi*102;
xi_ z = 0.01;
wp = 2*pi*239.4;
xi_ p = 0.015;
Giff = 2.2*(s^ 2 + 2*xi_ z*s*wz + wz^ 2)/(s^ 2 + 2*xi_ p*s*wp + wp^ 2) * ... % Dynamics
10*(s/3/pi/(1 + s/3/pi)) * ... % Low pass filter and gain of the voltage amplifier
exp(-Ts*s); % Time delay induced by ADC/DAC
\end { minted}
The comparison of the identified dynamics and the developed model is done in Figure \ref { fig:iff_ plant_ model} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/iff_ plant_ model.png}
\caption { \label { fig:iff_ plant_ model} IFF Plant + Model}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Root Locus and optimal Controller}
\label { sec:org3cf82e8}
2021-02-02 18:54:35 +01:00
Now, the root locus for the Integral Force Feedback strategy is computed and shown in Figure \ref { fig:iff_ root_ locus} .
Note that the controller used is not a pure integrator but rather a first order low pass filter with a cut-off frequency set at 2Hz.
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/iff_ root_ locus.png}
\caption { \label { fig:iff_ root_ locus} Root Locus for the IFF control}
\end { figure}
The controller that yield maximum damping (shown by the red cross in Figure \ref { fig:iff_ root_ locus} ) is:
\begin { minted} []{ matlab}
Kiff_ opt = 102/(s + 2*pi*2);
\end { minted}
2021-02-02 19:11:20 +01:00
\section { Verification of the achievable damping}
\label { sec:org32eb089}
2021-02-02 18:54:35 +01:00
A new identification is performed with the IFF control strategy applied to the system.
Data is loaded and offset removed.
\begin { minted} []{ matlab}
id_ cl = load('identification_ noise_ iff_ bis.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
\end { minted}
\begin { minted} []{ matlab}
id_ cl.d = detrend(id_ cl.d, 0);
id_ cl.acc_ 1 = detrend(id_ cl.acc_ 1, 0);
id_ cl.acc_ 2 = detrend(id_ cl.acc_ 2, 0);
id_ cl.geo_ 1 = detrend(id_ cl.geo_ 1, 0);
id_ cl.geo_ 2 = detrend(id_ cl.geo_ 2, 0);
id_ cl.f_ meas = detrend(id_ cl.f_ meas, 0);
id_ cl.u = detrend(id_ cl.u, 0);
\end { minted}
The open-loop and closed-loop dynamics are estimated.
\begin { minted} []{ matlab}
[tf_ G_ ol_ est, f] = tfestimate(id_ ol.u, id_ ol.d, win, [], [], 1/Ts);
[co_ G_ ol_ est, ~] = mscohere( id_ ol.u, id_ ol.d, win, [], [], 1/Ts);
[tf_ G_ cl_ est, ~] = tfestimate(id_ cl.u, id_ cl.d, win, [], [], 1/Ts);
[co_ G_ cl_ est, ~] = mscohere( id_ cl.u, id_ cl.d, win, [], [], 1/Ts);
\end { minted}
The obtained coherence is shown in Figure \ref { fig:Gd_ identification_ iff_ coherence} and the dynamics in Figure \ref { fig:Gd_ identification_ iff_ bode_ plot} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/Gd_ identification_ iff_ coherence.png}
\caption { \label { fig:Gd_ identification_ iff_ coherence} Coherence for the transfer function from F to d, with and without IFF}
\end { figure}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/Gd_ identification_ iff_ bode_ plot.png}
\caption { \label { fig:Gd_ identification_ iff_ bode_ plot} Coherence for the transfer function from F to d, with and without IFF}
\end { figure}
2021-02-02 19:11:20 +01:00
\chapter { Generate the excitation signal}
\label { sec:org45b7b7f}
2021-02-02 18:54:35 +01:00
\label { sec:optimal_ excitation}
In order to properly estimate the dynamics of the inertial sensor, the excitation signal must be properly chosen.
The requirements on the excitation signal is:
\begin { itemize}
\item General much larger motion that the measured motion during the huddle test
\item Don't damage the actuator
\end { itemize}
To determine the perfect voltage signal to be generated, we need two things:
\begin { itemize}
\item the transfer function from voltage to mass displacement
\item the PSD of the measured motion by the inertial sensors
\item not saturate the sensor signals
\item provide enough signal/noise ratio (good coherence) in the frequency band of interest (\textasciitilde { } 0.5Hz to 3kHz)
\end { itemize}
2021-02-02 19:11:20 +01:00
\section { Transfer function from excitation signal to displacement}
\label { sec:org91cbb90}
2021-02-02 18:54:35 +01:00
Let's first estimate the transfer function from the excitation signal in [V] to the generated displacement in [m] as measured by the inteferometer.
\begin { minted} []{ matlab}
id_ cl = load('identification_ noise_ iff_ bis.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
\end { minted}
\begin { minted} []{ matlab}
Ts = id_ cl.t(2) - id_ cl.t(1);
win = hann(ceil(10/Ts));
\end { minted}
\begin { minted} []{ matlab}
[tf_ G_ cl_ est, f] = tfestimate(id_ cl.u, id_ cl.d, win, [], [], 1/Ts);
[co_ G_ cl_ est, ~] = mscohere( id_ cl.u, id_ cl.d, win, [], [], 1/Ts);
\end { minted}
Approximate transfer function from voltage output to generated displacement when IFF is used, in [m/V].
\begin { minted} []{ matlab}
G_ d_ est = -5e-6*(2*pi*230)^ 2/(s^ 2 + 2*0.3*2*pi*240*s + (2*pi*240)^ 2);
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/Gd_ plant_ estimation.png}
\caption { \label { fig:Gd_ plant_ estimation} Estimation of the transfer function from the excitation signal to the generated displacement}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Motion measured during Huddle test}
\label { sec:org75df0bc}
2021-02-02 18:54:35 +01:00
We now compute the PSD of the measured motion by the inertial sensors during the huddle test.
\begin { minted} []{ matlab}
ht = load('huddle_ test.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
ht.d = detrend(ht.d, 0);
ht.acc_ 1 = detrend(ht.acc_ 1, 0);
ht.acc_ 2 = detrend(ht.acc_ 2, 0);
ht.geo_ 1 = detrend(ht.geo_ 1, 0);
ht.geo_ 2 = detrend(ht.geo_ 2, 0);
\end { minted}
\begin { minted} []{ matlab}
[p_ d, f] = pwelch(ht.d, win, [], [], 1/Ts);
[p_ acc1, ~] = pwelch(ht.acc_ 1, win, [], [], 1/Ts);
[p_ acc2, ~] = pwelch(ht.acc_ 2, win, [], [], 1/Ts);
[p_ geo1, ~] = pwelch(ht.geo_ 1, win, [], [], 1/Ts);
[p_ geo2, ~] = pwelch(ht.geo_ 2, win, [], [], 1/Ts);
\end { minted}
Using an estimated model of the sensor dynamics from the documentation of the sensors, we can compute the ASD of the motion in \( m / \sqrt { Hz } \) measured by the sensors.
\begin { minted} []{ matlab}
G_ acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
G_ geo = -120*s^ 2/(s^ 2 + 2*0.7*2*pi*2*s + (2*pi*2)^ 2); % [V/(m/s)]
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/huddle_ test_ psd_ motion.png}
\caption { \label { fig:huddle_ test_ psd_ motion} ASD of the motion measured by the sensors}
\end { figure}
From the ASD of the motion measured by the sensors, we can create an excitation signal that will generate much motion motion that the motion under no excitation.
We create \texttt { G\_ exc} that corresponds to the wanted generated motion.
\begin { minted} []{ matlab}
G_ exc = 0.2e-6/(1 + s/2/pi/2)/(1 + s/2/pi/50);
\end { minted}
And we create a time domain signal \texttt { y\_ d} that have the spectral density described by \texttt { G\_ exc} .
\begin { minted} []{ matlab}
Fs = 1/Ts;
t = 0:Ts:180; % Time Vector [s]
u = sqrt(Fs/2)*randn(length(t), 1); % Signal with an ASD equal to one
y_ d = lsim(G_ exc, u, t);
\end { minted}
\begin { minted} []{ matlab}
[pxx, ~] = pwelch(y_ d, win, 0, [], Fs);
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/comp_ huddle_ test_ excited_ motion_ psd.png}
\caption { \label { fig:comp_ huddle_ test_ excited_ motion_ psd} Comparison of the ASD of the motion during Huddle and the wanted generated motion}
\end { figure}
We can now generate the voltage signal that will generate the wanted motion.
\begin { minted} []{ matlab}
y_ v = lsim(G_ exc * ... % from unit PSD to shaped PSD
(1 + s/2/pi/50) * ... % Inverse of pre-filter included in the Simulink file
1/G_ d_ est * ... % Wanted displacement => required voltage
1/(1 + s/2/pi/5e3), ... % Add some high frequency filtering
u, t);
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/optimal_ exc_ signal_ time.png}
\caption { \label { fig:optimal_ exc_ signal_ time} Generated excitation signal}
\end { figure}
2021-02-02 19:11:20 +01:00
\chapter { Identification of the Inertial Sensors Dynamics}
\label { sec:orgcc056c7}
2021-02-02 18:54:35 +01:00
\label { sec:inertial_ sensor_ dynamics}
Using the excitation signal generated in Section \ref { sec:optimal_ excitation} , the dynamics of the inertial sensors are identified.
2021-02-02 19:11:20 +01:00
\section { Load Data}
\label { sec:orge1aa902}
2021-02-02 18:54:35 +01:00
Both the measurement data during the identification test and during an ``huddle test'' are loaded.
\begin { minted} []{ matlab}
id = load('identification_ noise_ opt_ iff.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
ht = load('huddle_ test.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
\end { minted}
\begin { minted} []{ matlab}
ht.d = detrend(ht.d, 0);
ht.acc_ 1 = detrend(ht.acc_ 1, 0);
ht.acc_ 2 = detrend(ht.acc_ 2, 0);
ht.geo_ 1 = detrend(ht.geo_ 1, 0);
ht.geo_ 2 = detrend(ht.geo_ 2, 0);
ht.f_ meas = detrend(ht.f_ meas, 0);
\end { minted}
\begin { minted} []{ matlab}
id.d = detrend(id.d, 0);
id.acc_ 1 = detrend(id.acc_ 1, 0);
id.acc_ 2 = detrend(id.acc_ 2, 0);
id.geo_ 1 = detrend(id.geo_ 1, 0);
id.geo_ 2 = detrend(id.geo_ 2, 0);
id.f_ meas = detrend(id.f_ meas, 0);
\end { minted}
2021-02-02 19:11:20 +01:00
\section { Compare PSD during Huddle and and during identification}
\label { sec:org5e81ca2}
2021-02-02 18:54:35 +01:00
The Power Spectral Density of the measured motion during the huddle test and during the identification test are compared in Figures \ref { fig:comp_ psd_ huddle_ test_ identification_ acc} and \ref { fig:comp_ psd_ huddle_ test_ identification_ geo} .
\begin { minted} []{ matlab}
Ts = ht.t(2) - ht.t(1);
win = hann(ceil(10/Ts));
\end { minted}
\begin { minted} []{ matlab}
[p_ id_ d, f] = pwelch(id.d, win, [], [], 1/Ts);
[p_ id_ acc1, ~] = pwelch(id.acc_ 1, win, [], [], 1/Ts);
[p_ id_ acc2, ~] = pwelch(id.acc_ 2, win, [], [], 1/Ts);
[p_ id_ geo1, ~] = pwelch(id.geo_ 1, win, [], [], 1/Ts);
[p_ id_ geo2, ~] = pwelch(id.geo_ 2, win, [], [], 1/Ts);
\end { minted}
\begin { minted} []{ matlab}
[p_ ht_ d, ~] = pwelch(ht.d, win, [], [], 1/Ts);
[p_ ht_ acc1, ~] = pwelch(ht.acc_ 1, win, [], [], 1/Ts);
[p_ ht_ acc2, ~] = pwelch(ht.acc_ 2, win, [], [], 1/Ts);
[p_ ht_ geo1, ~] = pwelch(ht.geo_ 1, win, [], [], 1/Ts);
[p_ ht_ geo2, ~] = pwelch(ht.geo_ 2, win, [], [], 1/Ts);
[p_ ht_ fmeas, ~] = pwelch(ht.f_ meas, win, [], [], 1/Ts);
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/comp_ psd_ huddle_ test_ identification_ acc.png}
\caption { \label { fig:comp_ psd_ huddle_ test_ identification_ acc} Comparison of the PSD of the measured motion during the Huddle test and during the identification}
\end { figure}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/comp_ psd_ huddle_ test_ identification_ geo.png}
\caption { \label { fig:comp_ psd_ huddle_ test_ identification_ geo} Comparison of the PSD of the measured motion during the Huddle test and during the identification}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Compute transfer functions}
\label { sec:org8059f91}
2021-02-02 18:54:35 +01:00
The transfer functions from the motion as measured by the interferometer (and that should represent the absolute motion of the mass) to the inertial sensors are estimated:
\begin { minted} []{ matlab}
[tf_ acc1_ est, f] = tfestimate(id.d, id.acc_ 1, win, [], [], 1/Ts);
[co_ acc1_ est, ~] = mscohere( id.d, id.acc_ 1, win, [], [], 1/Ts);
[tf_ acc2_ est, ~] = tfestimate(id.d, id.acc_ 2, win, [], [], 1/Ts);
[co_ acc2_ est, ~] = mscohere( id.d, id.acc_ 2, win, [], [], 1/Ts);
[tf_ geo1_ est, ~] = tfestimate(id.d, id.geo_ 1, win, [], [], 1/Ts);
[co_ geo1_ est, ~] = mscohere( id.d, id.geo_ 1, win, [], [], 1/Ts);
[tf_ geo2_ est, ~] = tfestimate(id.d, id.geo_ 2, win, [], [], 1/Ts);
[co_ geo2_ est, ~] = mscohere( id.d, id.geo_ 2, win, [], [], 1/Ts);
\end { minted}
The obtained coherence are shown in Figure \ref { fig:id_ sensor_ dynamics_ coherence} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/id_ sensor_ dynamics_ coherence.png}
\caption { \label { fig:id_ sensor_ dynamics_ coherence} Coherence for the estimation of the sensor dynamics}
\end { figure}
We also make a simplified model of the inertial sensors to be compared with the identified dynamics.
\begin { minted} []{ matlab}
G_ acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
G_ geo = -1200*s^ 2/(s^ 2 + 2*0.7*2*pi*2*s + (2*pi*2)^ 2); % [[V/(m/s)]
\end { minted}
The model and identified dynamics show good agreement (Figures \ref { fig:id_ sensor_ dynamics_ accelerometers} and \ref { fig:id_ sensor_ dynamics_ geophones} .)
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/id_ sensor_ dynamics_ accelerometers.png}
\caption { \label { fig:id_ sensor_ dynamics_ accelerometers} Identified dynamics of the accelerometers}
\end { figure}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/id_ sensor_ dynamics_ geophones.png}
\caption { \label { fig:id_ sensor_ dynamics_ geophones} Identified dynamics of the geophones}
\end { figure}
2021-02-02 19:11:20 +01:00
\chapter { Inertial Sensor Noise and the \( \mathcal { H } _ 2 \) Synthesis of complementary filters}
\label { sec:orgbc08fc3}
2021-02-02 18:54:35 +01:00
\label { sec:inertial_ sensor_ noise}
In this section, the noise of the inertial sensors (geophones and accelerometers) is estimated.
2021-02-02 19:11:20 +01:00
\section { Load Data}
\label { sec:org967725e}
2021-02-02 18:54:35 +01:00
As before, the identification data is loaded and any offset if removed.
\begin { minted} []{ matlab}
id = load('identification_ noise_ opt_ iff.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
\end { minted}
\begin { minted} []{ matlab}
id.d = detrend(id.d, 0);
id.acc_ 1 = detrend(id.acc_ 1, 0);
id.acc_ 2 = detrend(id.acc_ 2, 0);
id.geo_ 1 = detrend(id.geo_ 1, 0);
id.geo_ 2 = detrend(id.geo_ 2, 0);
id.f_ meas = detrend(id.f_ meas, 0);
\end { minted}
2021-02-02 19:11:20 +01:00
\section { ASD of the Measured displacement}
\label { sec:org46e2c6e}
2021-02-02 18:54:35 +01:00
The Power Spectral Density of the displacement as measured by the interferometer and the inertial sensors is computed.
\begin { minted} []{ matlab}
Ts = id.t(2) - id.t(1);
win = hann(ceil(10/Ts));
\end { minted}
\begin { minted} []{ matlab}
[p_ id_ d, f] = pwelch(id.d, win, [], [], 1/Ts);
[p_ id_ acc1, ~] = pwelch(id.acc_ 1, win, [], [], 1/Ts);
[p_ id_ acc2, ~] = pwelch(id.acc_ 2, win, [], [], 1/Ts);
[p_ id_ geo1, ~] = pwelch(id.geo_ 1, win, [], [], 1/Ts);
[p_ id_ geo2, ~] = pwelch(id.geo_ 2, win, [], [], 1/Ts);
\end { minted}
Let's use a model of the accelerometer and geophone to compute the motion from the measured voltage.
\begin { minted} []{ matlab}
G_ acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
G_ geo = -1200*s^ 2/(s^ 2 + 2*0.7*2*pi*2*s + (2*pi*2)^ 2); % [[V/(m/s)]
\end { minted}
The obtained ASD in \( m / \sqrt { Hz } \) is shown in Figure \ref { fig:measure_ displacement_ all_ sensors} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/measure_ displacement_ all_ sensors.png}
\caption { \label { fig:measure_ displacement_ all_ sensors} ASD of the measured displacement as measured by all the sensors}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { ASD of the Sensor Noise}
\label { sec:org39f2454}
2021-02-02 18:54:35 +01:00
The noise of a sensor can be estimated using two identical sensors by computing:
\begin { itemize}
\item the Power Spectral Density of the measured motion by the two sensors
\item the Cross Spectral Density between the two sensors (coherence)
\end { itemize}
This technique to estimate the sensor noise is described in \cite { barzilai98_ techn_ measur_ noise_ sensor_ presen} .
The Power Spectral Density of the sensor noise can be estimated using the following equation:
\begin { equation}
|S_ n(\omega )| = |S_ { x_ 1} (\omega )| \Big ( 1 - \gamma _ { x_ 1 x_ 2} (\omega ) \Big )
\end { equation}
with \( S _ { x _ 1 } \) the PSD of one of the sensor and \( \gamma _ { x _ 1 x _ 2 } \) the coherence between the two sensors.
The coherence between the two accelerometers and between the two geophones is computed.
\begin { minted} []{ matlab}
[coh_ acc, ~] = mscohere(id.acc_ 1, id.acc_ 2, win, [], [], 1/Ts);
[coh_ geo, ~] = mscohere(id.geo_ 1, id.geo_ 2, win, [], [], 1/Ts);
\end { minted}
Finally, the Power Spectral Density of the sensors is computed and converted in \( [ m ^ 2 / Hz ] \) .
\begin { minted} []{ matlab}
pN_ acc = p_ id_ acc1.*(1 - coh_ acc) .* ... % [V^2/Hz]
1./abs(squeeze(freqresp(G_ acc*s^ 2, f, 'Hz'))).^ 2; % [(m/V)^2]
pN_ geo = p_ id_ geo1.*(1 - coh_ geo) .* ... % [V^2/Hz]
1./abs(squeeze(freqresp(G_ geo*s, f, 'Hz'))).^ 2; % [(m/V)^2]
\end { minted}
The ASD of obtained noises are compared with the ASD of the measured signals in Figure \ref { fig:noise_ inertial_ sensors_ comparison} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/noise_ inertial_ sensors_ comparison.png}
\caption { \label { fig:noise_ inertial_ sensors_ comparison} Comparison of the computed ASD of the noise of the two inertial sensors}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Noise Model}
\label { sec:orgc2ba23b}
2021-02-02 18:54:35 +01:00
Transfer functions are adjusted in order to fit the ASD of the sensor noises (expressed in \( [ m / s / \sqrt { Hz } ] \) for more easy fitting).
These transfer functions are defined below and compared with the measured ASD in Figure \ref { fig:noise_ models_ velocity} .
\begin { minted} []{ matlab}
N_ acc = 1*(s/(2*pi*2000) + 1)^ 2/(s + 0.1*2*pi)/(s + 1e3*2*pi); % [m/sqrt(Hz)]
N_ geo = 4e-4*(s/(2*pi*200) + 1)/(s + 1e3*2*pi); % [m/sqrt(Hz)]
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/noise_ models_ velocity.png}
\caption { \label { fig:noise_ models_ velocity} ASD of the velocity noise measured by the sensors and the noise models}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { \( \mathcal { H } _ 2 \) Synthesis of the Complementary Filters}
\label { sec:org79826d9}
2021-02-02 18:54:35 +01:00
We now wish to synthesize two complementary filters to merge the geophone and the accelerometer signal in such a way that the fused signal has the lowest possible RMS noise.
To do so, we use the \( \mathcal { H } _ 2 \) synthesis where the transfer functions representing the noise density of both sensors are used as weights.
The generalized plant used for the synthesis is defined below.
\begin { minted} []{ matlab}
P = [0 N_ acc 1;
N_ geo -N_ acc 0];
\end { minted}
And the \( \mathcal { H } _ 2 \) synthesis is done using the \texttt { h2syn} command.
\begin { minted} []{ matlab}
[H_ geo, ~, gamma] = h2syn(P, 1, 1);
H_ acc = 1 - H_ geo;
\end { minted}
The obtained complementary filters are shown in Figure \ref { fig:complementary_ filters_ velocity_ H2} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/complementary_ filters_ velocity_ H2.png}
\caption { \label { fig:complementary_ filters_ velocity_ H2} Obtained Complementary Filters}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Results}
\label { sec:org1260a50}
2021-02-02 18:54:35 +01:00
Finally, the signals of both sensors are merged using the complementary filters and the super sensor noise is estimated and compared with the individual sensor noises in Figure \ref { fig:super_ sensor_ noise_ asd_ velocity} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/super_ sensor_ noise_ asd_ velocity.png}
\caption { \label { fig:super_ sensor_ noise_ asd_ velocity} ASD of the super sensor noise (velocity)}
\end { figure}
Finally, the Cumulative Power Spectrum is computed and compared in Figure \ref { fig:super_ sensor_ noise_ cas_ velocity} .
\begin { minted} []{ matlab}
[~, i_ 1Hz] = min(abs(f - 1));
\end { minted}
\begin { minted} []{ matlab}
CPS_ acc = 1/pi*flip(-cumtrapz(2*pi*flip(f), flip((pN_ acc.*(2*pi*f)).^ 2)));
CPS_ geo = 1/pi*flip(-cumtrapz(2*pi*flip(f), flip((pN_ geo.*(2*pi*f)).^ 2)));
CPS_ SS = 1/pi*flip(-cumtrapz(2*pi*flip(f), flip((pN_ acc.*(2*pi*f)).^ 2.*abs(squeeze(freqresp(H_ acc, f, 'Hz'))).^ 2 + (pN_ geo.*(2*pi*f)).^ 2.*abs(squeeze(freqresp(H_ geo, f, 'Hz'))).^ 2)));
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/super_ sensor_ noise_ cas_ velocity.png}
\caption { \label { fig:super_ sensor_ noise_ cas_ velocity} Cumulative Power Spectrum of the Sensor Noise (velocity)}
\end { figure}
2021-02-02 19:11:20 +01:00
\chapter { Inertial Sensor Dynamics Uncertainty and the \( \mathcal { H } _ \infty \) Synthesis of complementary filters}
\label { sec:orgf41d7c1}
2021-02-02 18:54:35 +01:00
\label { sec:inertial_ sensor_ uncertainty}
When merging two sensors, it is important to be sure that we correctly know the sensor dynamics near the merging frequency.
Thus, identifying the uncertainty on the sensor dynamics is quite important to perform a robust merging.
2021-02-02 19:11:20 +01:00
\section { Load Data}
\label { sec:org5adcc7b}
2021-02-02 18:54:35 +01:00
Data is loaded and offset is removed.
\begin { minted} []{ matlab}
id = load('identification_ noise_ opt_ iff.mat', 'd', 'acc_ 1', 'acc_ 2', 'geo_ 1', 'geo_ 2', 'f_ meas', 'u', 't');
\end { minted}
\begin { minted} []{ matlab}
id.d = detrend(id.d, 0);
id.acc_ 1 = detrend(id.acc_ 1, 0);
id.acc_ 2 = detrend(id.acc_ 2, 0);
id.geo_ 1 = detrend(id.geo_ 1, 0);
id.geo_ 2 = detrend(id.geo_ 2, 0);
id.f_ meas = detrend(id.f_ meas, 0);
\end { minted}
2021-02-02 19:11:20 +01:00
\section { Compute the dynamics of both sensors}
\label { sec:org74a5ed1}
2021-02-02 18:54:35 +01:00
The dynamics of inertial sensors are estimated (in \( [ V / m ] \) ).
\begin { minted} []{ matlab}
Ts = id.t(2) - id.t(1);
win = hann(ceil(10/Ts));
\end { minted}
\begin { minted} []{ matlab}
[tf_ acc1_ est, f] = tfestimate(id.d, id.acc_ 1, win, [], [], 1/Ts);
[co_ acc1_ est, ~] = mscohere( id.d, id.acc_ 1, win, [], [], 1/Ts);
[tf_ acc2_ est, ~] = tfestimate(id.d, id.acc_ 2, win, [], [], 1/Ts);
[co_ acc2_ est, ~] = mscohere( id.d, id.acc_ 2, win, [], [], 1/Ts);
[tf_ geo1_ est, ~] = tfestimate(id.d, id.geo_ 1, win, [], [], 1/Ts);
[co_ geo1_ est, ~] = mscohere( id.d, id.geo_ 1, win, [], [], 1/Ts);
[tf_ geo2_ est, ~] = tfestimate(id.d, id.geo_ 2, win, [], [], 1/Ts);
[co_ geo2_ est, ~] = mscohere( id.d, id.geo_ 2, win, [], [], 1/Ts);
\end { minted}
The (nominal) models of the inertial sensors from the absolute displacement to the generated voltage are defined below:
\begin { minted} []{ matlab}
G_ acc = 1/(1 + s/2/pi/2000)
G_ geo = -1200*s^ 2/(s^ 2 + 2*0.7*2*pi*2*s + (2*pi*2)^ 2);
\end { minted}
These models are very simplistic models, and we then take into account the un-modelled dynamics with dynamical uncertainty.
2021-02-02 19:11:20 +01:00
\section { Dynamics uncertainty estimation}
\label { sec:org5ba11de}
2021-02-02 18:54:35 +01:00
Weights representing the dynamical uncertainty of the sensors are defined below.
\begin { minted} []{ matlab}
w_ acc = createWeight('n', 2, 'G0', 10, 'G1', 0.2, 'Gc', 1, 'w0', 6*2*pi) * ...
createWeight('n', 2, 'G0', 1, 'G1', 5/0.2, 'Gc', 1/0.2, 'w0', 1300*2*pi);
w_ geo = createWeight('n', 2, 'G0', 0.6, 'G1', 0.2, 'Gc', 0.3, 'w0', 3*2*pi) * ...
createWeight('n', 2, 'G0', 1, 'G1', 10/0.2, 'Gc', 1/0.2, 'w0', 800*2*pi);
\end { minted}
The measured dynamics are compared with the modelled one as well as the modelled uncertainty in Figure \ref { fig:dyn_ uncertainty_ acc} for the accelerometers and in Figure \ref { fig:dyn_ uncertainty_ geo} for the geophones.
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/dyn_ uncertainty_ acc.png}
\caption { \label { fig:dyn_ uncertainty_ acc} Modeled dynamical uncertainty and meaured dynamics of the accelerometers}
\end { figure}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/dyn_ uncertainty_ geo.png}
\caption { \label { fig:dyn_ uncertainty_ geo} Modeled dynamical uncertainty and meaured dynamics of the geophones}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { \( \mathcal { H } _ \infty \) Synthesis of Complementary Filters}
\label { sec:org719855a}
2021-02-02 18:54:35 +01:00
A last weight is now defined that represents the maximum dynamical uncertainty that is allowed for the super sensor.
\begin { minted} []{ matlab}
wu = inv(createWeight('n', 2, 'G0', 0.7, 'G1', 0.3, 'Gc', 0.4, 'w0', 3*2*pi) * ...
createWeight('n', 2, 'G0', 1, 'G1', 6/0.3, 'Gc', 1/0.3, 'w0', 1200*2*pi));
\end { minted}
This dynamical uncertainty is compared with the two sensor uncertainties in Figure \ref { fig:uncertainty_ weight_ and_ sensor_ uncertainties} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/uncertainty_ weight_ and_ sensor_ uncertainties.png}
\caption { \label { fig:uncertainty_ weight_ and_ sensor_ uncertainties} Individual sensor uncertainty (normalized by their dynamics) and the wanted maximum super sensor noise uncertainty}
\end { figure}
The generalized plant used for the synthesis is defined:
\begin { minted} []{ matlab}
P = [wu*w_ acc -wu*w_ acc;
0 wu*w_ geo;
1 0];
\end { minted}
And the \( \mathcal { H } _ \infty \) synthesis using the \texttt { hinfsyn} command is performed.
\begin { minted} []{ matlab}
[H_ geo, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
\end { minted}
\begin { verbatim}
Test bounds: 0.8556 <= gamma <= 1.34
gamma X>=0 Y>=0 rho(XY)<1 p/f
1.071e+00 0.0e+00 0.0e+00 0.000e+00 p
9.571e-01 0.0e+00 0.0e+00 9.436e-16 p
9.049e-01 0.0e+00 0.0e+00 1.084e-15 p
8.799e-01 0.0e+00 0.0e+00 1.191e-16 p
8.677e-01 0.0e+00 0.0e+00 6.905e-15 p
8.616e-01 0.0e+00 0.0e+00 0.000e+00 p
8.586e-01 1.1e-17 0.0e+00 6.917e-16 p
8.571e-01 0.0e+00 0.0e+00 6.991e-17 p
8.564e-01 0.0e+00 0.0e+00 1.492e-16 p
Best performance (actual): 0.8563
\end { verbatim}
The complementary filter is defined as follows:
\begin { minted} []{ matlab}
H_ acc = 1 - H_ geo;
\end { minted}
The bode plot of the obtained complementary filters is shown in Figure
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/h_ infinity_ obtained_ complementary_ filters.png}
\caption { \label { fig:h_ infinity_ obtained_ complementary_ filters} Bode plot of the obtained complementary filters using the \( \mathcal { H } _ \infty \) synthesis}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Obtained Super Sensor Dynamical Uncertainty}
\label { sec:org467a55a}
2021-02-02 18:54:35 +01:00
The obtained super sensor dynamical uncertainty is shown in Figure \ref { fig:super_ sensor_ uncertainty_ h_ infinity} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/super_ sensor_ uncertainty_ h_ infinity.png}
\caption { \label { fig:super_ sensor_ uncertainty_ h_ infinity} Obtained Super sensor dynamics uncertainty}
\end { figure}
2021-02-02 19:11:20 +01:00
\chapter { Optimal and Robust sensor fusion using the \( \mathcal { H } _ 2 / \mathcal { H } _ \infty \) synthesis}
\label { sec:org37f1a9c}
2021-02-02 18:54:35 +01:00
\label { sec:optimal_ sensor_ fusion}
2021-02-02 19:11:20 +01:00
\section { Noise and Dynamical uncertainty weights}
\label { sec:orgbc4a9e6}
2021-02-02 18:54:35 +01:00
\begin { minted} []{ matlab}
N_ acc = (s/(2*pi*2000) + 1)^ 2/(s + 0.1*2*pi)/(s + 1e3*2*pi)/(1 + s/2/pi/1e3); % [m/sqrt(Hz)]
N_ geo = 4e-4*((s + 2*pi)/(2*pi*200) + 1)/(s + 1e3*2*pi)/(1 + s/2/pi/1e3); % [m/sqrt(Hz)]
\end { minted}
\begin { minted} []{ matlab}
w_ acc = createWeight('n', 2, 'G0', 10, 'G1', 0.2, 'Gc', 1, 'w0', 6*2*pi) * ...
createWeight('n', 2, 'G0', 1, 'G1', 5/0.2, 'Gc', 1/0.2, 'w0', 1300*2*pi);
w_ geo = createWeight('n', 2, 'G0', 0.6, 'G1', 0.2, 'Gc', 0.3, 'w0', 3*2*pi) * ...
createWeight('n', 2, 'G0', 1, 'G1', 10/0.2, 'Gc', 1/0.2, 'w0', 800*2*pi);
\end { minted}
\begin { minted} []{ matlab}
wu = inv(createWeight('n', 2, 'G0', 0.7, 'G1', 0.3, 'Gc', 0.4, 'w0', 3*2*pi) * ...
createWeight('n', 2, 'G0', 1, 'G1', 6/0.3, 'Gc', 1/0.3, 'w0', 1200*2*pi));
\end { minted}
\begin { minted} []{ matlab}
P = [wu*w_ acc -wu*w_ acc;
0 wu*w_ geo;
N_ acc -N_ acc;
0 N_ geo;
1 0];
\end { minted}
And the mixed \( \mathcal { H } _ 2 / \mathcal { H } _ \infty \) synthesis is performed.
\begin { minted} []{ matlab}
[H_ geo, ~] = h2hinfsyn(ss(P), 1, 1, 2, [0, 1], 'HINFMAX', 1, 'H2MAX', Inf, 'DKMAX', 100, 'TOL', 1e-3, 'DISPLAY', 'on');
\end { minted}
\begin { minted} []{ matlab}
H_ acc = 1 - H_ geo;
\end { minted}
2021-02-02 19:11:20 +01:00
\section { Obtained Super Sensor Noise}
\label { sec:org074377b}
2021-02-02 18:54:35 +01:00
\begin { minted} []{ matlab}
freqs = logspace(0, 4, 1000);
PSD_ Sgeo = abs(squeeze(freqresp(N_ geo, freqs, 'Hz'))).^ 2;
PSD_ Sacc = abs(squeeze(freqresp(N_ acc, freqs, 'Hz'))).^ 2;
PSD_ Hss = abs(squeeze(freqresp(N_ acc*H_ acc, freqs, 'Hz'))).^ 2 + ...
abs(squeeze(freqresp(N_ geo*H_ geo, freqs, 'Hz'))).^ 2;
\end { minted}
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/psd_ sensors_ htwo_ hinf_ synthesis.png}
\caption { \label { fig:psd_ sensors_ htwo_ hinf_ synthesis} Power Spectral Density of the Super Sensor obtained with the mixed \( \mathcal { H } _ 2 / \mathcal { H } _ \infty \) synthesis}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Obtained Super Sensor Dynamical Uncertainty}
\label { sec:org3d908c9}
2021-02-02 18:54:35 +01:00
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/super_ sensor_ dynamical_ uncertainty_ Htwo_ Hinf.png}
\caption { \label { fig:super_ sensor_ dynamical_ uncertainty_ Htwo_ Hinf} Super sensor dynamical uncertainty (solid curve) when using the mixed \( \mathcal { H } _ 2 / \mathcal { H } _ \infty \) Synthesis}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Experimental Super Sensor Dynamical Uncertainty}
\label { sec:orgfb4a53c}
2021-02-02 18:54:35 +01:00
The super sensor dynamics is shown in Figure \ref { fig:super_ sensor_ optimal_ uncertainty} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/super_ sensor_ optimal_ uncertainty.png}
\caption { \label { fig:super_ sensor_ optimal_ uncertainty} Inertial Sensor dynamics as well as the super sensor dynamics}
\end { figure}
2021-02-02 19:11:20 +01:00
\section { Experimental Super Sensor Noise}
\label { sec:orgaa2a032}
2021-02-02 18:54:35 +01:00
The obtained super sensor noise is shown in Figure \ref { fig:super_ sensor_ optimal_ noise} .
\begin { figure} [htbp]
\centering
\includegraphics [scale=1] { figs/super_ sensor_ optimal_ noise.png}
\caption { \label { fig:super_ sensor_ optimal_ noise} ASD of the super sensor obtained using the \( \mathcal { H } _ 2 / \mathcal { H } _ \infty \) synthesis}
\end { figure}
2021-02-02 19:11:20 +01:00
\chapter { Matlab Functions}
\label { sec:org895c065}
2021-02-02 18:54:35 +01:00
\label { sec:matlab_ functions}
2021-02-02 19:11:20 +01:00
\section { \texttt { createWeight} }
\label { sec:org94e9cbb}
2021-02-02 18:54:35 +01:00
\label { sec:createWeight}
This Matlab function is accessible \href { src/createWeight.m} { here} .
\begin { minted} []{ matlab}
function [W] = createWeight(args)
% createWeight -
%
% Syntax: [in_data] = createWeight(in_data)
%
% Inputs:
% - n - Weight Order
% - G0 - Low frequency Gain
% - G1 - High frequency Gain
% - Gc - Gain of W at frequency w0
% - w0 - Frequency at which |W(j w0)| = Gc
%
% Outputs:
% - W - Generated Weight
arguments
args.n (1,1) double { mustBeInteger, mustBePositive} = 1
args.G0 (1,1) double { mustBeNumeric, mustBePositive} = 0.1
args.G1 (1,1) double { mustBeNumeric, mustBePositive} = 10
args.Gc (1,1) double { mustBeNumeric, mustBePositive} = 1
args.w0 (1,1) double { mustBeNumeric, mustBePositive} = 1
end
mustBeBetween(args.G0, args.Gc, args.G1);
s = tf('s');
W = (((1/args.w0)*sqrt((1-(args.G0/args.Gc)^ (2/args.n))/(1-(args.Gc/args.G1)^ (2/args.n)))*s + (args.G0/args.Gc)^ (1/args.n))/((1/args.G1)^ (1/args.n)*(1/args.w0)*sqrt((1-(args.G0/args.Gc)^ (2/args.n))/(1-(args.Gc/args.G1)^ (2/args.n)))*s + (1/args.Gc)^ (1/args.n)))^ args.n;
end
% 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 G1.';
throwAsCaller(MException(eid,msg))
end
end
\end { minted}
2021-02-02 19:11:20 +01:00
\section { \texttt { plotMagUncertainty} }
\label { sec:orgc1c3060}
2021-02-02 18:54:35 +01:00
\label { sec:plotMagUncertainty}
This Matlab function is accessible \href { src/plotMagUncertainty.m} { here} .
\begin { minted} []{ matlab}
function [p] = plotMagUncertainty(W, freqs, args)
% plotMagUncertainty -
%
% Syntax: [p] = plotMagUncertainty(W, freqs, args)
%
% Inputs:
% - W - Multiplicative Uncertainty Weight
% - freqs - Frequency Vector [Hz]
% - args - Optional Arguments:
% - G
% - color_i
% - opacity
%
% Outputs:
% - p - Plot Handle
arguments
W
freqs double { mustBeNumeric, mustBeNonnegative}
args.G = tf(1)
args.color_ i (1,1) double { mustBeInteger, mustBePositive} = 1
args.opacity (1,1) double { mustBeNumeric, mustBeNonnegative} = 0.3
args.DisplayName char = ''
end
% Get defaults colors
colors = get(groot, 'defaultAxesColorOrder');
p = patch([freqs flip(freqs)], ...
[abs(squeeze(freqresp(args.G, freqs, 'Hz'))).*(1 + abs(squeeze(freqresp(W, freqs, 'Hz')))); ...
flip(abs(squeeze(freqresp(args.G, freqs, 'Hz'))).*max(1 - abs(squeeze(freqresp(W, freqs, 'Hz'))), 1e-6))], 'w', ...
'DisplayName', args.DisplayName);
p.FaceColor = colors(args.color_ i, :);
p.EdgeColor = 'none';
p.FaceAlpha = args.opacity;
end
\end { minted}
2021-02-02 19:11:20 +01:00
\section { \texttt { plotPhaseUncertainty} }
\label { sec:orgf8a75f6}
2021-02-02 18:54:35 +01:00
\label { sec:plotPhaseUncertainty}
This Matlab function is accessible \href { src/plotPhaseUncertainty.m} { here} .
\begin { minted} []{ matlab}
function [p] = plotPhaseUncertainty(W, freqs, args)
% plotPhaseUncertainty -
%
% Syntax: [p] = plotPhaseUncertainty(W, freqs, args)
%
% Inputs:
% - W - Multiplicative Uncertainty Weight
% - freqs - Frequency Vector [Hz]
% - args - Optional Arguments:
% - G
% - color_i
% - opacity
%
% Outputs:
% - p - Plot Handle
arguments
W
freqs double { mustBeNumeric, mustBeNonnegative}
args.G = tf(1)
args.color_ i (1,1) double { mustBeInteger, mustBePositive} = 1
args.opacity (1,1) double { mustBeNumeric, mustBePositive} = 0.3
args.DisplayName char = ''
end
% Get defaults colors
colors = get(groot, 'defaultAxesColorOrder');
% Compute Phase Uncertainty
Dphi = 180/pi*asin(abs(squeeze(freqresp(W, freqs, 'Hz'))));
Dphi(abs(squeeze(freqresp(W, freqs, 'Hz'))) > 1) = 360;
% Compute Plant Phase
G_ ang = 180/pi*angle(squeeze(freqresp(args.G, freqs, 'Hz')));
p = patch([freqs flip(freqs)], [G_ ang+Dphi; flip(G_ ang-Dphi)], 'w', ...
'DisplayName', args.DisplayName);
p.FaceColor = colors(args.color_ i, :);
p.EdgeColor = 'none';
p.FaceAlpha = args.opacity;
end
\end { minted}
\printbibliography
\end { document}
