chapitre3.tex

% !TEX encoding = UTF-8 Unicode
\chapter{Simulation results} \label{chap-simulation}
%
This chapter describes numerical simulations carried out to demonstrate and evaluate the disturbance models and observer designs described in Chapter \ref{chap-methods}. Section \ref{sec:sim-RODDs} presents a few simulated examples of different \gls{RODD}s. Section \ref{sec:sim-obs-lin} describes experiments to evaluate the multiple-model observers in estimating the states of two simulated systems, a \acrlong{SISO} (\acrshort{SISO}) linear system with one \gls{RODD} step disturbance, and a 2-input, 2-output linear system with two \gls{RODD} step disturbances. Section \ref{sec:sim-ore-SISO} describes an experiment to evaluate the observer performance on the simulated grinding circuit model with one output variable and a step disturbance in the ore feed mix. %\textcolor{red}{Section \ref{sec:sim-ore-mimo-ctrl} describes an experiment to evaluate the observers in a multi-variable feedback control scenario with the simulated grinding circuit model.}


\section{Generating RODD disturbances} \label{sec:sim-RODDs}

\gls{RODD}s are easy to generate by numerical simulation.  Plot (a) in Figure \ref{fig:rodd-sim-plots} shows a random shock sequence defined by \eqref{eq:wpk1} of length 1000 samples generated using a pseudo-random number generator. Plots (b), (c), and (d) show step, ramp, and exponential change disturbances generated with this random shock sequence using the \gls{RODD} models in \eqref{eq:RODD-step}, \eqref{eq:RODD-ramp}, and \eqref{eq:RODD-exp}. Figure \ref{fig:rodd-sim-plot2} shows a combined \gls{RODD} with abruptly changing ramps and steps defined by \eqref{eq:RODD-step-ramp}.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rodd_sim_plots.pdf}
	\caption{Examples of \gls{RODD}s}
	\label{fig:rodd-sim-plots}
\end{figure}
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rodd_sim_plot2.pdf}
	\caption{A \gls{RODD} with steps and ramps}
	\label{fig:rodd-sim-plot2}
\end{figure}


\section{Observer evaluation with linear systems} \label{sec:sim-obs-lin}

\subsection{SISO linear system} \label{sec:sim-obs-lin-1}

To demonstrate state estimation in the presence of \gls{RODD}s, a single \gls{RODD} was simulated at the input to a SISO process represented by a discrete-time linear model, as shown in the functional diagram in Figure \ref{fig:sim-sys-diag-siso}. In addition to the unmeasured \gls{RODD}, $p(k)$, the process has a known input, $u(k)$, and a measured output, $y_M(k)$. The measurements are simulated by adding a random noise, $v(k)$, with zero mean and standard deviation, $\sigma_M$, to the output of the process. At each sample time, the input and measured output are passed to an observer, which calculates estimates of the process states, $\hat{\mathbf{x}}(k|k)$, and an estimate of the process output, $\hat{y}(k|k)$. % TODO: add glossary items for SISO variables.
\nomenclature{$\sigma_M$}{standard deviation of the measurement noise}
\begin{figure}[htp]
	\centering
	\includegraphics[width=11.5cm]{images/sim-sys-diag-siso.pdf}
	\caption{Functional diagram of the simulated \acrshort{SISO} system with observer}
	\label{fig:sim-sys-diag-siso}
\end{figure}

For the purposes of this experiment, the linear model used to represent the process was the stable first order system,
\begin{equation}
	\frac{B(q^{-1})}{A(q^{-1})} = \frac{0.3q^{-1}}{1-0.7q^{-1}}
\end{equation}
with a sampling period, $T_s$, of 0.5.
\nomenclature{$T_s$}{sampling interval in units of time}%

The \gls{RODD} was a step disturbance created by setting
\begin{equation}
	\frac{C(q^{-1})}{D(q^{-1})} = \frac{1}{1-q^{-1}}.
\end{equation}
The random shock, $w_p(k)$, was defined by \eqref{eq:wpik2} with $\varepsilon=0.01$, $\sigma_{w_p}=0.01$, and $b=100$.

The state-space model used to simulate the augmented system was
\begin{equation} \label{eq:sim-sys-siso-ss-aug}
	\begin{split}
	\mathbf{x}_{a}(k+1) & =\left[\begin{array}{cc}
		0.7 & 1 \\
		0 & 1
	\end{array}\right] \mathbf{x}_{a}(k)+\left[\begin{array}{l}
		1 \\
		0
	\end{array}\right] u(k) + \mathbf{w}_{a}(k) \\
	y_M(k) & =\left[\begin{array}{cc}
	0.3 & 0
\end{array}\right] \mathbf{x}_{a}(k) + v(k)
\end{split}
\end{equation}
where
\begin{equation} \label{eq:sim-sys-siso-ss-aug2}
		\mathbf{x}_{a}(k) = \left[\begin{array}{l}
			x_{a,1}(k) \\
			x_{a,2}(k)
		\end{array}\right] = \left[\begin{array}{l}
		x_{1}(k) \\
		p(k)
	\end{array}\right], \mathbf{w}_{a}(k) = \left[\begin{array}{l}
	w(k) \\
	w_{p}(k)
\end{array}\right] .
\end{equation}

Note that with this representation, the second model state $x_{a,2}(k)$ corresponds exactly to the input disturbance $p(k)$.

\subsection{Analysis of sub-optimal estimators} \label{sec:sim-obs-lin-1-SKF-analysis}

To understand and investigate the behaviour of the observers, the system \eqref{eq:sim-sys-siso-ss-aug} was first simulated for 100 sample periods with two pre-determined shocks, no persistent disturbance ($\sigma_{w_p}=0$), and no measurement noise ($\sigma_M=0$). Figure \ref{fig:rod-obs-sim-test-ioplot-SF95} shows the simulation data. The two lower plots show the system inputs, including the random shock signal. The upper two plots show the system states and outputs as well as the estimates of a sub-optimal multiple-model observer using the sequence fusion algorithm described by \cite{robertson_detection_1995} with parameters $N_f=15$, $n_\text{max}=1$, and $d=5$.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_MKF_SF_test_sim_MKF_SF95_ioplot.pdf}
	\caption{Simulation of a \acrshort{SISO} linear system with a \acrshort{RODD} input disturbance}
	\label{fig:rod-obs-sim-test-ioplot-SF95}
\end{figure}

\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{images/rod_MKF_test_sim_MKF_SF95_prob.png}
	\caption{Multi-model observer probability estimates – MKF--SF95}
	\label{fig:rod-obs-sim-test-probs-SF95}
\end{figure}
The two plots in Figure \ref{fig:rod-obs-sim-test-probs-SF95} provide insight into the computations of the observer, which is labelled `MKF--SF95'. The upper plot shows four time series of the trace (i.e. the sum of the elements on the main diagonal) of the merged error covariance matrices, $\mathbf{P}_m(k)$ for $m=1,2,...,n_m$. The diagonal values of $\mathbf{P}_m(k)$ may be interpreted as indicators of the magnitude of the errors of the state estimates. Recall that $\mathbf{P}_m(k)$ is time-varying and influenced by the process noise covariance, $\mathcal{Q}(\gamma_f(k))$, which switches according to the hypothesis sequences \eqref{eq:Pfkp1}. In this simulation, the covariances at time $t=0$ were initialized to the identity matrix ($\mathbf{P}_m(0)=\mathbf{I}_2$). It can be seen that the trace values initially drop rapidly and tend towards zero, however, at certain times, they increase sharply to a peak before dropping again. The sharp peaks are caused by the shock hypotheses ($\gamma_m(k)=1$) of the sequences.

The lower plot shows the conditional likelihood estimates of the four hypotheses given the data up to time $k$ after the merging step \eqref{eq:xmkymk_hat_MKF}. The conditional likelihoods were initialized with equal values ($\operatorname{Pr}\left(\Gamma_m(0) \mid \mathbf{Y}_M(0)\right)=1/n_m$). It can be seen that the likelihoods of each hypothesis settle on a steady-state by $t=5$ and hypothesis 1 dominates the probability distribution until $t=10$. This means that the merged estimates of the observer are determined almost completely by hypothesis 1 during this period. This makes sense since hypothesis 1 represents a sequence containing no shocks and the first shock does not occur until $t=9.5$. A short time after the first shock, the probability density shifts to hypothesis sequence 3 which happens to assume a shock at $t=10$. This is evident in the upper plot where it can be seen that the error covariance of hypothesis 3 increases and peaks at $t=10.5$. The probability switches from hypothesis 1 to 3 by $t=11.5$ and remains there until $t=17$ when it shifts back to hypothesis 1. This final shift can be explained by the fusion horizon, $N_f$, which in this case is 15 sample periods or 7.5 time units. Hypothesis 3 assumes another shock occurs at this time, which is not the case. Therefore, hypothesis 1 becomes the most likely hypothesis after that. The second shock occurs at $t=29$, which is not as favourable since it does not align well with the shocks in any of the hypothesis sequences. Initially, the likelihood switches to hypotheses 4, which assumes a shock at $t=27.5$ but then it switches to hypothesis 2, which assumes a shock at $t=30$. Again, the probability remains high until the next shock is scheduled, at which point hypothesis 1 becomes the most likely again.

\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_MKF_test_sim_MKF_SF95_x_est.pdf}
	\caption{Multi-model observer state estimates – MKF--SF95}
	\label{fig:rod-obs-sim-test-x_est-SF95}
\end{figure}
Figure \ref{fig:rod-obs-sim-test-x_est-SF95} shows the merged state estimates associated with each hypothesis (thin coloured lines), as well as the overall state estimate (thick blue lines) and the true states of the system (dashed lines). These plots reveal the role that each hypothesis plays in the overall estimates. The estimates of hypothesis 3 (purple lines) are the first to respond to the first shock at time $t=9.5$ and the overall estimate closely follows this hypothesis. The other estimates (green, yellow and orange lines) respond at later intervals and do not play a role in the overall estimate. After the second shock occurs, the overall estimate follows that of hypothesis 4 at first, then it switches to 2, however, it lags the true system state and the response of the output estimate is slower and overshoots the true system output for a few time steps. These simulation results demonstrate that the performance of the sequence fusion algorithm is somewhat sensitive to the timing of the true shocks. Nevertheless, the output estimates are quite good and would likely be satisfactory in most process control applications. 

As described in Section \ref{sec:fusion}, Robertson and coworkers proposed an alternative implementation of the sequence fusion algorithm \citep{robertson_method_1998}. There are two differences between this and the previous algorithm. Firstly, shocks are assumed to act over the duration of the detection interval rather than at one sample time within it. Secondly, the variance of the shocks, $b^2\sigma_{w_p}^2$, is reduced according to the interval length \eqref{eq:wpdk}.
%% TODO: add parameters to above? with $f=15$ $m=1$ and $d=5$.

\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{images/rod_MKF_test_sim_MKF_SF98_prob.png}
	\caption{Multi-model observer probability estimates – MKF--SF98}
	\label{fig:rod-obs-sim-test-probs-SF98}
\end{figure}
Figures \ref{fig:rod-obs-sim-test-probs-SF98} and \ref{fig:rod-obs-sim-test-x_est-SF98} show the results of simulating the 1998 version of the observer, labelled `MKF--SF98', with the same simulation data. The effect of the modifications on the error covariance is clear. The traces of the error covariance peak at a lower value and the peaks are sustained for the duration of the detection intervals, unlike in the 1995 version. The switching of the hypothesis probabilities is also different. After the first shock at $t=9.5$ the probability shifts to hypothesis 2 rather than 3. This may be due to the fact that the shock occurred in the middle of the detection interval of hypothesis 2. After the second shock at $t=29$ the hypothesis 4 becomes the most probable.

Figure \ref{fig:rod-obs-sim-test-x_est-SF98} shows the merged state estimates of each hypothesis, the overall estimates, and the true system states for the 1998 version of the algorithm. The estimates responding quickly to the shocks and there is noticeably less overshoot of the estimates in this simulation compared to the 1995 version.   
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_MKF_test_sim_MKF_SF98_x_est.pdf}
	\caption{Multi-model observer state estimates – MKF--SF98}
	\label{fig:rod-obs-sim-test-x_est-SF98}
\end{figure}

\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{images/rod_MKF_test_sim_MKF_SP_prob.png}
	\caption{Multi-model observer probability estimates – MKF--SP}
	\label{fig:rod-obs-sim-test-probs-SP}
\end{figure}
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_MKF_test_sim_MKF_SP_x_est.pdf}
	\caption{Multi-model observer state estimates – MKF--SP}
	\label{fig:rod-obs-sim-test-x_est-SP}
\end{figure}
Figures \ref{fig:rod-obs-sim-test-probs-SP} and \ref{fig:rod-obs-sim-test-x_est-SP} show the corresponding results of simulating the sequence pruning multiple-model observer, labelled `MKF--SP', with $n_h=5$ and $N_\text{min}=2$, as described by \cite{eriksson_classification_1996}. Recall that in this algorithm, hypotheses are adaptively pruned and replaced at every sample time. From these plots it can be seen that the hypotheses include more possible shocks than those modelled by the sequence fusion observers with detection intervals. Also, the switching of hypothesis probabilities is swifter and distinct, and the estimation errors are visibly lower than those of both versions of the sequence fusion observer in these simulations.

To effectively evaluate and compare the performance of these observers, longer simulations with pseudo-random shocks and measurement noise are needed. Figure \ref{fig:rod-obs-sim1-ioplot} shows the first 600 input-output samples from a simulation of the system with a total duration of 2500 (5000 samples). The lower plot shows the input signals, $u(k)$ and $p(k)$, from which it can be seen that two significant random shocks occurred during this period (at $t=98$ and $220.5$). At other times, the value of $p(k)$ is a random walk due to the persistent component of the \gls{RODD} disturbance \eqref{eq:wpk2}. In addition, a step change in $u(k)$ of magnitude 1 occurred at time $t=5$. The upper plot shows the simulated output measurements. The standard deviation of the measurement noise, $\sigma_M$, was $0.1$. No disturbances other than $w_p(k)$ were simulated (i.e. $w_1(k)=0$).
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim1_all_seed_ioplot.pdf}
	\caption{Simulation of linear system with a \acrshort{RODD} input disturbance}
	\label{fig:rod-obs-sim1-ioplot}
\end{figure}

\subsubsection{Tuning of Kalman filters} \label{sec:sim-obs-lin-1-KF-tuning}

Standard Kalman filters were used to develop performance baselines with which to evaluate the sub-optimal multiple-model observers. The question of how to set the gain of a Kalman filter when there is a \gls{RODD} affecting the system is open to interpretation. One naive approach is to tune the Kalman filter using the variance of the persistent component of the disturbance, $\sigma_{w_p}^2$, if it is known. In this case, the Kalman filter, which is labelled `KF1', has a noise covariance matrix
\nomenclature{$\mathbf{Q}_0$}{noise covariance matrix used in observer simulations}%
\begin{equation} \label{eq:sim-sys-siso-KF1-Q}
	\begin{aligned}
		\mathbf{Q}_{\text{KF1}}=\mathbf{Q}_0=\begin{bmatrix}
			\sigma_{x}^2 & 0 \\
			0 &  \sigma_{w_p}^2 \\
		\end{bmatrix}=\begin{bmatrix}
		0.1^2 & 0 \\
		0 & 0.01^2 \\
	\end{bmatrix}
	\end{aligned},
\end{equation}
\nomenclature{$\sigma_{x}^2$}{variance of the errors of the process model state(s)}
where $\sigma_{x}^2$ is a parameter representing the variance of the errors of the process model state, $x(k)$. However, this leads to a filter with a slow response to shocks. The choice of $\sigma_{x}=0.1$ was somewhat arbitrary. In these simulations, the process is simulated by a model that is identical to the model used by the observers, so there is no model error in steady-state. However, in transient periods such as after a shock disturbance, and in real applications, errors between the observer model and the true states are expected.

A second option is to tune a Kalman filter to the variance of the infrequent shocks, $b^2\sigma_{w_p}^2$.  In this case, the Kalman filter, labelled `KF2', has a state disturbance covariance
\nomenclature{$\mathbf{Q}_1$}{noise covariance matrix used in observer simulations}%
\begin{equation} \label{eq:sim-sys-siso-KF2-Q}
	\begin{aligned}
		\mathbf{Q}_{\text{KF2}}=\mathbf{Q}_1=\begin{bmatrix}
			\sigma_{x}^2 & 0 \\
			0 & b^2\sigma_{w_p}^2 \\
		\end{bmatrix}=\begin{bmatrix}
			0.1^2 & 0 \\
			0 & 1^2 \\
		\end{bmatrix}.
	\end{aligned}
\end{equation}
The problem with this filter is that it is very sensitive to the measurement noise, $v(k)$. 

As explained by \cite{andersson_adaptive_1985} and \cite{robertson_detection_1995}, a trade-off is needed. One approach to making this trade-off is to minimize the average estimation errors over a suitably long set of input-output data. In this case, 5000 samples were generated from the simulated system (\ref{fig:sim-sys-diag-siso}) and then a Kalman filter was simulated multiple times, each time with a different tuning. 
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim1_3KF_Q_seed_6.pdf}
	\caption{Tuning of Kalman filter KF3 – \acrshort{SISO} system}
	\label{fig:sim-sys-siso-KF3-tuning}
\end{figure}
Figure \ref{fig:sim-sys-siso-KF3-tuning} shows the results of these simulations, from which it can be seen that the lowest \gls{RMSE}s of the state estimates were achieved when the variance parameter $\sigma_{w_p}^2$ was close to $0.01$. Using this result, a third Kalman filter labelled `KF3' was designed with the state disturbance covariance matrix
\nomenclature{$\mathbf{Q}_{\text{opt}}$}{noise covariance matrix used in the Kalman filter tuned to minimize overall \acrshort{RMSE} of the output estimates}%
\begin{equation} \label{eq:sim-sys-siso-KF3-Q}
	\begin{aligned}
		\mathbf{Q}_{\text{KF3}}=\mathbf{Q}_{\text{opt}}=\begin{bmatrix}
			\sigma_{x}^2 & 0 \\
			0 & \sigma_{w_p,\text{opt}}^2 \\
		\end{bmatrix}=\begin{bmatrix}
			0.1^2 & 0 \\
			0 & 0.1^2 \\
		\end{bmatrix}.
	\end{aligned}
\end{equation}

Figure \ref{fig:sim-sys-siso-KF123-est} shows the estimates of the three Kalman filters described above when simulated on the input-output data. The upper plot shows the output estimates and the lower plot shows the estimates of the model state $x_{a,2}(k)$, which corresponds to the estimate of the \gls{RODD}, $p(k)$. As expected, the response of KF1 to the infrequent step disturbances is noticeably slower than that of the other two. The response of KF2 is fast, but it is sensitive to the measurement noise. KF3 is a compromise between the two—a relatively fast response to the infrequent steps and less sensitivity to the measurement noise.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim1_all_seed_y_est1_KF123.pdf}
	\caption{Kalman filter estimates}
	\label {fig:sim-sys-siso-KF123-est}
\end{figure}

One way to visualise the observer performance is to plot the square of the output estimation errors, which are the differences between the estimates, $\hat{y}(k \mid k)$, and the true values of the system outputs, $y(k)$. The upper plot in Figure \ref{fig:sim-sys-siso-KF123-cumerr} shows the squared errors in the output estimates of the three Kalman filters. The lower plot shows the cumulative sum of the squared errors. This clearly shows the differences in performance. The estimation errors of KF1 are small when the system is in steady-state but increase dramatically after each step disturbance. The magnitude of the errors of KF2 during steady-state periods is high but roughly constant throughout the simulation. KF3 achieves the lowest overall cumulative sum-of-squared-errors by the end of the simulation.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim1_all_seed_cum_err_y1_KF123.pdf}
	\caption{Cumulative errors of Kalman filter estimates}
	\label{fig:sim-sys-siso-KF123-cumerr}
\end{figure}

The RMSEs of the state estimates in these simulations are shown in Figure \ref{fig:rod-obs-sim1-KF123-xest-RMSE-bar}. These plots show the total RMSEs divided into two components. The blue part, labelled `no noise' is the RMSE when the system was simulated without measurement noise, (i.e. $y_M(k)=y(k)$). The orange part, labelled `due to noise' represents the difference in the RMSE between simulating with measurement noise and without. This distinction is only possible with a simulation model---in real applications it is not possible to obtain the true system outputs.
\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{images/rod_obs_sim1_all_seed_x_err_bar_KF123.pdf}
	\caption{Root-mean-squared errors of state estimates – \acrshort{SISO} system}
	\label{fig:rod-obs-sim1-KF123-xest-RMSE-bar}
\end{figure}

As well as allowing the overall RMSE values of the three Kalman filters to be compared, these plots also illustrate the relative magnitudes of the sources of the error. As expected, the errors of KF1 are almost all due to its slow dynamic response to the shocks, whereas very little is due to the measurement noise. The opposite is true of KF2, which has a very fast dynamic response to the shocks. The trade-off between the two sources of errors is clearly illustrated for KF3. It should be noted that this may not be the best trade-off in every case. In some applications, the speed of response may be more important than the sensitivity to noise, or vice-versa. The trade-off resulting from minimizing the overall RMSE is used in this work.

\subsubsection{Selection of multiple-model observer parameters} \label{sec:sim-obs-lin-1-MKF-tuning}

Parameter values for the sub-optimal multiple-model observer algorithms described in Section \ref{sec:fusion} were chosen by a similar method. The sequence fusion algorithms proposed by \cite{robertson_detection_1995} and \cite{robertson_method_1998} have three parameters, $N_f$, $n_\text{max}$, and $d$. Candidate values for these parameters were generated by considering every combination satisfying

\begin{equation} \label{eq:sim-sys-siso-MKF-SF-param-values}
	\begin{aligned}
		\frac{N_f}{N_d} &\in \left\{3, 5, 10, 15, 20\right\},  \\
		n_\text{max} &\in \left\{1, 2, 3\right\},  \\
		N_d &\in \left\{1, 2, 3, 5, 10\right\}.
	\end{aligned}
\end{equation}

Note that $\frac{N_f}{d}$ represents the number of detection intervals within the fusion horizon and is always a whole number. For example, when $d=3$, the lengths of the fusion were 9, 15, 30, 45, and 60. However, not all combinations were simulated. Any combination with more than 200 hypotheses ($n_h>200$) or with a total probability modelled, $\beta$ \eqref{eq:p_gamma}, less than 0.85 was rejected. The remaining 58 combinations were then evaluated by simulating the observer and calculating the \gls{RMSE} of the output estimates on the 5000 data samples from the simulated system. Tables \ref{tb:obs-sim1-popt-SF95} and \ref{tb:obs-sim1-popt-SF98} show the resulting \gls{RMSE} results for the 10 best combinations of parameter values (those with the lowest output \gls{RMSE}) for the 1995 and 1998 variants of this algorithm.
\begin{table}[ht]
	\begin{center}
		\caption{Observer parameter search results, \acrshort{SISO} linear system, MKF--SF95.} \label{tb:obs-sim1-popt-SF95}
		% See: https://texblog.org/2019/06/03/control-the-width-of-table-columns-tabular-in-latex/
		\begin{tabular}{p{0.05\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.15\textwidth}}
			$N_f$ & $n_\text{max}$ & $d$ & $n_h$ & $\beta$ & $\text{\gls{RMSE}}(\hat{Y},Y)$  \\
			\hline
			% See script rod_obs_sim1_MKF_SF95_popt_table.m
			% 26-Apr-2023 13:31:25 results with seed = 6, sigma_M = 0.1, d_min = 1
			15 &   2 &   1 & 151 & 0.9996 & 0.0411 \\
			10 &   2 &   1 &  76 & 0.9999 & 0.0411 \\
			5 &   1 &   1 &   8 & 0.9990 & 0.0415 \\
			5 &   2 &   1 &  26 & 1.0000 & 0.0415 \\
			5 &   3 &   1 &  48 & 1.0000 & 0.0415 \\
			15 &   1 &   1 &  18 & 0.9904 & 0.0418 \\
			10 &   1 &   1 &  13 & 0.9957 & 0.0419 \\
			20 &   2 &   2 &  56 & 0.9991 & 0.0425 \\
			10 &   2 &   2 &  16 & 0.9999 & 0.0426 \\
			10 &   3 &   2 &  26 & 1.0000 & 0.0426 \\
			\hline
		\end{tabular}
	\end{center}
\end{table}
\begin{table}[ht]
	\begin{center}
		\caption{Observer parameter search results, \acrshort{SISO} linear system, MKF--SF98.} \label{tb:obs-sim1-popt-SF98}
		% See: https://texblog.org/2019/06/03/control-the-width-of-table-columns-tabular-in-latex/
		\begin{tabular}{p{0.05\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.15\textwidth}}
			$N_f$ & $n_\text{max}$ & $d$ & $n_h$ & $\beta$ & $\text{\gls{RMSE}}(\hat{Y},Y)$  \\
			\hline
			% See script rod_obs_sim1_MKF_SF98_popt_table.m
			% 26-Apr-2023 12:30:13 results with seed = 6, sigma_M = 0.1, d_min = 1
			30 &   2 &   2 & 151 & 0.9970 & 0.0414 \\
			20 &   2 &   2 &  76 & 0.9991 & 0.0414 \\
			10 &   2 &   2 &  26 & 0.9999 & 0.0414 \\
			10 &   3 &   2 &  48 & 1.0000 & 0.0414 \\
			6 &   1 &   2 &   6 & 0.9988 & 0.0414 \\
			6 &   2 &   2 &  13 & 1.0000 & 0.0414 \\
			6 &   3 &   2 &  16 & 1.0000 & 0.0414 \\
			10 &   1 &   2 &   8 & 0.9962 & 0.0416 \\
			9 &   1 &   3 &   6 & 0.9974 & 0.0418 \\
			9 &   2 &   3 &  13 & 1.0000 & 0.0418 \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

From a practical standpoint, it is not desirable to simulate a large number of hypotheses. On this basis, the third combination in Table \ref{tb:obs-sim1-popt-SF95} ($N_f=10,n_\text{max}=2,d=1$, requiring 76 hypotheses) or the fifth ($N_f=5,n_\text{max}=1,d=1$, requiring only 8 hypotheses) are good choices for the 1995 version. The fifth combination was chosen for the simulations and the observer was labelled `MKF--SF95'. 

In terms of overall estimation errors, the two sequence fusion algorithms are very similar. Note that in the case of the 1998 version, it produces identical results to the MKF--SF95 observer when the spacing parameter $d$ is 1. The differences between the two algorithms are only in the way the estimates are updated during the steps of the detection interval. Therefore, results for the 1998 algorithm with $d=1$ were discarded. \cite{robertson_method_1998} state that the approach is a means of obtaining a long detection horizon, which is necessary when many sampling intervals are required to discriminate among multiple abrupt changes that may occur. Since in this simulation there is only one \gls{RODD} disturbance, and the best results were achieved with a detection interval of 1, it was concluded that these features of the 1998 algorithm were not applicable to this problem.

Nevertheless, based on the results in Table \ref{tb:obs-sim1-popt-SF98}, an observer labelled `MKF--SF1' was selected for simulation with $N_f=6$, $n_\text{max}=1$, and $d=2$. This requires only 6 hypotheses and achieves as good a performance as `MKF--SF95'.

The sequence pruning algorithm proposed by \cite{eriksson_classification_1996} has two parameters that must be specified, $n_h$, and $N_\text{min}$. Candidate values for these parameters were generated by considering all combinations satisfying

\begin{equation} \label{eq:sim-sys-siso-MKF-SP-param-values}
	\begin{aligned}
		n_h &\in \left\{3, 4, 5, 7, 10, 14, 19, 25, 32\right\},  \\
		N_\text{min} &\in \left\{1, 2, 3, 4, 5, 6, 7, 9, 12, 16, 21\right\}.
	\end{aligned}
\end{equation}

Combinations where rejected when $n_h - N_\text{min} < 1$, or in other words, when the minimum life, $N_\text{min}$, was greater than or equal to the number of hypotheses, $n_h$. This is a necessary condition since there must be enough filters to simulate the hypotheses until they reach the minimum life. Each remaining combination was then evaluated in the same way by calculating the \gls{RMSE} of the output estimates on the same 5000 samples used in the tuning of the sequence fusion observers described above. Table \ref{tb:obs-sim1-popt-SP} shows the top 10 combinations of parameter values in terms of lowest \gls{RMSE}.

\begin{table}[ht]
	\begin{center}
		\caption{Observer parameter search results, \acrshort{SISO} linear system, MKF--SP.} \label{tb:obs-sim1-popt-SP}
		% See: https://texblog.org/2019/06/03/control-the-width-of-table-columns-tabular-in-latex/
		\begin{tabular}{p{0.05\textwidth}>{\centering\arraybackslash}p{0.07\textwidth}>{\centering\arraybackslash}p{0.15\textwidth}}
			$n_h$ & $N_\text{min}$ & $\text{\gls{RMSE}}(\hat{Y},Y)$  \\
			\hline
			% See script rod_obs_sim1_MKF_SP_popt_table.m
			% 26-Apr-2023 14:07:02 results with seed = 6, sigma_M = 0.1
			10 &   7 & 0.0409  \\
			10 &   6 & 0.0410  \\
			19 &  16 & 0.0411  \\
			14 &  12 & 0.0411  \\
			25 &  21 & 0.0411  \\
			14 &   7 & 0.0411  \\
			19 &   7 & 0.0412  \\
			19 &   6 & 0.0412  \\
			25 &  16 & 0.0412  \\
			25 &  12 & 0.0412  \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

As above, the differences between the \gls{RMSE}s of the best 10 combinations are small, indicating that the choice of parameters from these combinations is not critical. The first combination ($n_h=10,N_\text{min}=7$) was chosen and the observer was labelled `MKF--SP1' in the simulations below. Table \ref{tb:obs-params-sim1} is a summary of the parameters of the observers simulated.
\begin{table}[ht]
	\begin{center}
		\caption{Observer parameters for \acrshort{SISO} linear system.} \label{tb:obs-params-sim1}
		% See: https://texblog.org/2019/06/03/control-the-width-of-table-columns-tabular-in-latex/
		\begin{tabular}{p{0.16\textwidth}>{\centering\arraybackslash}p{0.12\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}}
			& KF3 & MKF--SF95 & MKF--SF1 & MKF--SP \\
			\hline
			Type & Kalman filter & Multi-model Kalman filter & Multi-model Kalman filter & Multi-model Kalman filter \\
			Algorithm & - & Sequence fusion 95 & Sequence fusion 98 & Sequence pruning \\
			\hline
			Parameters &  &  &  &  \\
			$\mathcal{Q}$ & $\mathbf{Q}_{opt}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ \\
			$\mathbf{R}$ & $0.1^2$ & $0.1^2$ & $0.1^2$ & $0.1^2$ \\
			$\mathbf{P}(0)$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ \\
			$n_h$ & 1 & 8 & 6 & 10 \\
			$N_f$ & - & 5 & 6 & - \\
			$n_\text{max}$ & - & 1 & 1 & - \\
			$d$ & - & 1 & 2 & - \\
			$N_\text{min}$ & - & - & - & 7 \\
			$\varepsilon$ & - & 0.01 & 0.01 & 0.01 \\
			$\sigma_{w_p}$ & - & 0.01 & 0.01 & 0.01 \\
			$b$ & - & 100 & 100 & 100 \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

\subsubsection{Scheduled Kalman filter} \label{sec:sim-obs-lin-1-SKF}

To determine the best performance that could be achieved by an ideal multiple-model observer, a hypothetical observer that has access to the true shock indicator signal, $gamma_k$, was simulated. This `scheduled' Kalman filter, labelled ‘SKF’, has only one shock sequence---the correct one---and therefore the process noise covariance is automatically switched to match the known shock signal. The \gls{SKF} is hypothetical because it could not be implemented in a real application where the random shocks are not known.

\subsubsection{Simulation results} \label{sec:sim-obs-lin-1-results}

Figure \ref{fig:rod-obs-sim1-yest-1-SF} shows the estimates of the sequence fusion observers, MKF--SF95 and MKF--SF1, compared to the scheduled Kalman Filter, \gls{SKF}, over the first 600 time steps of the simulation. From these results it is apparent that the estimates of the sequence fusion observers are generally very close to those of the \gls{SKF}. The estimates respond quickly when the random shocks occur and their variance at other times is relatively small despite the measurement noise. A few spontaneous errors of short duration are also noticeable in both the \gls{SKF} and MKF--SF1 estimates at certain times. However, when the variances are compared to those of the Kalman filters in Figure \ref{fig:sim-sys-siso-KF123-est}, it is clear that these multiple-model observers have an advantage.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim1_all_seed_y_est1_SF95_SF1.pdf}
	\caption{Estimates by sequence fusion observers – \acrshort{SISO} system}
	\label{fig:rod-obs-sim1-yest-1-SF}
\end{figure}

Figure \ref{fig:rod-obs-sim1-yest-1-SP} shows the estimates of the sequence pruning observer, MKF--SP1, compared to the scheduled Kalman Filter. It also responds quickly to the shocks and is insensitive to the measurement noise, similar to the sequence fusion observers.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim1_all_seed_y_est1_SP1.pdf}
	\caption{Estimates by sequence pruning observer – \acrshort{SISO} system}
	\label{fig:rod-obs-sim1-yest-1-SP}
\end{figure}

The plots in Figure \ref{fig:rod-obs-sim1-cum-err-y1-all} show the squared output errors and the cumulative squared-errors of the five observers over the full duration of the simulation. The differences between the multiple-model observers are almost indistinguishable. By the end of the simulation, the cumulative errors are significantly lower than those of KF3 and slightly higher than the \gls{SKF}.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim1_all_seed_cum_err_y1.pdf}
	\caption{Cumulative errors of observer estimates – \acrshort{SISO} system}
	\label{fig:rod-obs-sim1-cum-err-y1-all}
\end{figure}

Figure \ref{fig:rod-obs-sim1-xest-RMSE-bar} compares the \gls{RMSE}s of the state estimates of the five observers over the full simulation. From this, it can be seen that the multiple-model observers achieve lower \gls{RMSE}s than KF3, both in simulations with and without noise. The reduction in the errors due to the lower sensitivity to measurement noise (orange bars) is a more significant contribution to the improved performance of the multiple-model observers. It is noteworthy that most of the errors in the estimates of state 2, which corresponds to the input disturbance, are not due to measurement noise. This may partly be attributed to the delay in the process, which makes it impossible for the observers to immediately detect an unmeasured disturbance. They are also due to errors in tracking the persistent component of the disturbance, which can be seen in Figures \ref{fig:rod-obs-sim1-yest-1-SF} and \ref{fig:rod-obs-sim1-yest-1-SP}.
\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{images/rod_obs_sim1_all_seed_x_err_bar.pdf}
	\caption{Root-mean-squared errors of state estimates – \acrshort{SISO} system}
	\label{fig:rod-obs-sim1-xest-RMSE-bar}
\end{figure}

Because these processes are stochastic, it is important to run simulations of sufficient duration to be confident that the results are representative of average performance. Figure \ref{fig:rod-obs-sim1-yest-all-seed-RMSE-box} provides an indication of the sensitivity of the \gls{RMSE} results to the initialization of the pseudo-random number generator used in these simulations. The data for this figure are from the results of 10 simulations, each of 5000 samples and each with a different pseudo-random initialization. It can be seen that the \gls{RMSE} values are somewhat sensitive to the random initialization but the variation is an order of magnitude lower than the difference between the multiple-model observers and KF3.
\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{images/rod_obs_sim1_all_seed_y_err_box.pdf}
	\caption{Sensitivity of output estimation errors to random initialization – \acrshort{SISO} system}
	\label{fig:rod-obs-sim1-yest-all-seed-RMSE-box}
\end{figure}
Figure \ref{fig:rod-obs-sim1-yest-all-seed-RMSE-box} confirms that the multiple-model observers are a significant improvement on the best single Kalman filter. However, there is no significant difference between the \gls{MKF}s in this case.
% Plots for analysis of SF performance
%
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=15cm]{images/rod-obs-sim-1-4-wfplot-DRAFT.png}
%	\caption{Evolution of conditional probabilities with sequence fusion}
%	\label{fig:rod-obs-sim1-4-wfplot}
%\end{figure}
%
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=14cm]{images/rod-obs-sim-1-4-est-MKF-SF-plot-DRAFT.pdf}
%	\caption{Filter estimates of sequence fusion observer during step disturbances}
%	\label{fig:rod-obs-sim1-4-est-MKF-SF-plot-DRAFT}
%\end{figure}
%
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=14cm]{images/rod_obs_sim2_7_mse_ashocks.pdf}
%	\caption{Average errors of observer estimates after shock disturbances}
%	\label{fig:rod_obs_sim2_7_mse_ashocks}
%\end{figure}

\subsection{MIMO linear system} \label{sec:sim-obs-lin-2}

\cite{robertson_method_1998} stated that in the case of more than one \gls{RODD} disturbance, it often takes many sampling intervals to discriminate among possible disturbances. The main motivation for the detection intervals used in their approach was to achieve a longer \textit{detection horizon} without requiring a very large number of hypotheses. To evaluate this feature, a \gls{MIMO} system with a 2-input, 2-output ($2\times2$) process model was chosen. The system has two independent \gls{RODD} step disturbances added to each manipulated input, as shown in the functional diagram in Figure \ref{fig:sim-sys-diag-2x2}. 
\begin{figure}[htp]
	\centering
	\includegraphics[width=11.5cm]{images/sim-sys-diag-2x2.pdf}
	\caption{Functional diagram of the simulated \acrshort{MIMO} system}
	\label{fig:sim-sys-diag-2x2}
\end{figure}

The process model is a discrete-time linear model derived from the $2\times2$ continuous-time system,
% 2x2 system #1 - see sys_rodin_step_2x2sym.m
%\begin{equation} \label{eq:sim-sys-mimo-ct}
%	\mathbf{G}(s) = \left[\begin{array}{cc}
%		G_{11}(s) & G_{12}(s)  \\
%		G_{21}(s) & G_{22}(s)  \\
%	\end{array}\right] = \left[\begin{array}{cc}
%		\frac{1}{1+8.5s} & \frac{-0.5}{1+8.5s}  \\
%		\frac{-0.5}{1+8.5s} & \frac{1}{1+8.5s}  \\
%	\end{array}\right]
%\end{equation}
% 2x2 system #2 - see sys_rodin_step_2x2sym2.m
\begin{equation} \label{eq:sim-sys-mimo-ct}
	\mathbf{G}(s) = \left[\begin{array}{cc}
		G_{11}(s) & G_{12}(s)  \\
		G_{21}(s) & G_{22}(s)  \\
	\end{array}\right] = \left[\begin{array}{cc}
		\frac{1}{1+8.5s} & \frac{0.5}{1+8.5s}  \\
		\frac{0.5}{1+8.5s} & \frac{1}{1+8.5s}  \\
	\end{array}\right],
\end{equation}
where $s$ is the Laplace variable.  This model was chosen because the dynamics of each input-output combination are similar (they have the same time constant and the sign of the gain is the same). The only differences are in the magnitude of the static gains. The model was converted to a discrete-time model with a sampling period of one ($T_s=1$). The discrete-time state space model of the augmented system used by the observers was
\begin{equation} \label{eq:sim-sys-2x2-ss-aug}
	\begin{split}
		\mathbf{x}_{a}(k+1) & =\left[\begin{array}{cccc}
			0.8890 & 0 & 1 & 0.5 \\
			0 & 0.8890 & 0.5 & 1 \\
			0 & 0 & 1 & 0 \\
			0 & 0 & 0 & 1
		\end{array}\right] \mathbf{x}_{a}(k) + \left[\begin{array}{cc}
			1 & 0.5 \\
			0.5 & 1
		\end{array}\right] \mathbf{u}(k) + \mathbf{w}_{a}(k) \\
		\mathbf{y}_M(k) & =\left[\begin{array}{cccc}
			0.1110 & 0 & 0 & 0 \\
			0 & 0.1110 & 0 & 0
		\end{array}\right] \mathbf{x}_{a}(k) + \mathbf{v}(k)
	\end{split}
\end{equation}
where
\begin{equation} \label{eq:sim-sys-2x2-ss-aug2}
	\mathbf{x}_{a}(k) = \left[\begin{array}{l}
		x_{1}(k) \\
		x_{2}(k) \\
		p_{1}(k) \\
		p_{2}(k)
	\end{array}\right], \mathbf{w}_{a}(k) = \left[\begin{array}{l}
		w_1(k) \\
		w_2(k) \\
		w_{p,1}(k) \\
		w_{p,2}(k)
	\end{array}\right].
\end{equation}

Note that with this representation, the third and fourth model states, $x_{a,3}(k)$ and $x_{a,4}(k)$, correspond exactly to the two input disturbances, $p_1(k)$ and $p_2(k)$.

The two random shocks were simulated with the same parameters as those in the \gls{SISO} simulation, $\varepsilon_i=0.01, \sigma_{w_{p,i}}=0.01$, and $b_i=100$. However, the standard deviations of the measurement noises, $\sigma_{M,1}$ and $\sigma_{M,2}$, were 0.2, which is twice as high as they were in the \gls{SISO} simulations.

Figure \ref{fig:rod-obs-sim-2-ioplot} shows the first 300 input-output samples from a simulation of the system with a total duration of 5000. The lower of the two plots shows the input signals, $u_1(k)$, $u_2(k)$, $p_1(k)$, and $p_2(k)$. In this simulation, there were no changes to the known input variables. The upper plot shows the simulated output measurements, $y_1(k)$ and $y_2(k)$. No input disturbances other than $w_{p,1}(k)$ and $w_{p,2}(k)$ were simulated (i.e. $w_1(k)=w_2(k)=0$).
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim2_all_seed_ioplot.pdf}
	\caption{Simulation of linear \acrshort{MIMO} system with two \acrshort{RODD}s}
	\label{fig:rod-obs-sim-2-ioplot}
\end{figure}

As in the previous experiment, a standard Kalman filter was tuned experimentally to achieve a low \gls{RMSE} of the output estimates. The details of this tuning process are included in Annex \ref{sec:annex-sim-2-KF-tuning} including the resulting value of $\mathbf{Q}_{\text{KF3}}$ \eqref{eq:sim-sys-sim2-KF3-Q}.


\subsubsection{Selection of multiple-model observer parameters} \label{sec:sim-obs-lin-2-MKF-tuning}

Parameter values for the sub-optimal multiple-model observers were chosen by a similar method to that used in the SISO experiment. The details of the parameter search and selection process are in Annex \ref{sec:annex-sim-2-MKF-tuning}. Table \ref{tb:obs-params-sim2} summarises the observer parameters used in the simulations.
\begin{table}[ht]
	\begin{center}
		\caption{Observer parameters for $2\times2$ linear system.} \label{tb:obs-params-sim2}
		% See: https://texblog.org/2019/06/03/control-the-width-of-table-columns-tabular-in-latex/
		\begin{tabular}{p{0.16\textwidth}>{\centering\arraybackslash}p{0.12\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}}
			& KF3 & MKF--SF95 & MKF--SF1 & MKF--SP1 \\
			\hline
			Type & Kalman filter & Multi-model Kalman filter & Multi-model Kalman filter & Multi-model Kalman filter \\
			Algorithm & - & Sequence fusion 95 & Sequence fusion 98 & Sequence pruning \\
			\hline
			Parameters &  &  & &  \\
			$\mathcal{Q}$ & $\mathbf{Q}_{opt}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ \\
			$\mathbf{R}$ & $\left[\begin{smallmatrix}0.2^2 & 0 \\ 0 & 0.2^2\end{smallmatrix}\right]$
			& $\left[\begin{smallmatrix}0.2^2 & 0 \\ 0 & 0.2^2\end{smallmatrix}\right]$
			& $\left[\begin{smallmatrix}0.2^2 & 0 \\ 0 & 0.2^2\end{smallmatrix}\right]$
			& $\left[\begin{smallmatrix}0.2^2 & 0 \\ 0 & 0.2^2\end{smallmatrix}\right]$ \\
			$\mathbf{P}(0)$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ \\
			$n_h$ & 1 & 116 & 58 & 61 \\
			$N_f$ & - & 15 & 15 & - \\
			$n_\text{max}$ & - & 2 & 2 & - \\
			$d$ & - & 3 & 5 & - \\
			$N_\text{min}$ & - & - & - & 29 \\
			$\varepsilon$ & - & $\begin{bsmallmatrix}0.01 & 0.01\end{bsmallmatrix}^\intercal$
			& $\begin{bsmallmatrix}0.01 & 0.01\end{bsmallmatrix}^\intercal$ 
			& $\begin{bsmallmatrix}0.01 & 0.01\end{bsmallmatrix}^\intercal$ \\
			$\sigma_{w_p}$ & - & $\begin{bsmallmatrix}0.01 & 0.01\end{bsmallmatrix}^\intercal$
			& $\begin{bsmallmatrix}0.01 & 0.01\end{bsmallmatrix}^\intercal$ 
			& $\begin{bsmallmatrix}0.01 & 0.01\end{bsmallmatrix}^\intercal$ \\
			$b$ & - & $\begin{bsmallmatrix}100 & 100\end{bsmallmatrix}^\intercal$
			& $\begin{bsmallmatrix}100 & 100\end{bsmallmatrix}^\intercal$ 
			& $\begin{bsmallmatrix}100 & 100\end{bsmallmatrix}^\intercal$ \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

\subsubsection{Simulation results} \label{sec:sim-obs-lin-2-results}

Figure \ref{fig:rod-obs-sim2-yest-1-SF} shows the estimates of both sequence fusion observers (MKF--SF95 and MKF--SF1) compared to the scheduled Kalman Filter (\acrshort{SKF}) over the first 300 time steps of the simulation. The upper two plots show the estimates of outputs 1 and 2 compared to the true system outputs. The lower two plots show the estimates of states 3 and 4, which correspond to the input disturbances. The largest errors occur in the periods following the shocks. At these times, the estimates of MKF--SF95 and MKF--SF1 lag those of the \gls{SKF}. At other times, the estimates of all three observers are quite similar. There is no noticeable difference between MKF--SF95 and MKF--SF1.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim2_all_seed_y_est1_SF95_SF1.pdf}
	\caption{Estimates by sequence fusion observer –  $2\times2$ system}
	\label{fig:rod-obs-sim2-yest-1-SF}
\end{figure}

Figure \ref{fig:rod-obs-sim2-yest-1-SP} shows the estimates of the sequence pruning observer (MKF--SP) compared to the scheduled Kalman Filter on the same simulation data. The behaviour is similar to that of the sequence fusion observers, matching the \gls{SKF} during periods between shocks and lagging it in the periods immediately after.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim2_all_seed_y_est1_SP1.pdf}
	\caption{Estimates by sequence pruning observer –  $2\times2$ system}
	\label{fig:rod-obs-sim2-yest-1-SP}
\end{figure}

%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=13cm]{images/rod_obs_sim2_all_seed_y_est2_MKF.pdf}
%	\caption{Comparison of multiple-model observer estimates –  $2\times2$ system}
%	\label{fig:rod-obs-sim2-yest-2-MKF}
%\end{figure}
Figure \ref{fig:sim-sys-sim2-MKF-cumerr} shows the squared output estimation errors and the cumulative squared-errors of the five observers over the full simulation. This reveals that on average, the sequence fusion observers (MKF--SF95 and MKF--SF1) perform better than the sequence pruning observer. The cumulative errors of all the multiple-model observers are significantly lower than those of KF3 by the end of the simulation, but not as close to those of the \gls{SKF} as they were in the SISO system simulation. This may be attributed to the delays in detecting the shocks in this simulation.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/rod_obs_sim2_all_seed_cum_err_y2.pdf}
	\caption{Cumulative sum of squared errors of output estimates –  $2\times2$ system}
	\label{fig:sim-sys-sim2-MKF-cumerr}
\end{figure}

% Replaced below plots with above single plot to save space
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=13cm]{images/rod_obs_sim2_cum_err2.pdf}
%	\caption{Cumulative errors of multiple-model observer output estimates –  $2\t_y2mes2$ system}
%	\label{fig:sim-sys-sim2-MKF-cumerr2}
%\end{figure}
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=11cm]{images/rod_obs_sim2_all_seed_x_err_bar.pdf}
%	\caption{Root-mean-squared errors of state estimates – $2\times2$ system}
%	\label{fig:rod-obs-sim2-xest-RMSE-bar}
%\end{figure}
As in the \gls{SISO} case, 10 simulations were carried out with different random initializations to check the sensitivity of the results. Figure \ref{fig:rod-obs-sim2-yest-all-seed-RMSE-box} shows the variation in the \gls{RMSE}s. This confirms that the performance of the three multiple-model observers is quite similar but MKF--SF95 and MKF--SF1 achieve a slightly lower error on average.
\begin{figure}[htp]
	\centering
	\includegraphics[width=12cm]{images/rod_obs_sim2_all_seed_y_err_box.pdf}
	\caption{Root-mean-squared errors of output estimates – $2\times2$ system}
	\label{fig:rod-obs-sim2-yest-all-seed-RMSE-box}
\end{figure}

\section{Grinding process simulations} \label{sec:sim-ore-SISO}

In this section, results of an experiment to test the observers with data from the grinding process simulator are presented. Unlike the simulations presented in section \ref{sec:sim-obs-lin}, the grinding simulator is a non-linear model with a high number of state variables and complex dynamics. As described in Section {\ref{sec:grinding-simulator}}, it is a phenomenological model of the grinding process with ore streams modelled as populations of ore particles in discrete size intervals, a recycle stream with a transport delay, and simulated {\gls{PI}} controllers to regulate key process variables. When the {\gls{PSD}} of the ore feed switches from one mixture to the other, the process responds dynamically and will eventually reach a new steady-state, if the ore does not change again. However, when successive ore changes occur within short periods, the internal state of the simulated process at the time of the change has an influence on the subsequent response. Thus, the dynamic response of the measured output variable is non-linear and difficult to predict accurately.

This work is an extension of previously published work \citep{tubbs_observer_2022}. In the previous work, the same data from the grinding simulator was used to test a multiple-model observer with sequence pruning and similar results were presented. However, the work here has some differences. Here, two types of observer with sequence fusion and sequence pruning were evaluated and compared. Also, note that in the previous work, the prediction form of the Kalman filter \eqref{eq:xkp1_hat_p} was used and the mean-squared errors (MSEs) of the estimates were compared, rather than the \gls{RMSE}s \eqref{eq:rmse-calc-yest}.

As in the previous work, the particle size (\acrshort{P80}) of the ground product leaving the cyclone overflow was sampled at 3 minute intervals ($T_s=0.05$ hours) and a zero-mean Gaussian measurement noise with a standard deviation of 5 $\mu\text{m}$ was added directly to the simulator output. No disturbances, other than the switching ore feed were simulated and all manipulatable inputs other than those manipulated by the {\gls{PI}}s, such as ore feed rate, rotation speed, ... etc., were held constant during the simulations. In total, data from 13 different runs of the simulation model were used. Table \ref{tb:grind1-sims} summarises the type of input signal used in each simulation, the duration, number of samples, and the simulation outputs that were used. The use of this data is described in the subsequent sections.
\begin{table}[ht]
	\begin{center}
		\caption{Grinding process simulations} \label{tb:grind1-sims}
		% See: https://texblog.org/2019/06/03/control-the-width-of-table-columns-tabular-in-latex/
		\begin{tabular}{
				>{\raggedleft\arraybackslash}p{0.32in}
				>{\centering\arraybackslash}p{0.8in}
				>{\centering\arraybackslash}p{0.6in}
				>{\centering\arraybackslash}p{0.6in}
				>{\centering\arraybackslash}p{0.8in}
				>{\raggedright\arraybackslash}p{1.8in}}
			\# & Input type & Duration (hours) & Samples ($N$) & Output data & Purpose \\
			\hline
			1 & \gls{PRBS} & 15 & 300 & $p(k),y_M(k)$ & Estimation of process model parameters \\
			2 & \gls{PRBS} & 15 & 300 & $p(k),y_M(k)$ & Validation of process models and model selection  \\
			3 & \gls{PRBS} & 15 & 300 & $p(k)$, $y(k)$, $y_M(k)$ & Figure \ref{fig:grind1_rod_obs_est_MKF_SF1}  \\
			6 & Random switching & 123 & 2460 & $y(k)$, $y_M(k)$ & Observer parameter selection \\
			6--15 & Random switching & 123 & 2460 & $y(k)$, $y_M(k)$ & Observer simulation and evaluation. \\
			\hline
		\end{tabular}
	\end{center}
\end{table}


\subsection{System identification} \label{sec:grind1-sysid}

The system model consists of a single-input, single-output (SISO) process model with a \gls{RODD} step disturbance, $p(k)$, at the input, representing step-changes in the ore mix factor, and a measured output of the P80 of the ground product. The functional diagram of this system is shown in Figure \ref{fig:grind1_obs_model}. Note that there is no known input, $u(k)$, available to the observer.
\begin{figure}[ht]
	\centering
	\includegraphics[width=9.5cm]{images/grind1-obs-model-diag.pdf}
	\caption{Observer model structure}
	\label{fig:grind1_obs_model}
\end{figure}

To construct multiple-model observers for this system, the structure and parameters of the disturbance model and process model are needed. As mentioned in Section \ref{sec:sys-id}, standard system identification techniques are not suitable for estimating hybrid dynamical systems. To avoid this problem, it is assumed here that the time and magnitude of the disturbances (i.e. changes in the ore mix) have been determined for a period of time sufficient for identification of a process model. In practice, this might be possible by sampling the ore and carrying out a particle size analysis, or in the case where stocks of ores with known properties can be fed to the mill in a controlled manner for the duration of the data acquisition.

To mimic the system identification process as realistically as possible, three simulations were carried out (numbered 1--3 in Table \ref{tb:grind1-sims}) with a \textit{pseudo-random binary sequence} (\glsadd{PRBS}\gls{PRBS}) of step-changes in the ore mix factor and the process output measurements were recorded. Each simulation was started with the same initial state but run with a different disturbance input sequence, $p(0), p(1), ..., p(N)$. The simulations were of 15 hours duration ($N=300$) and the process was kept at steady-state for the first 3 hours of each simulation. The initial 3-hour period of the first simulation was used to estimate the normal operating point of the process. Each step in the \gls{PRBS} was one hour in duration. The lower plot in Figure \ref{fig:grind1-sim_ioplots} shows the values of $p(k)$ used in simulation 1. The simulated measurements, $y_M(k)$, are shown in the upper plot.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/grind1_rod_obs_sim_1_ioplot_P2DcTd4.pdf}
	\caption{Simulation inputs and outputs}
	\label{fig:grind1-sim_ioplots}
\end{figure}

The Mathworks System Identification Toolbox™ was used to carry out the system identification. The data from simulation 1 was used to estimate the parameters of candidate process models. A model structure was selected after fitting a number of different models and calculating the prediction errors on the input-output data from simulation 2. Both discrete and continuous time models were fitted and compared. To ensure that the identified models were realistic, it was necessary to constrain some of the parameters, and it was found that the continuous time models were easier to work with in this respect. After checking the significance of the estimated model parameters, a simple second-order model with a delay was selected. The transfer function of this model is
% Model name P2DcTd4
\begin{equation} \label{eq:grind1-id-model-ctf}
	Y(s)= \frac{-32.4\exp(-0.2s)}{(1 + 0.106s)^2}U(s),
\end{equation}
\nomenclature{$s$}{Laplace variable}%
\nomenclature{$Y(s)$}{Laplace transform of the process output signal}%
\nomenclature{$U(s)$}{Laplace transform of the process input signal}%
where $Y(s)$ and $U(s)$ here, are the Laplace transforms of the continuous-time output and input signals.

This model has a delay of 0.2 hours (i.e. 4 sample periods, or 12 minutes) and a settling time (to within $\pm5\%$ of steady-state) of approximately 42 minutes. The upper plot in Figure \ref{fig:grind1-sim_ioplots} also shows the predictions of this model, $y_{model}(k)$, compared to the true output, $y(k)$. Note that there is a visible mismatch between the predictions of the identified linear model and the outputs of the grinding process simulator. This may be attributed partly to the measurement errors, and partly to the non-linear nature of the true response.

It was also assumed that the structure of the \gls{RODD} model \eqref{eq:RODD} was known, in this case a \gls{RODD} step disturbance \eqref{eq:RODD-step}, and that values of the parameters $\sigma_{w_p}$, $b$, and $\varepsilon$ were available. In practice, it should be possible to estimate $\varepsilon$ from historic operating data by simply counting the number of ore changes in a sufficiently long period and dividing by the number of sample periods, although this method is somewhat inaccurate because small shocks are undetectable in noisy data. Estimating the magnitude of the shocks, $b\sigma_{w_p}$, could be more difficult. It may require prior knowledge of the true size distributions of the ores, perhaps from geological data or test work.

\subsection{Observer design} \label{sec:grind1-obs-design}

The discrete-time state-space representation of the augmented model (including \gls{RODD} disturbance) used by the observers was
% P2DcTd4
\begin{multline} \label{eq:grind1-obs_ss_model}
	\mathbf{x}(k+1) = 
	\begin{bmatrix}
		1.248 & -0.7786 & 4 & 0 & 0 & 0 & 0 \\
		0.5   &  0      & 0 & 0 & 0 & 0 & 0 \\
		0     &  0      & 0 & 1 & 0 & 0 & 0 \\
		0     &  0      & 0 & 0 & 1 & 0 & 0 \\
		0     &  0      & 0 & 0 & 0 & 1 & 0 \\
		0     &  0      & 0 & 0 & 0 & 0 & 1 \\
		0     &  0      & 0 & 0 & 0 & 0 & 1 \\
	\end{bmatrix} \mathbf{x}(k)
	+ \begin{bmatrix}
		\mathbf{w}(k) \\
		w_p(k) \\
	\end{bmatrix}, \\
	y_M(k) = 
	\begin{bmatrix}
		-0.66214 & -0.96672 & 0 & 0 & 0 & 0 & 0 \\
	\end{bmatrix} \mathbf{x}(k) + \mathbf{v}(k) \\
\end{multline}
with a sampling period of 0.05 hours.

The noise covariance \eqref{eq:init_Q_R} switched between two values,
\begin{equation} \label{eq:Q0}
	\begin{aligned}
		\mathbf{Q}_0=\begin{bmatrix}
			\sigma_{x_1}^2 & 0 & \cdots & 0 \\
			0 & \sigma_{x_2}^2 & \cdots & 0 \\
			\vdots & \vdots & \ddots & \vdots \\
			0 & 0 & \cdots & \sigma_{w_p}^2 \\
		\end{bmatrix}=\begin{bmatrix}
			0.01^2 & 0 & \cdots & 0 \\
			0 & 0.01^2 & \cdots & 0 \\
			\vdots & \vdots & \ddots & \vdots \\
			0 & 0 & \cdots & 0.002717^2
		\end{bmatrix}, \\
	\end{aligned}
\end{equation}
and 
\begin{equation} \label{eq:Q1}
	\begin{aligned}
		\mathbf{Q}_1=\begin{bmatrix}
			\sigma_{x_1}^2 & 0 & \cdots & 0 \\
			0 & \sigma_{x_2}^2 & \cdots & 0 \\
			\vdots & \vdots & \ddots & \vdots \\
			0 & 0 & \cdots & b^2\sigma_{w_p}^2 \\
		\end{bmatrix}=\begin{bmatrix}
			0.01^2 & 0 & \cdots & 0 \\
			0 & 0.01^2 & \cdots & 0 \\
			\vdots & \vdots & \ddots & \vdots \\
			0 & 0 & \cdots & 0.2717^2
		\end{bmatrix},
	\end{aligned}
\end{equation}
where $\sigma_{x_i}^2$ ($i=1,2,...,6$) are noise variance parameters of the process model states. Since model errors are not known in practice, an arbitrary value of $0.01^2$ was chosen, although it did not have a significant impact on the results. In this work, the random shock variances, $\sigma_{w_p}^2$ and $b^2\sigma_{w_p}^2$, were set according to the simulated ore changes, as described below.

To evaluate the observers, longer simulations with randomly-occurring step changes in the ore mix factor were simulated. Each simulation consisted of an initial period of 3 hours at steady-state followed by 120 hours (equivalent to 5 days) of randomly-occurring step changes. Note that this is not a \gls{RODD} according to the definition because the magnitude of the shocks is not random. \gls{RODD}s are unbounded and therefore not realistic representations of actual process disturbances, which are unlikely to be Gaussian random variables. Instead, a simple switching ore mix factor, $p_\text{ore}(k)$, was used. $p_\text{ore}(k)$ is a discrete random variable with two possible values, that switches according to the shock probability $\varepsilon$. It can be represented by a \gls{HMM} with one binary state, $r(k)$, such that
\nomenclature{$p_\text{ore}(k)$}{simulated ore disturbance signal at time $k$}%
\begin{equation} \label{eq:pk-grind}
	p_\text{ore}(k) = \begin{cases*}
		0.2283 & if $r(k)=1$, \\
		0.5 & if $r(k)=2$,
	\end{cases*} \\
\end{equation}
and with the transition probability matrix
\begin{equation} \label{eq:Pi-rk-grind}
	\Pi_{p_\text{ore}} = \begin{bmatrix}
		1-\varepsilon & \varepsilon \\
		\varepsilon & 1-\varepsilon
	\end{bmatrix}.
\end{equation}

The parameters of the \gls{RODD} model used in the observer design were chosen in the following way. Although the probability density of $p_\text{ore}(k)$ is obviously not Gaussian, it has a conditional standard deviation of either 0.2717 (the difference between 0.5 and 0.2283) or 0. Based on this, $b\sigma_{w_p}$ was set to 0.2717 and $b$ was arbitrarily set to 100, meaning that $\sigma_{w_p}$ was 0.002717. $\varepsilon$ was set to 0.01, which was the true value used in \eqref{eq:Pi-rk-grind} to simulate $p_\text{ore}(k)$. This means that, on average, an ore change occurred every 100 sample periods, or every 5 hours.

Three multiple-model observers were designed using a similar method to that described in Section \ref{sec:sim-obs-lin-1-MKF-tuning}. It was found that large data sets from long simulations were necessary to obtain good combinations of parameter values for these observers. When this was attempted with shorter data sets like simulations 1--3, the results were variable, and the performance of the observers was often poor. This can be explained by the fact that shocks are infrequent, and also that the sequence fusion algorithms in particular, are somewhat sensitive to the exact timing of shocks. The simulations must be long enough to provide a sufficient number and diversity of shock sequences to avoid `overfitting' the parameters to the data.

For this reason, data from one of the longer simulations was used (number 6) to obtain the parameter values of the observers evaluated below. This simulation was 5 days in duration (2460 samples), simulated with $p_\text{ore}(k)$ as the input, and included 32 step changes. This was sufficient to obtain good parameter values. Note that the true values of the disturbance are not needed for observer parameter selection.

However, in practice it would not be possible to use the true process outputs because they are unknown. In that case, the parameter tuning would have to be done using the \gls{RMSE} of the measured outputs, $\text{RMSE}(\hat{Y}(N), Y_M(N))$. The true outputs from the simulation model were used in this experiment to ensure the best possible parameter values were obtained and thus allow a more reliable comparison of the possible performance of each observer type.

The performance of the multiple-model observers was again compared to a standard Kalman filter labelled `KF3', with a manually-tuned covariance matrix,
\begin{equation} \label{eq:Q_opt}
	\mathbf{Q}_{opt}=\begin{bmatrix}
		0.01^2 & \cdots & 0 & 0 \\
		\vdots & \ddots & \vdots & \vdots \\
		0 & \cdots & 0.01^2 & 0 \\
		0 & \cdots & 0 & 0.0234^2
	\end{bmatrix}.
\end{equation}

$\mathbf{Q}_{opt}$ was found to produce the lowest average error of the output estimates, $\gls{RMSE}(\hat{Y}(N),Y(N))$, when simulated on the data from simulation number 6.

All three observers used the same initial covariance of the state estimates:
\begin{equation} \label{eq:P0}
	\mathbf{P}_0=\begin{bmatrix}
		0.1 & \cdots & 0 & 0 \\
		\vdots & \ddots & \vdots & \vdots \\
		0 & \cdots & 0.1 & 0 \\
		0 & \cdots & 0 & 0.01
	\end{bmatrix}.
\end{equation}

Table \ref{tb:grind1-obs-params-sim1} is a summary of the parameters of the four observers chosen for evaluation.
\begin{table}[ht]
	\begin{center}
		\caption{Observer parameters for grinding process simulation.} \label{tb:grind1-obs-params-sim1}
		% See: https://texblog.org/2019/06/03/control-the-width-of-table-columns-tabular-in-latex/
		\begin{tabular}{p{0.16\textwidth}>{\centering\arraybackslash}p{0.12\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}>{\centering\arraybackslash}p{0.16\textwidth}}
			& KF3 & MKF--SF95 & MKF--SF1 & MKF--SP \\
			\hline
			Type & Kalman filter & Multi-model Kalman filter & Multi-model Kalman filter & Multi-model Kalman filter \\
			Algorithm & - & Sequence fusion 95 & Sequence fusion 98 & Sequence pruning \\
			\hline
			Parameters &  &  &  &  \\
			$\mathcal{Q}$ & $\mathbf{Q}_{opt}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ & $\{\mathbf{Q}_0,\mathbf{Q}_1\}$ \\
			$\mathbf{R}$ & $5^2$ & $5^2$ & $5^2$ & $5^2$ \\
			$\mathbf{P}(0)$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ & $\mathbf{P}_0$ \\
			$n_h$ & 1 & 34 & 26 & 25 \\
			$N_f$ & - & 60 & 60 & - \\
			$n_\text{max}$ & - & 2 & 2 & - \\
			$d$ & - & 10 & 12 & - \\
			$N_\text{min}$ & - & - & - & 23 \\
			$\varepsilon$ & - & 0.01 & 0.01 & 0.01 \\
			$\sigma_{w_p}$ & - & 0.002717 & 0.002717 & 0.002717 \\
			$b$ & - & 100 & 100 & 100 \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

\subsection{Simulation results} \label{sec:grind1-sim-results}

Figure \ref{fig:grind1_rod_obs_est_MKF_SF1} shows the estimates of the sequence fusion observer, MKF--SF1, compared to the Kalman filter and the true system outputs. For clarity, only the estimates of MKF--SF1 are shown here (it had the lowest estimation errors). The estimates of the other two MKF observers were very similar. Note that the data used in this simulation was not used for either system identification or observer parameter selection. From the upper plot it is clear that the measurement noise level is higher than in the previous simulations, relative to the magnitude of the changes in the system outputs. The estimates of both observers lag the true outputs after step changes in the disturbance. It is unclear from this plot whether the multiple-model observer has a significant performance advantage over the Kalman filter. The \gls{MKF} observer estimates seem to fluctuate less during steady-state periods, immediately before the first step change. However, the differences are not as obvious as they were in the linear system simulations.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/grind1_rod_obs_sim_3_est_MKF_SF1.pdf}
	\caption{Observer estimates -- MKF--SF1}
	\label{fig:grind1_rod_obs_est_MKF_SF1}
\end{figure}
Occasionally, rapid changes occur in the estimates of both observers even when there is no disturbance. This may be attributed to extreme values in the measurement noise which occurred at certain times.

As in the previous experiments, longer simulations are needed to evaluate the observers. This was achieved using data from 10 simulations, similar to and including simulation number 6 but generated from unique pseudo-random initializations. These simulations were numbered 6 to 15 as shown in Table \ref{tb:grind1-sims}. Two additional \gls{RMSE} metrics were used in this experiment. As explained in section {\ref{sec:evaluation}}, the {\gls{RMSE}} during transient periods captures the prediction errors when the system was responding to a disturbance. The longest settling time was found to be approximately 1.2 hours when the process was switching from ore mix \#2 to mix \#1. The {\gls{RMSE}} in steady-state captures the errors in sample times occurring more than 1.2 hours after one disturbance and before the next.
\nomenclature{$N_{\pm5\%}$}{settling time of the process to within $\pm5$ percent of steady-state, in number of sampling intervals}%
\begin{table}[ht]
	\begin{center}
		\caption{Observer performance evaluation metrics.} \label{tb:grind1-RMSE-results}
		% See: https://texblog.org/2019/06/03/control-the-width-of-table-columns-tabular-in-latex/
		\begin{tabular}{p{0.29\textwidth}>{\centering\arraybackslash}p{0.08\textwidth}>{\centering\arraybackslash}p{0.14\textwidth}>{\centering\arraybackslash}p{0.12\textwidth}>{\centering\arraybackslash}p{0.12\textwidth}>{\centering\arraybackslash}p{0.08\textwidth}}
			Metric & KF3 & MKF--SF95 & MKF--SF1 & MKF-SP & \gls{SKF} \\
			\hline
			% See script rod_obs_calc_metrics.m
			% 13-Dec-2022 09:56:07 results with system P2DcTd4, sigma_M = 5, tau_ss = 1.2
			$\mathrm{RMSE}(\hat{Y},Y)$ overall & 1.81 & 1.70 & 1.66 & 1.67 & 1.24 \\
			$\mathrm{RMSE}(\hat{Y},Y)$ transient & 2.69 & 2.95 & 2.83 & 2.93 & 2.04 \\
			$\mathrm{RMSE}(\hat{Y},Y)$ steady-state & 1.46 & 1.12 & 1.12 & 1.08 & 0.90 \\
			$\mathrm{Var}(\hat{Y})$ steady-state & 1.61 & 0.66 & 0.70 & 0.58 & 0.25 \\
			$\mathrm{RMSD}(\hat{Y})$ steady-state &0.62 & 0.44 & 0.45 & 0.43 & 0.17 \\
			\hline
		\end{tabular}
	\end{center}
\end{table}
Table \ref{tb:grind1-RMSE-results} shows the values of the five metrics for the five observers evaluated on the data from the 10 longer simulations. In terms of the overall \gls{RMSE}s, the multiple-model observers offer only a marginal improvement compared to the Kalman filter (which was tuned to minimize overall \gls{RMSE}). MKF--SF1 has the lowest overall \gls{RMSE}, which is 8.3\% lower than that of KF3. However, it is 39\% higher than that of the scheduled Kalman filter (\gls{SKF}). The greater difference between the performance of the \gls{SKF} and the other observers in this experiment may be attributed to the discrete time delays in the process model \eqref{eq:grind1-obs_ss_model} (it is impossible for any observer to respond to an input disturbance before its effect has reached the measured outputs, leading to inevitable estimation errors).

The other metrics in Table \ref{tb:grind1-RMSE-results} allow the specific strengths and weaknesses of each observer to be compared. In terms of transient errors, the multiple-model observers are inferior to KF3, with 5.2 to 9.7\% higher \gls{RMSE}s. However, in steady-state periods, they perform significantly better, with 23.3 to 26.0\% lower \gls{RMSE}s. The \gls{RMSD}s of the multiple-model observers are 27.4 to 30.7\% lower, suggesting that insensitivity to measurement noise may be their main advantage. A reduced sensitivity to noise, combined with low transient errors, would be beneficial in feedback control applications where the closed-loop response time is limited by the speed of the observer and where sensitivity to noise causes problems---for example wear or damage to actuators.

Recall that the observers in this evaluation were selected using the overall \gls{RMSE} metric \eqref{eq:rmse-calc-yest}. As mentioned in Section \ref{sec:sim-obs-lin-1-KF-tuning}, this is only one of many possible ways to make a compromise between responsiveness and sensitivity to measurement noise. The multiple-model observers in Table \ref{tb:grind1-RMSE-results} all have higher \gls{RMSE}s in transient periods than the Kalman filter and lower ones in steady-state, which makes it difficult to compare them precisely. Figure \ref{fig:grind1-RMSE-scatter} allows their performance to be compared in both dimensions---transient and steady-state performance. The plot shows the results of all the parameter search simulations in yellow and blue, as well as the final observers chosen for evaluation, which are shown in orange and labelled `min RMSE'. Note that these simulations were run using data from grinding simulation 6 only. Therefore, the \gls{RMSE} values in this plot are not the same as those in Table \ref{tb:grind1-RMSE-results}, which were based on averages of 10 simulations.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/grind1_rod_obs_sim_popt_RMSE_scatter.pdf}
	\caption{\acrshort{RMSE} trade-off comparison}
	\label{fig:grind1-RMSE-scatter}
\end{figure}
The lines connecting the different \gls{KF} tunings give an indication of the \textit{Pareto front} of the performance of a Kalman filter. Any multiple-model result to the left or below the Pareto front of the \gls{KF} is an improvement since it has a lower \gls{RMSE} in transient or steady-state, or both. From this it can be seen that the best multiple-model observers are clearly a significant improvement over any \gls{KF} tuning.

Figure \ref{fig:grind1-RMSE-scatter} also shows that the combined Pareto front of the performance of all observers switches from \gls{KF} for low transient errors and high steady-state errors, to the multiple-model observers for low steady-state errors and higher transient errors. Thus, if low transient errors are more important than low steady-state errors, a Kalman filter is the best solution. However, at a certain point along the Pareto front, it becomes better to use a multiple-model observer and at this point, a significant drop in the steady-state error is achievable. Again, this suggests that multiple-model observers may be advantageous in applications where sensitivity to measurement noise during steady-state periods is important and where the best transient performance is not essential. The plot also reveals that even lower steady-state errors can be achieved using sequence pruning observers with different parameter-settings. Some of these have steady-state \gls{RMSE}s up to 27\% lower than \gls{KF}, although these settings result in high transient errors.

To provide further insight into the transient behaviour of one of the observers, 100 responses of the sequence fusion observer MKF-SF1 to step disturbances were selected from simulations 6 to 15. Step changes were selected only when the system was in steady-state with ore mix \#1 when a step change occurred to mix \#2 and the process subsequently reached steady-state. The response sequences were aligned by subtracting the step time, $t_\text{step}$, from the times of each response.\nomenclature{$t_\text{step}$}{time when a step change in the ore feed occurred in the grinding process simulations} Then the mean, minimum, and maximum values of the observer output estimates were calculated at each sample time since the step change. The results are shown in Figure \ref{fig:sim_resp_plot_MKF_SF1} as a thick orange line representing the average of the 100 responses, and a pale orange area indicating the minima and maxima. Similar data for the Kalman filter, KF3, and the true system output are superimposed in the plot. This allows both the average (i.e. expected) observer response and the range of the variation that occurred to be visualised and compared.
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=12cm]{images/grind1_rod_obs_sim_resp_plot1_MKF_SF95.pdf}
%	\caption{Average observer responses to step disturbances -- MKF--SF95}
%	\label{fig:sim_resp_plot_MKF_SF95}
%\end{figure}
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/grind1_rod_obs_sim_resp_plot1_MKF_SF1.pdf}
	\caption{Average observer responses to step disturbances -- MKF--SF1}
	\label{fig:sim_resp_plot_MKF_SF1}
\end{figure}
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=12cm]{images/grind1_rod_obs_sim_resp_plot1_MKF_SP1.pdf}
%	\caption{Average observer responses to step disturbances -- MKF--SP}
%	\label{fig:sim_resp_plot_MKF_SP}
%\end{figure}
Figure \ref{fig:sim_resp_plot2_MKF_SF1} is a similar plot produced in the same way for reverse steps when the process switched from ore mix \#2 to \#1.
\begin{figure}[htp]
	\centering
	\includegraphics[width=13cm]{images/grind1_rod_obs_sim_resp_plot2_MKF_SF1.pdf}
	\caption{Average observer responses to step disturbances -- MKF--SF1}
	\label{fig:sim_resp_plot2_MKF_SF1}
\end{figure}

Figures \ref{fig:sim_resp_plot_MKF_SF1} and \ref{fig:sim_resp_plot2_MKF_SF1} reveal a few interesting characteristics of the multiple-model observer. Firstly, the average response has a characteristic `S'-shape when compared to that of KF3. It is initially slightly slower to respond, then overtakes the response of KF3 and reaches the vicinity of the true system output significantly earlier. This can be explained by the bi-modal behaviour of the multiple-model observer, which is insensitive to noise most of the time but has a high correction gain once a shock is detected. However, the average response is not representative of individual responses. From the minimum and maximum values it can be seen that during the transition period, the multiple-model observer exhibits high variation and often larger estimation errors than the Kalman filter. Around 1.5 hours after the shocks, the variation in the estimates of MKF--SF1 reduces and becomes less than that of KF3.

\subsection{Sensitivity analysis} \label{sec:grind1-sensitivity-analysis}
Finally, an analysis was carried out to investigate the sensitivity of the observer performance to selected parameters. Firstly, the sensitivity of the overall \gls{RMSE}s to the process model parameters was analysed. The continuous-time linear process model identified from simulation data \eqref{eq:grind1-id-model-ctf} was generalised by introducing two adjustable parameters, a static gain, $K_p$, and a time constant, $T_{p,1}$,
\nomenclature{$T_{p,1}$}{time constant of the process model used in sensitivity analysis}%
\nomenclature{$K_p$}{gain of the process model used in sensitivity analysis}%
\begin{equation} \label{eq:grind1-sens-model-ctf}
	Y(s)= \frac{K_p \exp(-0.2s)}{(1 + T_{p,1})^2}U(s).
\end{equation}

The process delay of 0.2 hours remained fixed. Each parameter was than adjusted according to a set of 7 ratios from 0.5 to 2 (including 1) and the adjusted models were used in the observers and simulated using data set number 6. Since all possible combinations of both parameter values were considered, there were $7\times7=49$ simulations. After each simulation, the ratios of the overall \gls{RMSE} of each observer, to that with the original parameter values were calculated. The results are shown in Figure \ref{fig:grind1-obs-sim-sens-model-y-est} as percentage differences, with a colour-mapping to indicate the magnitude of the changes.
\begin{figure}[htp]
	\centering
	\includegraphics[width=15cm]{images/rod_obs_sim_sens_sys_4obs_RMSE_y_est.pdf}
	\caption{Percentage changes in \acrshort{RMSE}s due to changes in process model parameters}
	\label{fig:grind1-obs-sim-sens-model-y-est}
\end{figure}

From these results, it can be seen that the average magnitude of the estimation errors of the multiple-model observers varied by $-5.2$ to $+11.3$\% over the range of parameter values simulated. The highest errors occurred when the process gain, $K_p$ was doubled. The effect of changes in the time constant, $T_{p,1}$, was less significant.

By comparing these variations to those of KF3, it can be seen that the multiple-model observers are less sensitive to errors in the model parameters. The few cases where the changes in the performance of the multiple-model observers are larger are all cases where the estimation errors were lower than the KF3, therefore better performance. Based on this, there does not appear to be a significant downside to the multiple-model observers when it comes to the accuracy of the process model. In fact, it shows that better overall estimation performance may have been possible with different parameter values.

In the second part of this analysis, the sensitivity of the observer performance to the \gls{RODD} model parameters, $\sigma_{w_p}$ and $\varepsilon$, was investigated. The two parameters were adjusted in the same way as in the previous analysis and the changes in the overall \gls{RMSE}s calculated. The results are shown in Figure \ref{fig:grind1-obs-sim-sens-rod-y-est}. Note that the Kalman filter does not include a \gls{RODD} model and was therefore excluded. Instead, the sensitivity of the scheduled Kalman filter, \gls{SKF}, was included for comparison.
\begin{figure}[htp]
	\centering
	\includegraphics[width=15cm]{images/rod_obs_sim_sens_rod_4obs_RMSE_y_est.pdf}
	\caption{Percentage changes in \acrshort{RMSE}s due to changes in \acrshort{RODD} parameters}
	\label{fig:grind1-obs-sim-sens-rod-y-est}
\end{figure}
Note that the \gls{SKF} has the true shock sequence and therefore does not depend on $\varepsilon$. From these results it can be seen that the average magnitude of the estimation errors of the multiple-model observers varied by $-0.8$ to $+5.8$\% over the range of parameter values simulated. The highest errors occurred when both $\sigma_{w_p}$ and $\varepsilon$ were double their original values.
%\section{Grinding circuit control simulation} \label{sec:sim-ore-mimo-ctrl} 
%
%Outline notes:
%\begin{itemize}
%	\item Grinding simulation model in closed loop with \gls{MPC} controller.
%	\item Diagram of feedback system – Figure \ref{fig:sim-mpc-diag}
%	\item Table of results - Performance metrics — e.g. tracking error.
%	\item Robustness?  E.g. stability margins.
%\end{itemize}
%
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=11cm]{images/sim-mpc-diag.pdf}
%	\caption{Functional diagram of the simulated feedback control system}
%	\label{fig:sim-mpc-diag}
%\end{figure}
%
%\begin{figure}[htp]
%	\centering
%	\includegraphics[width=15.5cm]{images/mpc4x4_stepresp_matrix.pdf}
%	\caption{Step responses of identified input-output model}
%	\label{fig:mpc4x4-stepresp-matrix}
%\end{figure}

\section{Discussion} \label{chap-results-discussion}

In this section, the results presented in sections \ref{sec:sim-obs-lin} and \ref{sec:sim-ore-SISO} are discussed to identify common findings and general conclusions.

\subsection{Observer parameters} \label{sec-disc-obs-params}

By comparing the best observer parameter values for each of the three simulation experiments, it is possible to identify relations between certain parameters and relevant characteristics of the simulated systems. These correlations should align with intuitive expectations based on the design of the algorithms. Selected details of the three simulated systems are listed on the left of Table \ref{tb:param-summary-all-sims}, including the type of process model (linear or non-linear), the input-output dimensions of the process, the probability of the simulated shocks, $\varepsilon$, the \textit{signal-to-noise ratio} (SNR) of the output measurements, and the settling time of the process in number of sampling intervals, $N_{\pm5\%}=t_{\pm 5\%} / T_s$, where $t_{\pm5\%}$ is the time taken for the process output, $y(k)$, to reach within $\pm5\%$ of the steady-state value. The signal-to-noise ratio is defined 
\begin{equation} \label{eq:SNR}
	\text{SNR}_i=\frac{\max{Y_i(N)}-\min{Y_i(N)}}{\sigma_M}.
\end{equation}%
Selected parameter values of the multiple-model observers are included on the right of the table.
\begin{table}[ht]
	\begin{center}
		\caption{Summary of experiments -- observer parameters} \label{tb:param-summary-all-sims}
		\begin{tabular}{
				l
				|
				c
				c
				c
				c
				c
				|
				>{\centering\arraybackslash}p{0.9in}
				>{\centering\arraybackslash}p{0.8in}
				>{\centering\arraybackslash}p{0.8in}
			}
			\multirow{2}{*}{Sec.}
			& \multicolumn{5}{c|}{System}
			& \multicolumn{3}{c}{Observer parameters} \\             
			& Type & Dim. & $\varepsilon$ & \acrshort{SNR} & $N_{\pm5\%}$ & MKF--SF95 & MKF--SF1 & MKF-SP1 \\
			\hline
			\ref{sec:sim-obs-lin-1} & linear & \gls{SISO} & 0.01 & 10 & 9 & 
			$\begin{aligned} n_h & :8 \\ N_f & :5 \\ n_\text{max} & :1 \\ d & :1 \end{aligned}$  & 
			$\begin{aligned} n_h & :6 \\ N_f & :6 \\ n_\text{max} & :1 \\ d & :2 \end{aligned}$  & 
			$\begin{aligned} n_h & :10 \\ N_\text{min} & :7 \end{aligned}$ \\
			\hline
			\ref{sec:sim-obs-lin-2}  &  linear & $2 \times 2$ & 0.01 & 5  & 26  &
			$\begin{aligned} n_h & :116 \\ N_f & :15 \\ n_\text{max} & :2 \\ d & :3 \end{aligned}$  & 
			$\begin{aligned} n_h & :58 \\ N_f & :18 \\ n_\text{max} & :2 \\ d & :5 \end{aligned}$  & 
			$\begin{aligned} n_h & :61 \\ N_\text{min} & :29 \end{aligned}$ \\
			\hline
			\ref{sec:sim-ore-SISO}  &  non-linear  &  \gls{SISO}  &  0.01  &  1.9  &  24 &
			$\begin{aligned} n_h & :34 \\ N_f & :60 \\ n_\text{max} & :2 \\ d & :10 \end{aligned}$ & 
			$\begin{aligned} n_h & :26 \\ N_f & :60 \\ n_\text{max} & :2 \\ d & :12 \end{aligned}$ & 
			$\begin{aligned} n_h & :25 \\ N_\text{min} & :23 \end{aligned}$ \\
			\hline
		\end{tabular}
	\end{center}
\end{table}

By comparing the number of hypotheses required by each observer, $n_h$, it is clear that the $2 \times 2$ (\gls{MIMO}) system simulations in Section \ref{sec:sim-obs-lin-2} required a significantly higher number than the \gls{SISO} system simulations. This makes sense since the number of possible shock combinations is an order of magnitude higher with two independent \gls{RODD}s than it is with one.

Nevertheless, there is a difference in the number of hypotheses required between the linear \gls{SISO} system in Section \ref{sec:sim-obs-lin-1} and the \gls{SISO} non-linear grinding model simulations in Section \ref{sec:sim-ore-SISO}. There are various possible explanations for this difference. It could be the fact that the true system in the latter case is non-linear and there was a mismatch between the prediction model used by the observer and the true system dynamics. However, there is another significant difference between the two systems. The \gls{SNR} was over 5 times less in the grinding simulations. Also, the settling time of the grinding process model was 24 sample intervals compared to only 9 for the linear \gls{SISO} model. Both these attributes increase the \textit{estimation horizon}---the number of samples needed to achieve good estimates (i.e. a low error covariance). When the estimation horizon is long, hypotheses must be maintained longer before being evaluated and merged or pruned. In the case of the sequence fusion observers, MKF--SF95 and MKF--SF1, this is achieved by a longer detection interval. From the values of $d$ in Table \ref{tb:param-summary-all-sims}, it can be seen that this is indeed the case. In the case of the sequence pruning observer, MKF--SP1, a long estimation horizon is achieved by a longer minimum life of the hypotheses, $N_\text{min}$. The results show that this is also the case in these simulations.

Something worth noting about the sequence pruning algorithm is the relationship between the best values of the two parameters, $n_h$ and $N_\text{min}$. In the case of the $2 \times 2$ system, two new hypotheses are generated every sample time, whereas in the case of the \gls{SISO} system only one is. In general, $n_pN_\text{min}$ hypotheses have not reached minimum life at any point in time, where $n_p$ is the number of independent shock variables. Therefore, the number of hypotheses which have reached or exceeded the minimum life is $n_h - n_p N_\text{min}$. This number is 2 or 3 in all three cases, indicating that for some reason, there is no benefit in maintaining more than a small number of hypotheses beyond the minimum life. Thus, the minimum life parameter, $N_\text{min}$, is effectively the sole decision variable, and a good value for $n_h$ may be found using a simple heuristic such as
\begin{equation} \label{eq:nh-sp}
	n_h = n_p N_{\text{min}} + 2.
\end{equation}

One thing that is contrary to expectations is the number of hypotheses needed by the sequence fusion algorithms compared to the sequence pruning algorithm. \cite{robertson_method_1998} stated that the sequence fusion algorithm with detection intervals is designed to reduce the number of hypotheses required when the detection horizon is long, especially when there is more than one \gls{RODD} to distinguish between. In both the $2 \times 2$ simulation and the grinding process simulations, MKF--SP1 requires about the same number as MKF--SF1. However, the number of hypotheses used by MKF--SF1 is significantly less than that of MKF--SF95, so it is an improvement in this respect.

The shock probability, $\varepsilon$, was the same in all experiments, so its effect on the observer parameters cannot be inferred from these results.

\subsection{Advantages and disadvantages} \label{sec-disc-obs-performance}

The relative performance of the multiple-model observers compared to the best-tuned Kalman filter varied between the three experiments. Table \ref{tb:rmse-summary-all-sims} shows the percent differences of the overall RMSEs of the multiple model observers, compared to those of the Kalman filters, labelled KF3. The percent differences were calculated as
\begin{equation} \label{eq:pct-diff}
	\text{\% difference} = 100 \left( \frac{\text{RMSE}_\text{MKF}}{\text{RMSE}_\text{KF3}} - 1 \right).
\end{equation}
\begin{table}[ht]
	\begin{center}
		\caption{Summary of experiments -- observer performance} \label{tb:rmse-summary-all-sims}
		\begin{tabular}{
				l
				|
				c
				c
				c
				c
				c 
				|
				>{\centering\arraybackslash}p{0.9in}
				>{\centering\arraybackslash}p{0.8in}
				>{\centering\arraybackslash}p{0.8in}
			}
			\multirow{2}{*}{Sec.}
			& \multicolumn{5}{c|}{System}
			& \multicolumn{3}{c}{\% difference in RMSEs vs. KF3} \\             
			& Type & Dim. & $\varepsilon$ & \acrshort{SNR} & $N_{\pm5\%}$ & MKF--SF95 & MKF--SF1 & MKF-SP1 \\
			\hline
			\ref{sec:sim-obs-lin-1} & linear & \gls{SISO}              & 0.01 & 10 & 9 & -35.7 & -35.5 & -36.1 \\
			\ref{sec:sim-obs-lin-2} & linear & $2 \times 2$           & 0.01 & 5 & 26 & -21.2 & -21.2 & -18.4 \\ 
			\ref{sec:sim-ore-SISO} & non-linear & \gls{SISO} & 0.01 & 1.9 & 24 & -6.0 & -8.3 & -7.4 \\      
			\hline
		\end{tabular}
	\end{center}
\end{table}

The differences in overall performance between the multiple-model observers (\acrshort{MKF}s) in each simulation experiment are quite small. The sequence pruning observer was slightly better in the first simulation. In the other two, the sequence fusion observers were better. However, the differences between the \gls{MKF}s and the Kalman filters (KF3) are significant and depend on the simulation. In the first simulation, the \gls{MKF}s have approximately 36\% lower overall \gls{RMSE}s, whereas in the third experiment, they achieve only a 6 to 8\% reduction. This can be partly explained by the delay in the response of the process in this experiment. When the system has a delay between the unmeasured disturbance and the output, estimation errors are unavoidable. The greater the total error, the less the relative reduction due to any improvement in state estimation is. Another factor is the longer estimation horizon needed by the multiple model observers in the third experiment, as explained in Section \ref{sec-disc-obs-params}. Estimation errors are higher when the response of the observers is slower. Although KF3 is not ideal, it has a faster response than a Kalman filter tuned to the persistent disturbance (KF1 in sections \ref{sec:sim-obs-lin-1} and \ref{sec:sim-obs-lin-2}). The longer the detection interval, the less advantage the MKFs have over the response of KF3. As discussed in Section \ref{sec:grind1-sim-results}, the main advantages of the \gls{MKF}s were found to be in steady-state, when the \gls{RMSE}s were up to 30\% less than those of KF3.

While the multiple-model observers have a number of clear advantages over a single Kalman filter, their disadvantages are important to recognize. The complexity of these algorithms is a significant impediment from the perspective of practical implementation. The computer code to implement each of them is extensive compared to that of a single Kalman filter with a time-varying gain. To justify the higher implementation costs, the benefit of any improvement in state estimation would need to be quantified, and this would depend on the intended application. Finally, the MKFs evaluated here are of course limited to the infrequently-occurring disturbances that can be represented by the \gls{RODD} model. However, this is a versatile model capable of representing a wide range of possible disturbances.

In this chapter, the results of numerical simulations to evaluate the disturbance models and observer designs were presented and discussed. The main conclusions, limitations of the work, and recommendations for future research are described in the next section.