%===============================================================================
% $Id: ifacconf.tex 19 2011-10-27 09:32:13Z jpuente $  
% Template for IFAC meeting papers
% Copyright (c) 2007-2008 International Federation of Automatic Control
%===============================================================================
\documentclass{ifacconf}

\usepackage{graphicx}      % include this line if your document contains figures
\usepackage{natbib}        % required for bibliography
\usepackage{epsfig} % for postscript graphics files\usepackage{subcaption}

%===============================================================================
\begin{document}
\begin{frontmatter}

\title{Fusing Physical Process Models with Measurement Data Using FIR Calibration\thanksref{footnoteinfo}} 
% Title, preferably not more than 10 words.

\thanks[footnoteinfo]{The work in this paper was partly funded by the H2020 project COGNITWIN (grant number 870130). }

\author[First]{Enso Ikonen} 
\author[First]{Istv\'an Selek} 

\address[First]{ Intelligent Machines and Systems,
        University of Oulu, FIN-90014 Oulun yliopisto, Finland (e-mail: Enso.Ikonen@oulu.fi; Istvan.Selek@oulu.fi).}
%\address[Second]{ Intelligent Machines and Systems,
%        University of Oulu, FIN-90014 Oulun yliopisto, Finland (e-mail: Istvan.Selek@oulu.fi).}
%\address[Third]{Electrical Engineering Department, 
%   Seoul National University, Seoul, Korea, (e-mail: author@snu.ac.kr)}

\begin{abstract}                % Abstract of not more than 250 words.

Tuning of physical plant models is considered in the context of monitoring and control of industrial processes. The problem of calibrating physical models is discussed. A method is proposed, consisting extending the physical model with finite impulse response (FIR) filters at the output. The approach is illustrated in extensive simulations and applied in state estimation using a nonlinear chemical continuous stirred tank reactor benchmark. 

\end{abstract}

\begin{keyword}
dynamic models; machine learning; physical models; process control; process models; state estimation 
\end{keyword}

\end{frontmatter}
%===============================================================================

\section{INTRODUCTION}

Modern design of process control and monitoring relies heavily on dynamic simulation models. On-line simulations can be necessary for evaluating model predictive controllers, model-based state-estimators, or in model-based fault detection and isolation. Optimization of plant short term dynamic operation or integrated plant and control design are other examples of applications requiring a proper dynamic process simulator. A modern digital twin gets is essential behavioral contents from a plant model, linking  simulations with on-line data and various services such as data validation, visualization, process analytics, what-if analysis and 3D plant animations, for example. In this paper, the focus is on process monitoring and control. 

In the era of IoT (Internet of Things) and ML/AI (Machine Learning / Artificial Intelligence), the role of data-driven modeling has been emphasized. A vast literature exists on proposed data-driven approaches to process modeling as well as experiences on their application to simulated and/or real-life plants (\cite{Sjoberg}; \cite{HastieESLII}; \cite{APIC}). The approaches include various basis-function constructions proposed suitable for process modelling tasks, statistical analytics on conclusions that can be drawn from data, and algorithms providing optimal active exploration, just to mention few. Data-driven modeling is plagued by the bias-variance dilemma, related to supression of noise in measurements and extracting a proper representation of the underlying nonlinear mapping. Both problems are impossible to solve using data alone. The impact of noise can be reduced by taking repeated samples, but the uncertainty will remain. If only sampled data is availabe, there is no information about the mapping between samples, in time nor space. Consequently, interpolation/extrapolation, i.e. generalization, is always based on some additional assumptions outside of the data. 

The application of physical models is widespread in many fields of engineering (\cite{Skogestad}). Physical models provide justified and transparent predictions and explanations for the outcomes, based on mass, momentum and energy balances, laws of mechanics, chemistry, etc. The role of physical models is pronounced in the heavy process industry, where models are extensively used in plant design and experimentation and exploration with plants is severly complicated in the industrial practice, due to costs, slowness of processes and safety issues. However, construction of physical models may also require significant resources and some tuning of such models to local site conditions is always necessary. 

The essential component of ML/AI is that of learning from examples, such as process measurement data. Therefore, the topics of model tuning and learning are closely entwined. A recent trend in process control literature is to fuse machine learning with physical models (\cite{Bikmukhametov}; \cite{Qin}), with a recognized aim to simultaneously improve the automatic adaptation of physical process models to local conditions and the explainability of data-driven models. This is also the topic of this paper. Compared to existing works, our focus is in this paper to maximally exploit the knowledge in physical models/simulators, and to emphasize the importance of system dynamics. Less weight is put on validating new inventive ML/AI-driven methodologies, as simplicity often comes with robustness and ease of implementation.

This paper considers the problem of identifying an I/O model for a MIMO system, such as
\begin{equation}
\begin{array}{ll}
\mathbf{x}\left(k+1\right)=f\left(\mathbf{x}\left(k\right),\mathbf{u}\left(k\right),\mathbf{w}\left(k\right)\right) \\
\mathbf{y}\left(k\right)=g\left(\mathbf{x}\left(k\right),\mathbf{v}\left(k\right)\right)
\end{array}
\end{equation}
where $f:R^{n_x}\times R^{n_u}\times R^{n_{w}}\rightarrow R^{n_x}$ is the state propagation function and $g:R^{n_x}\times R^{n_{v}}\rightarrow R^{n_y}$ the measurement function. For brevity, the sampling time is assumed to be a constant and the important topics of input/output selection and design of sampling are not considered here. The inputs $\mathbf{u}$ are assumed to be known and measured, possible uncertainties in inputs can be included in the random state noise $\mathbf{w}$. States $\mathbf{x}$ need not be known, measured or estimated. The plant outputs are measured, corrupted by random noise $\mathbf{v}$. Noise is assumed to have finite variance, smoothness or existance of derivatives of $f$ or $g$ is not assumed. It is assumed that the system is time-invariant, or that it varies with time very slowly so that (some kind of) ergodic assumptions hold. 

This paper proposes an approach based on the use of physical modeling and FIR (Finite Impulse Response) filters. The main benefits are in that the plant models (simulators) can be incorporated into the proposed structure and that the parameter estimation in the data-driven components is particularily simple and robust, enabling automatic procedures. The few parameters to be set by the user are very transparent and their tuning does not require an experienced process control engineer. Section 2 considers calibration of physical plant models and leads to propose a PMFIR structure for model calibaration in Section 3. Applications in state estimation and control are briefly discussed in Section 4. Section 5 illustrates and approach using the well known van der Vusse CSTR benchmark. Discussion and conclusions end the paper.

\section{Calibration of physical models}

Suppose that a given plant has a dynamic I/O model, such as a state space model:
\begin{equation}
\begin{array}{ll}
\mathbf{x}\left(k+1\right)=f_{PM}\left(\mathbf{x}\left(k\right),\mathbf{u}\left(k\right)\right)
\\
\mathbf{y}_{PM}\left(k\right)=g_{PM}\left(\mathbf{x}\left(k\right)\right)
\label{eq:PM}
\end{array}
\end{equation}
$f_{PM}:R^{n_x}\times R^{n_u}\rightarrow R^{n_x}$ and $g_{PM}:R^{n_x}\rightarrow R^{n_y}$.
This model, called a ``physical model'', PM, is likely to have been derived based on 1st principles or detailed engineering knowledge. It may be complex e.g. due to size or amount of detail, and nothing can be assumed about its internal structure or parameters. It is convenient to simulate the propagation of $\mathbf{y}_{PM}$ in real time in parallel with the plant, or using
historical data. The plant model is assumed to have been calibrated as found appropriate and feasible given the constraints on knowledge and resources. It is provided ``as is''. This kind of models commonly emerge in industrial plant and/or control design. 

Suppose also, that corresponding I/O measurements from a real plant are available: $\mathbf{u}\left(i\right)$ and $\mathbf{y}\left(i\right)$, $i=1,2,...,k$. This data is to be used to calibrate the process model, so that the predictions would be consistent with the data. The need for tuning may be due to various reasons, such as correction of simplifications/idealisations made during modeling or adaptation to systematic errors in plant measurements. The purpose of adjustments is to match the physical model and measurements from the local site, in this sense the true state of the plant is irrelevant. 

This wish to calibrate the physical model is constrained by a number of engineering requirements. i) Very often it is not desired tune the parameters inside $f_{PM}$ or $g_{PM}$, in order to preserve the original model. This is due to various reasons: these parameters may appear in a complex manner; the plant model is trusted/accepted by the users, or the users are familiar with the model with its defects. ii) The calibration must be automatic from data and should scale well. This requires needing minimal insight to the physical model, and robustness vs. measurement noise and structural errors. iii) The algorithm should be simple to implement and execute on-line in a distributed control system (DCS) environment or alike. iv) The approach should support the a priori assumption that industrial process models are often more accurate in their steady-state description, typically in the focus of plant design. v) The purpose of calibration is not to model phenomena/behaviour that is not described by the physical model. A more detailed plant modeling (of any unmodelled behaviour) is almost surely a project requiring an approach different from an automated procedure. However, it is clear that the borderline between modeling and tuning is somewhat vague. Therefore, vi) the procedure should give some means to analyze the adjustments made to the physical model mapping due to calibration. 

Joint physical/data-driven structures (see also \cite{Bikmukhametov}) for modeling dynamic processes can be roughly categorized
into the following:
\begin{itemize}
\item A data-driven model (DM) is placed in parallel with the physical model (PM). The outcomes of the two models can be compared by running models in like situations, and the multimodel information then presented to the end-users. Typically, data-driven approaches
can provide an estimate of the prediction uncertainty. If such uncertainty estimate is available also for the PM (e.g. via Monte Carlo), some sort of bayesian reasoning can be conducted to devise a more likely estimate.  
\item A DM is trained to minimize the residual between measurements and the PM output. The overall outcome is a superposition of the PM and DM. This approach is appealing in that the DM is used to compensate for the defects in PM predictions, the behaviour already modeled by PM need not be identified by DM. 
\item The PM is placed in series with the DM. The DM can precede or follow the PM. DM optimization aims at  tuning both static and dynamic behaviour of the PM input/output signals, minimizing the residual between the final prediction and output data. In the linear case this leads to the well known Wiener/Hammerstein stuctures, common in the industry due to their robustness (\cite{ProcME}).
\end{itemize}
A basic analysis of some of the differences of the structures is simple to derive in the linear case (\cite{ANZCC}). Suppose that the physical model dynamics are given by $\frac{num_{PM}}{den_{PM}}$ and the plant data is obtained from from $\frac{num_{D}}{den_{D}}$.
Let the DM minimize the difference between DM output and measured plant D output. The first approach leads to estimate $\frac{num_{DM}}{den_{DM}}\approx\frac{num_{D}}{den_{D}}$, while the second and third lead to estimate more complicated dynamics: $\frac{num_{DM}}{den_{DM}}\approx\frac{num_{D}den_{PM}-num_{PM}den_{D}}{den_{D}den_{PM}}$ and $\frac{num_{DM}}{den_{DM}}\approx\frac{num_{D}}{den_{D}}\frac{den_{PM}}{num_{PM}}$, respectively. However, the first approach
provides no remedy for the PM calibration problem, it merely provides a data-driven alternative DM. 

An alternative to the above ARX $\frac{num_{DM}}{den_{DM}}$ (autoregressive exogeneous) is to estimate only the $num_{DM}$ (finite impulse response, FIR) (\cite{APIC}):
\begin{equation}
\begin{array}{ll}
y\left(k\right)=b_{0}u\left(k\right)+b_{1}u\left(k-1\right)+...b_{n-1}u\left(k-n-1\right)
\end{array}
\end{equation}
The FIR parameters can be estimated using the least squares method (LS), a number of extensions have been proposed to handle the nonparsimonious nature of FIR (\cite{Nikolaou}). The FIR filters are known to be able to model complex dynamics. However, this comes at the price of a large number of parameters to estimate: for FIR the $dim(num)$ easily becomes much larger than $dim(num)+dim(den)$ for ARX. However, the number of parameters in ARX modeling can also be higher than expected at first glance, as was illustrated above. If, in addition, the models are to be used for predicting to the future, the ARX modeling is not feasible, but one needs to turn to output error (OE) modeling. Unfortunately, OE models can not be solved with LS, but require iterative approaches (RPEM, or alike). In conclusion, FIR modeling can be seen as an attractive structure for the problem at hand. This is further confirmed by the fact that a number of model-predictive control approaches use the FIR or finite step reseponse descriptions for the plant (\cite{APIC}). 

Many algorithms have been proposed for for data-driven modeling of nonlinear dynamic processes (\cite{Sjoberg}; \cite{APIC}; \cite{HastieESLII}). Most -if not all- model structures can be interpreted as basis-function approaches. 
A physical model provides one possible way to generate basis functions (e.g. in PM followed by DM structure). In fact, since the PM states $\mathbf{x}$ are often accessible as well, they could be used as a more extensive set of basis functions. This is likely to lead to collinearity problems, but they are routinely solved in machine learning (\cite{HastieESLII}). 
 However, this line of thinking is not now pursued further here but a simpler approach is proposed instead, where the PM is expected to provide the essential nonlinear gains/dynamics of the plant via its output predictions. This basic behaviour is then tuned by adjusting the parameters of additional linear FIR filters.

\section{PMFIR approach}


\begin{figure}
\begin{center}
\includegraphics[width=6 cm]{FigB.jpg}
\caption{PMFIR structure for physical model calibration using FIR filters.}
\label{fig:FigPMFIRStruct}
\end{center}
\end{figure}

The above discussion leads to consider the two-phase structure depicted in Fig. 1. The PMFIR model is evaluated in three phases. First, the physical model PM is simulated using the given input data (either on-line or a series of batch data), providing the PM predictions $\mathbf{y}_\mathrm{PM}$. In a second phase the CFIR filters, $z_i$, are evaluated, and the EFIR filters, $z_{\mathrm{u},i}$, are evaluated in parallel. Finally, the outputs of the two filters are summed, for each output, to provide calibrated the PMFIR predictions, $\hat{y}_i$.

The parameters for the CFIR and EFIR need to be determined by the proposed method. The parameter estimation starts by first evaluating the PM predictions $\mathbf{y}_\mathrm{PM}$ for the training data. Using these as inputs to the CFIR, the optimal CFIR parameters are solved by minimizing the squared residual between training data and CFIR predictions, $y_i-z_i$. The EFIR parameters are solved by minimizing the square of the remaining residual, $(y_i-z_i)-z_{\mathrm{u},i}$.

Next subsections detail the parameter estimation algorithm for batch data using LS. Algorithms for recursive LS with various forgetting schemes are well known and straightforward to devise.  

\subsection{Calibration FIR - CFIR}

Since the physical model $(f_{PM}, g_{PM})$ is to be evaluated ``as is'' using the known inputs $\mathbf{u}\left(k\right)$, the physical model is evaluated first. The physical model outputs $y_{PM}$ are calibrated by single-input single-output (SISO) FIR filters at each output of the PM
\begin{equation}
\begin{array}{ll}
z\left(k\right)=b_{0}+b_{1}y_{PM}\left(k\right)+...+b_{n}y_{PM}\left(k-n+1\right)
\end{array}
\end{equation}
where $y_{PM}\left(k\right)$ is one of the $n_y$ outputs in the vector $\mathbf{y}_{PM}$ (subscripting is omitted for simplicity). Note that the filter includes a bias term $b_{0}$. This, together with the fact the FIR gain is not constrained, implies that the output can be scaled and shifted in an affine manner.

Each filter has $n+1$ coefficients to estimate. If the physical model static performance is fully trusted (not to be calibrated), the steady state gain of the FIR can be constrained to 1. Since FIR parameter estimation is parsimonious (\cite{Nikolaou}), it can be justified to consider regularized approaches when $n$ is considerably large. The goal is then to minimize $\left(\mathbf{Y}-\Phi\theta\right)'\left(\mathbf{Y}-\Phi\theta\right)$ $+r\theta'\mathbf{Q}\theta$
where $r\geq0$ and $\mathbf{Q}$ is a positive definite matrix. The solution is given by
\begin{equation}
\begin{array}{ll}
\hat{\theta}=\left[\Phi_{0}'\Phi_{0}+r\mathbf{Q}\right]^{-1}\Phi_{0}'\mathbf{Y}_{0}
\end{array}
\end{equation}
An increasing weighing for deviations from zero for increasingly delayed elements is obtained e.g. by choosing $\mathbf{Q}=\tilde{N}\mathrm{diag}\left(\left[1,2,...,n\right]\right)$ (in the SISO case with no $b_{0}$) and $r$ is some small number (e.g., $r=0.1$). A recursive parameter estimation algorithm can be devised, e.g., from the Kalman filter.  

\subsection{Extended FIR - EFIR}

The CFIR-calibrated PM can be further extended with a linear
FIR model for the remaining residuals (EFIR, E for extended). Consider the following multiple-input MISO FIR structure
with filter length $n$ for each input ${u_1,u_2,...,u_p}$:
\begin{equation}
\begin{array}{ll}
z_{u}\left(k\right)=b_{u,0}+b_{u,1}u_{1}\left(k-1\right)+...+b_{u,n}u_{1}\left(k-n\right)+\\
\;\;\;\;\;\;\;\;\;\;\;\;...+b_{u,np}u_{p}\left(k-n\right) \nonumber
\end{array}
\end{equation}
Note that the model causality is preserved by the one sample delay in the
inputs. The batch and recursive parameter estimation can be conducted as discussed above, with the
exception that the targeted output is the residual error $y\left(k\right)-z\left(k\right)$,
instead of the measured output $y\left(k\right)$ (multiple-output
indexes omitted for brevity). 

\subsection{Physical model FIR - PMFIR}

The final calibrated PM outputs are then given by
\begin{equation}
\begin{array}{ll}
\hat{y}_{i}\left(k\right)=z_{i}\left(k\right)+z_{u,i}\left(k\right)
\end{array}
\end{equation}
$i=1,2,...,m$, where $m$ is the number of calibrated outputs.

The length of FIR filters $n$ is to be chosen by the modeller. This parameter can be set based on knowledge of plant settling time (e.g., from plant step responses), or optimized, e.g. by minimizing the Akaike Information Criterion (AIC)
\begin{equation}
\begin{array}{ll}
AIC = \sum_{k=1}^K(y(k)-\hat{y}(k))^2 + 2n' + \frac{2n'(n'+1)}{(K-n'-1)}
\end{array}
\end{equation}
where $n'$ is the number of parameters to estimate. 

The two-phase estimation emphasizes the PM-calibration aspect of the approach. All nonlinearities in PMFIR originate from the physical model PM. The CFIR provides a linear dynamic adjustment for the PM outputs. The EFIR extends the tuning by a linear dynamic mapping from past PM inputs, thereby enabling identification of linear components not modelled by the PM. If necessary for increased numerical precision, it would seem natural to devise further versions by replacing the EFIR (or CFIR) by a nonlinear time-series estimator. However, as discussed in previous Sections, the interpretability of the adjustements in such approaches may easily suffer.

\section{Process state estimation and control}

Model-based state estimation and control is regularily used in the advanced industrial practice. While it is not reasonable to discuss process control in detail in the context of this paper, it suffices to note that state estimation is an essential part of state feedback control, whether linear or nonlinear. Of course, monitoring applications are justified by their own right in the industrial practice, and can make great use of state estimation techniques.

Bayesian state estimation is a commonly used tool for estimating system states, for fusion of dynamic plant models and noisy measurements. The famous Kalman filter (\cite{Simon}) assumes a linear plant desciption. The extended Kalman filter (EKF) (\cite{Simon}) is based on on-line linearization via computation of Jacobians. The outcomes are optimal in an environment corrupted by Gaussian noise, but the Kalman filter provides the optimal linear filter also in a non-Gaussian context (see \cite{Simon}). Recent advances in state estimation can make a more direct use of nonlinear physical models. Despite of the advantages of particle filter (PF) (\cite{Gordon}\cite{Simon}), such as the ability to deal with complex unimodal densities, the computational load associated with PF can make it infeasible in practice. The unscented Kalman filter (UKF) (\cite{Julier}; \cite{Simon}) provides an approximation of the mean and covariance for nonlinear systems which is often sufficient and reasonable for practical purposes. The estimation can also be constrained so as to ensure feasible estimates at all times (\cite{Spivey}). 

The UKF is particularily useful in the PMFIR setup where the PM is given in a state space form, eq. (\ref{eq:PM}), as estimation of unknown states in PM is feasible. As one of the options, the calibrated model provides a means to back-calculate the equivalent measurements at the PM level. For the measurement correponding to $y_{PM}$ it is straightforward to derive
\begin{equation}
\begin{array}{ll}
q_{PM}\left(k\right)=\frac{1}{b_{1}}\left[y\left(k\right)-z_{u}\left(k\right)-z\left(k\right)\right]+y_{PM}\left(k\right)
\end{array}
\end{equation}
for the SISO case (indexes omitted for brevity). Note that this simple method requires that $b_1>0$, which may not always be the case. This approach aims at removing the impact of the calibration of the the PMFIR innovation signal. A Kalman filter (extended, unscented, or particle filter) type of estimator can be constructed on the basis of predictions by the process model $\mathbf{y}_{PM}$ and the back-calculated measurement $\mathbf{q}_{PM}$. An estimate of the states of PM  will then be adjusted with the available measurement, consistent with both the calibrated simulation and the physical model internal construction. 

The UKF algorithm is well known, and it is not repeated here (see e.g. \cite{Simon}). It is worthwhile to point out, however, that the one-step ahead simulations from given initial states (determined by the UKF sigma points) are straightforward to propagate and measure, given access to $(f_{PM}, g_{PM})$. The number of model evaluations required by UKF is relative to the number of states $n_x$ (not exponential) and the $2n_x$ one-step ahead simulations can be expected to remain feasible in many cases of industrial size.


\section{Numerical  CSTR benchmark example}

The approach was tested using the well known process control design benchmark of a highly
nonlinear CSTR with a cooling jacket (\cite{Chen}; \cite{CCE16}). The system exhibits characteristics
such as change of sign of gain in steady state behaviour and changes in zero
dynamics.

%$A\overset{k_{1}}{\rightarrow}B\overset{k_{2}}{\rightarrow}C,\:2A\overset{k_{3}}{\rightarrow}D$

The main reaction scheme is the van der Vusse reaction $A\rightarrow B\rightarrow C,\:2A\rightarrow D$. The control inputs are the normalized flow to the reactor $\dot{V}/V_{R}$ and heat $\dot{Q}_{K}$ withdrawn from the coolant by an external
heat exhanger. System dynamics are given by nonlinear ordinary differential equations obtained from component balances for the substances A and B, concentrations $c_{A}$ and $c_{B}$ [mol/l]; and energy balances for the reactor and cooling jacket in temperatures $\nu$ and $\nu_{K}$ [$^o$C]. The product concentration $c_{B}$ and reactor temperature $\nu$ are measured as outputs. A complete description of the equations and parameters can be found from \cite{Chen}. The system description is given by four ODE (\cite{Chen}; \cite{CCE16}):


\begin{equation}
\begin{array}{ll}
\dot{c}_{\textrm{A}} = \frac{\dot{V}}{V_{\textrm{R}}}\left(c_{\textrm{A0}}-c_{\textrm{A}}\right)-k_{1}c_{\textrm{A}}-k_{3}c_{\textrm{A}}^{2}\nonumber \\
\dot{c}_{\textrm{B}} = \frac{\dot{V}}{V_{\textrm{R}}}c_{\textrm{B}}+k_{1}c_{\textrm{A}}-k_{2}c_{\textrm{B}}  \\
\dot{v} = \frac{\dot{V}}{V_{\textrm{R}}} \left(v_{0}-v\right) \label{eq:VDVModelODE} \nonumber \\
-\frac{1}{\rho C_{\textrm{p}}}  \left( k_{1} c_{\textrm{A}}\triangle H_{R_{\textrm{AB}}} 
+k_{2}c_{\textrm{B}}\triangle H_{R_{\textrm{BC}}} 
+k_{3}c_{\textrm{A}}^{2}\triangle H_{R_{\textrm{AD}}}\right) \nonumber \\
+\frac{k_{\textrm{w}}A_{\textrm{R}}}{\rho C_{\textrm{p}}V_{\textrm{R}}}\left(v_{\textrm{K}}-v\right)\nonumber \\
\dot{v}_{\textrm{K}} = \frac{1}{m_{\textrm{K}}C_{\textrm{PK}}}\left(\dot{Q}_{\textrm{K}}+k_{\textrm{w}}A_{\textrm{R}}\left(v-v_{\textrm{K}}\right)\right)\nonumber 
\end{array}
\end{equation}
The benchmark uses a sampling interval of 20 seconds. 

All plant measurement data was generated using the basic parameter values (see \cite{Chen}, Table 1). The simulation inputs consisted of random ramps with random set point changes at random times at both inputs. The start of a ramp was a random event taking place with a probabilty of 20 seconds/5 hours. The length of the ramp ranged from 20 seconds to 20 minutes; the amplitude of the ramps ranged between minimum and maximum values, U(3,35) for  $\dot{V}/V_{R}$ and U(-9000,0) for $\dot{Q}_{K}$. The feed temperature was considered to be a constant at $\nu_0=104.9$ [$^o$C].  Data sets of length two days, two months and one year were considered for parameter estimation. A separate data set of 1 year was generated for testing purposes. In most of the experiments, the data was resampled to a sampling time of 2 min, with no averaging. Both output measurements were corrupted by additive zero mean normal noise with standard deviation corresponding to approximately 1\% of the range of the signal, $\sigma_{c_{B}}=0.01$ and $\sigma_{\nu}=1$, respectively. In all, The input data sequence can be considered as rich for an industrial plant, but does represent a possible scenario. 

The physical model PM was taken to be the worst case scenario 1 (see \cite{Chen}, Table 2), significantly different from the basic values. In some of the experiments the scenario 2 was considered (see \cite{Chen}, Table 2). The steady-state solution for the worst case 1 was considered as the PM in one setup. A FIR filter length of 20 min for $c_B$  and 2 hours for $\nu$ samples was used for all filters, leading to FIR orders $n=(10, 60)$ for $c_B$ and $\nu$, respectively. In some of the exeriments the AIC optima were used, $n=(1,38)$ for the 2 month data set and 2 min sampling time. With steady state PM, slightly longer windows $n=(30,90)$ were used. The set up with: PM equal to 2 month data from the plant sampled every 2 min, PM as the model for the worst case scenario 1, and filter orders $n=(10,60)$ will be referred to as the `nominal case'.

The number of PMFIR parameters to estimate is $1+np+(1+n)$, for each output. In the nominal set up, this results in 32 parameters for $c_B$  and 182 for $\nu$. All estimations were conducted by first normalizing the data to be zero mean with unit variance and using then a regularized LS with $r=0.1$ and weight equal to the delay in input (in samples). This choice of algorithm, hinted in \cite{Nikolaou}, together with the resampling from 20 sec to 2 min, resulted in a relatively smooth step responses in the estimates and very feasible computing times. Robust estimation of FIR/FSR models in the context of industrial procesesses is a rich topic in itself but not the content of this paper.

\subsection{Nominal case}

The nominal data consisted of 43200 sampled $(\mathbf{u},\mathbf{y})$ pairs of inputs $\dot{V}/V_{R}$ and $\dot{Q}_{K}$ and outputs $c_B$ and $\nu$, simulated with parameters from the nominal scenario. The plant model was taken to be the worst case scenario 1, and used for simulating the same input sequence, resulting in the PM prediction. The CFIR and EFIR parameters were then estimated using the data set. The Bode diagram (gain and phase) for the CFIR filters is shown in Fig \ref{fig:FigBode}. The correction of PM for $c_B$ involves a very slight lead, while for $\nu$ the filter is a low pass. The evaluation of the impact of CFIR tuning is immediate and precise via a Bode diagram. Similar plots can be drawn for EFIR's. 

Figures \ref{fig:FigFIRcoeffu1}--\ref{fig:FigFIRcoeffu2} illustrate an alternative view, by showing the unit step responses for the filters. For CFIR, the unit step is the reference for no tuning impact; for EFIR's a zero step indicates no impact. Again, the analysis of tuning outcomes is straightforward, albeit the impact of the scales of the input signals needs to be taken into account when assessing the unit steps. One can also analyze the estimated static affine mapping of the tuning filters by extracting $b_0$ and $b_{u_i,0}$ (shift) and summing the involved FIR coefficients. In the nominal case $y = 0.032+0.99y_p$ was obtained for $c_B$ and $y = -1.191+1.035y_p$ for $\nu$. The static gains were close to 1, as was already visible from the Bode, Fig \ref{fig:FigBode}, but these corrections in gains were significant. The shifts, over 0.03 mol/l in product concentration and one degree Celsius in reactor temperature, were also significant. Similar static summary info can be computed for EFIR's. 

Figure \ref{fig:FigPred} illustrates the predictions for  $\nu$ in a simulation over the test set. The measured data is noisy (standard deviation of 1 $^o$C). The shift error in PM is clearly visible, and it is corrected by the PMFIR. Corrections in the dynamic response are not visible by eye in the one-day simulation interval, these are better quantified from the Bode, Fig. \ref{fig:FigBode}.

\begin{figure}
\begin{center}
\includegraphics[width=8.4cm]{FigFIR_Bode.eps}
\caption{Bode diagrams for estimated CFIR filters.}
\label{fig:FigBode}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=8.4cm]{FigFIR_coeffu1.eps}
\caption{Estimated PMFIR coefficients for output $c_B$.}
\label{fig:FigFIRcoeffu1}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=8.4cm]{FigFIR_coeffu2.eps}
\caption{Estimated PMFIR coefficients for output $\nu$.}
\label{fig:FigFIRcoeffu2}
\end{center}
\end{figure}


\begin{figure}
\begin{center}
\includegraphics[width=8.4cm]{FigPred.eps}
\caption{Predictions of $\nu$ on a short interval of the test set.}
\label{fig:FigPred}
\end{center}
\end{figure}


\subsection{Further analysis}

A series of tests was conducted to better illustrate and quantify the impact of various components in the PMFIR approach and to compare it with alternative approaches. Starting from the nominal case, changes were made in the amount of data samples, sampling rate, PM and filter orders. The approach was also compared with a standard sigmoid neural network in a similar FIR time series setup. As a main metric, RMSE (root mean squared errors) were computed both on training and test data.

The test series consisted of 8 variations to the nominal setup, and three SNN topologies:
\begin{itemize}
\item nominal case: parameter estimation was conducted based on a 2 month data set with 2 min sampling interval. Physical model corresponds to scenario 1 (PM1).
\item nominal case, but using a model corresponding to scenario 2 as a physical model (PM2)
\item nominal case, but using a short data set of 2 days, only.
\item nominal case, but using a short data set and physical model of scenario 2.
\item nominal case, but using a long data set of 1 year data.
\item nominal case, but using a long data set and physical model of scenario 2.
\item nominal case, but using the AIC-optimal orders for the FIR-filters $n=(1,38)$
\item nominal case, but using a samping time of 20 sec and AIC optimal order for FIR-filters $n=(1,59)$
\item nominal case, but using  an approximate static physical model for scenario 1 (PM3) and a higher order of filters $n=(30,90)$
\item sigmoid neural network (SNN) with 8 nodes in the hidden layer, trained using 2 day data
\item SNN with 15 nodes in the hidden layer, trained using 2 month data
\item SNN with 15 nodes in the hidden layer, trained using 1 year data
\end{itemize}

For each of the above PMFIR setups, four RMSE were computed: 
\begin{itemize}
\item RMSE between measured data and the physical model, PM$\in$\{PM1,PM2,PM3\}
\item RMSE between measured data and a static affine model from the PM output: $\hat{y} = ay_p+b$ with coefficients $a$ and $b$. 
\item RMSE between measured data and the CFIR prediction
\item RMSE between measured data and the PMFIR prediction, i.e. CFIR + EFIR.
\end{itemize}
For SNN setups, only the RMSE between measured data and the SNN prediction was computed.

The test series allows to assess the impact of the amount of data, sampling, model order and modeling type. Looking at Figs. \ref{fig:FigRMSE1}--\ref{fig:FigRMSE2} it can be seen that on a short data set the CFIR and PMFIR can closely mimick the training set, on $\nu$ the RMSE on training data is even smaller than the noise component N(0,1). The RMSE increases on test set, however, indicating mild overlearning. This is not observed with 2 month or 1 year data sets, the RMSE of which are in many cases not distinguishable. There is a clear difference between the two extreme scenarios PM1 and PM2, which is visible in particular with $c_B$. A better PM results also to a better PMFIR RMSE. The RMSE are insensitive to differences in nominal case order and AIC-optimal orders. This is significant especially with $c_B$ where the AIC-optimized order was one. Hence, the FIR-window  was not beneficial from the point of view of improving the prediction RMSE. The result is confirmed by the small differences between the affine static modeling and CFIR modeling with both $c_B$  and $\nu$. This observation indicates that the dynamics of the worst case scenarios (used as plant models) are in fact not so far from the true ones. 

A completely different result is depicted by observing the RMSE for the case of static PM (PM3). A steady state PM is a special case of process models. In PM3, the physical model of the van der Vusse CSTR plant was replaced by its approximate steady state description, i.e. 
\begin{equation}
\begin{array}{ll}
\mathbf{x}\left(k+1\right)=f_{PM-ss}\left(\mathbf{u}\left(k\right)\right)
\\
\mathbf{y}_{PM}\left(k\right)=g_{PM}\left(\mathbf{x}\left(k\right)\right)
\end{array}
\end{equation}
where $f_{PM-ss}$ is the steady state solution of the van der Vusse CSTR model. In series with CFIR, it leads to Hammerstein structures, which are well-known and relatively commonly used in the industry due to their robustness. Due to lack of dynamics in the PM, larger filter orders were chosen: $n'=(30,90)$ corresponding to 1 and 3 hour windows (the AIC optimized orders were $n^*=(9,705)$). Observing Figs. \ref{fig:FigRMSE1}--\ref{fig:FigRMSE2}, while static affine mapping of the PM does reduce the RMSE from approximately 0.85 to 0.75 for $c_B$ (4.5 to 3.5 for $\nu$), the addition of FIR dynamics to the static mapping reduces it further to 0.5 (below 2). The overall RMSE's on the PMFIR based on steady state PM still remain rather unsatisfactory, however, as the van der Vusse dynamics are strongly nonlinear.

The static affine mapping can be considered as an alternative to PMFIR, with the advantages of being particularily robust to estimate and simple to implement. Clearly, it is very useful if the plant dynamics are well modelled by the physical model, and only a static adjustment is needed. Another alternative to nonlinear dynamic modeling is a standard feedforward artificial neural network. It is a well known representative of the long list of approaches of data-driven modeling. For comparison, FIR SNN time-series model structures were constructed, where the inputs consisted of the delayed $n=(10,60)$ inputs (normalized flows and heat withdrawn from the coolant), sampled at 2 min. The one-hidden layer structure consisted of $h$ hidden nodes, $h={3,8,15}$ with a bias term, and a linear output node with bias term. The SNN were trained using the default setup of Matlab Deep Learning Toolbox \textit{feedforwardnet} and \textit{train} functions (Levenberg--Marquardt for parameter estimation and reservation of 15 \% of training data for early stopping based on performance on validation set. The best-performing SNN structure (in terms of performance on test set) is reported for each data set. Also alternative FIR SNN structures were experimented, where inputs consisted of smaller subsets of the full window, but the conclusions remained the same. 

In general, the data-driven modelling easily suffers from lack of data when system dimensions increase. Looking at Figs \ref{fig:FigRMSE1}--\ref{fig:FigRMSE2}, the training on short data set seems to provide promising results, but the trained SNN perform poorly on test data. For two day data SNN case, the RMSE on test set is so high it goes outside the window (RMSE 0.26 and 7.7 for $c_B$ and $\nu$, respectively). Figure \ref{fig:FigPred} shows also a prediction given by the SNN. The figure seems to indicate that the static behaviour has been captured somewhat but that the combination of a large dimension of the delayed inputs and sparsity of the data did not lead to meaninful learning outcomes of system dynamics. These results are in line with those obtained using the finite state FIR (FFIR) (\cite{ANZCC}). This could be remedied to some extent by replacing the (nonlinear) FIR by an alternative time series approach. For example, the NOE is known to perform well (\cite{APIC}) but is complicated by the need to compute the impact of model dynamics to gradients, and the requirements it poses on data and system stability.


\begin{figure}
\begin{center}
\includegraphics[width=8.4cm]{FigRMSE1.eps}
\caption{RMSE between measured $c_B$ and predictions by various appoaches.}
\label{fig:FigRMSE1}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=8.4cm]{FigRMSE2.eps}
\caption{RMSE between measured $\nu$ and predictions by various appoaches.}
\label{fig:FigRMSE2}
\end{center}
\end{figure}


\subsection{State estimation}


In the plant data used in the above simulations, the feed temperature $v_0$ was a constant at 104.9 $^{\circ}$C. As there is no information (measured data) of the impact of feed temperature to system outputs, any data-driven approach based on the above data sets would fail in estimating the effect of a change in this temperature. Since this variable is available in the plant physical model  $(f_{PM},g_{PM})$, it is straightforward to extrapolate the impact of $v_0$ in a PM based approach. The extrapolations also carry to the calibrated model. The same applies to all variables/constants in the PM model. This simple case example again motivates taking full advantage of available physical models, as supported by the PMFIR approach.

An unscented Kalman filter (UKF) was constructed to illustrate the state estimation in the van der Vusse CSTR example. The unmeasured feed temperature was changed between $\{99.9,104.9,109.9\}$ $^o$C. An augmented state space was defined by $\mathbf{x}=[c_A,c_B,\nu,\nu_K,\nu_0]'$ and an UKF set to estimate all the states using the PMFIR estimated model for the nominal case. The simulations were started with an incorrect initial mean ($[2, 1, 100, 100, 100]'$) [mol/l and $^o$C] with covariance ($diag[.1, .1,  1, 1, 1]^2$). The covariances of the zero-mean noise components were picked as $Q=\frac{1}{100}diag[0.01, 0.01, 2, 2, 1]^2$ for system state propagation, measurement noise characteristics were assumed to be correctly known, $R=diag[0.01, 1]^2$. 

The estimated feed temperatures are show in Fig. \ref{fig:FigSEUKF}. Three stepwise disturbances impacted to $v_0$ during the simulation (see top plot). The estimation is started at 2 hours from initialization (a sufficient input history is required by PMFIR). All states were correctly estimated by UKF, Fig. \ref{fig:FigSEUKF} top plot shows the estimated feed temperature $v_0$. The innovation signal $(y_ {PM}-q_{PM})$ for both outputs is illustrated in the other plots, together with the PMFIR prediction and measured data.

\begin{figure}
\begin{center}
\includegraphics[width=8.4cm]{FigSEUKF.eps}
\caption{Estimation of feed temperature $\nu_0$ using UKF.}
\label{fig:FigSEUKF}
\end{center}
\end{figure}


\section{Discussion and conclusions}

The application of physical models is widespread in many fields of engineering. These models need to be tuned due to various reasons, such as correction of simplifications/idealisations made during modeling or adaptation to systematic errors in plant measurements. Calibration of physical models is constrained by several engineering requirements, such as emphasizing a clear distinction between the original physical model and its tuning, automaticity from data, and simplicity in implementation. In process control and state estimation, system dynamics are of particular interest. 

The proposed apporach can be seen as a machine learning approach fusing process models with measured data (\cite{Qin}; \cite{Bikmukhametov}). The requirements for such approaches in industrial process engineering were widely discussed in the paper, emphasizing the aspect of physical model tuning and that a model structure founded on a physical model can be expected to extrapolate much better to operation regions not visited in real life, i.e. for which no data is available. These are typically of interest when optimizing for improved plant performance or describing faulty or abnormal situations.

This work proposed an approach for calibrating an existing physical simulation model. In the suggested PMFIR approach, the main characterics of the mapping originate from the physical model. The CFIR provides a linear adjustment for the physical model outputs, the EFIR extends the tuning by a linear mapping from past inputs, thereby enabling identification of linear components not modelled by the physical model. The approach can be automatized with high robustness, and recursive implementations of the algorithms are readily available. 

A typical problem with physical models of industrial processes is that they corrupt over time. Consequently, the models eventually become useless unless updated by a specialized process modeling expert. A priori, identification / machine learning techniques can be readily used for adaptive data-driven modeling. However, it is rare that re-identification of such models could be constructed automatically, but a data scientist is required. If the plant varies (slowy) with time, the linear PMFIR corrections (given some forgetting scheme is used) will make global corrections to the response. It can be expected that such procedure could be automated in a robust fashion.

The suggested PMFIR approach was illustrated and experimented with using a well known process control design benchmark of a highly nonlinear CSTR with a cooling jacket. The incorrect (worst-case) models were successfully calibrated using nominal plant data and even a simple steady-state PM could be succesfully applied. The approach was compared with alternative data-driven methods, including static affine correction at the physical model output, and data-driven modeling of the dynamic process. The proposed approach provided improved performance and was shown to be much less sensitive to the amount of data than corresponding data-driven approaches. Finally, the application of the calibrated model in state estimation using UKF was succesfully demonstrated.

On-going work  focuses at robust estimation of linear dynamics and approaches for merging predictions from several plant models. In applied work, the focus is in gaining experiences on the potential in applications in the heavy process industry. 

%%

%\bibliography{ifacconf}  
\begin{thebibliography}{99}

%\bibitem{Bikmukhametov} T. Bikmukhametov and J. J{\"a}shke, Combining machine learning and process engineering physics towards enhanced accuracy and explainability of data-driven models, Computers and Chemical Engineering, vol. 138, pp. 1-27, 2020.
\bibitem[Bikmukhametov and J{\"a}chke (2020)]{Bikmukhametov} T. Bikmukhametov and J. J{\"a}shke. 
\newblock Combining machine learning and process engineering physics towards enhanced accuracy and explainability of data-driven models.
\newblock Computers and Chemical Engineering, vol. 138, pp. 1--27, 2020.

\bibitem[Chen et al., 1995)]{Chen} H. Chen, A. Gremling and F. Allg{\"o}wer.
\newblock Nonlinear predictive control of a benchmark CSTR. 
\newblock European Control Conference 1995, pp. 3247--3252.

\bibitem[Gordon et al. (1993)]{Gordon} N. Gordon, D. Salmond, and A. Smith.
\newblock Novel appoach to nonlinear/non-Gaussian Bayesian state estimation.
\newblock IEEE Proceedings-F, vol. 140, no 2, pp. 107--113, 1993.

\bibitem[Hastie et al. (2017)]{HastieESLII} T. Hastie, R Thibshirani, and J Friedman.
\newblock The elements of statistical learning - Data mining, inference and prediction. 
\newblock 2nd ed, Springer, 2017.

\bibitem[Ikonen and Najim (2001)]{ProcME} E. Ikonen, and K. Najim.
\newblock Nonlinear process modeling based on a Wiener approach.
\newblock Proceedings of the Institute of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, vol. 215, no 1, 15--27, 2001.

\bibitem[Ikonen and Najim (2002)]{APIC} E. Ikonen, and K. Najim.
\newblock Advanced Process Identification and Control.
\newblock New York: Marcel Dekker, 2002.

\bibitem[Ikonen et al. (2016)]{CCE16} E. Ikonen, I. Selek and K. Najim.
\newblock Process control using finite Markov chains with iterative clustering.
\newblock Computers and Chemical Engineering, vol. 93 pp. 293--308, 2016.

\bibitem[Ikonen and Selek (2020)]{ANZCC} E. Ikonen and I. Selek.
\newblock Calibration of physical models with process data using FIR filtering. 
\newblock Australian and New Zealand Control Conference, 2020.

\bibitem[Julier et al. (1995)]{Julier} S. Julier, J. Uhlmann, and H. Durrant-Whyte.
\newblock A new approach for filtering nonlinear systems.
\newblock American Control Conference, 1620-1632, 1995.

\bibitem[Nikolaou and Vuthandam (1998)]{Nikolaou} M. Nikolaou and  P. Vuthandam.
\newblock FIR model identification: Parsimony through Kernel compression with Wavelets.
\newblock AIChE Journa,l vol. 44 pp. 141-150, 1998.

\bibitem[Qin and Chiang (2019)]{Qin}  S. J. Qin and L. H. Chiang.
\newblock Advances and opportunities in machine learning for process data analytics.
\newblock Computers and Chemical Engineering, vol. 126, 465--473. 

\bibitem[Simon (2006)]{Simon} D. Simon.
\newblock Optimal State Estimation - Kalman, Hinf and nonlinear approaches. 
\newblock John Wiley, 2006.

\bibitem[Sjoberg et al. (1995)]{Sjoberg} J. Sj{\"o}berg, Q. Zhang, L. Ljung, A. Benveniste, B. Dylon, P.-Y. Glorennec, H. Hjalmarsson, and A. Juditsky.
\newblock Nonlinear black-box modeling in system identification: a unified overview.
\newblock Automatica, vol. 31, pp. 1691--1724, 1995.

\bibitem[Skogestad (2009)]{Skogestad} S. Skogestad.
\newblock Chemical and Energy Process Engineering.
\newblock CRC Press, 2009.

\bibitem[Spivey (2010)]{Spivey} B. J. Spivey, J. Hendegren, and T. Edgar.
\newblock  Constrained nonlinear estimation for industrial process fouling. 
\newblock Industrial and Engineering Chemistry Research, vol 49, pp. 7824--7831, 2010.

%\bibitem[Sjoberg (1995)]{Sjoberg} J. Sj{\"o}berg, Q. Zhang, L. Ljung, A. Benveniste, B. Dylon, P.-Y. Glorennec, H. %Hjalmarsson and A. Juditsky, Nonlinear black-box modeling in system identification: a unified overview, Automatica, vol. %31, pp. 1691-1724, 1995.

%\bibitem{APIC} E. Ikonen and K. Najim, Advanced Process Identification and Control, New York: Marcel Dekker, 2002.
% \bibitem{Skogestad} S. Skogestad, Chemical and Energy Process Engineering, CRC Press, 2009.
%\bibitem{ANZCC} E. Ikonen and I. Selek, Calibration of physical models with process data using FIR filtering. Submitted to ANZCC 2020.
%\bibitem{Nikolaou} M. Nikolaou and  P. Vuthandam, FIR model identification: Parsimony through Kernel compression with Wavelets, AIChE Journa,l vol. 44 pp. 141-150, 1998.
%\bibitem{Simon} D. Simon, Optimal State Estimation - Kalman, Hinf and nonlinear approaches. John Wiley, 2006.
%\bibitem{Gordon} N. Gordon, D. Salmond and A. Smith, Novel appoach to nonlinear/non-Gaussian Bayesian state estimation, IEEE Proceedings-F, vol. 140, no 2, pp. 107-113, 1993.
%\bibitem{Julier} S. Julier, J. Uhlmann and H. Durrant-Whyte.A new approach for filtering nonlinear systems, American Control Conference, 1620-1632, 1995.
%\bibitem{Spivey} B. J. Spivey, J. Hendegren and T. Edgar, Constrained nonlinear estimation for industrial process fouling. Industrial and Engineering Chemistry Research, vol 49, pp. 7824-7831, 2010.
%\bibitem{Chen} H. Chen, A. Gremling and F. Allg{\"o}wer, Nonlinear predictive control of a benchmark CSTR. In European Control Conference 1995, pp. 3247-3252.
%\bibitem{CCE16} E. Ikonen, I. Selek and K. Najim, Process control using finite Markov chains with iterative clustering, Computers and Chemical Engineering, vol. 93 pp. 293-308, 2016.
%\bibitem{Ungarala} S. Ungarala, Z. Chen and K. Li, Bayesian state estimation of nonlinear systems using approximate aggregate Markov chains, Ind. Eng. Chem. Res., vol. 45 pp. 4208-4221, 2006.
%\bibitem[Ungarala et al. (2006)]{Ungarala} S. Ungarala, Z. Chen and K. Li.
%\newblock Bayesian state estimation of nonlinear systems using approximate aggregate Markov chains.
%\newblock Ind. Eng. Chem. Res., vol. 45 pp. 4208-4221, 2006.
%%\bibitem{Powell} W. J. Powell. Approximate Dynamic Programming - Solving the curses of dimensionality. Wiley-Interscience, 2010.
%\bibitem[Powell (2010)]{Powell} W. J. Powell. 
%\newblock Approximate Dynamic Programming - Solving the curses of dimensionality. 
%\newblock Wiley-Interscience, 2010.
%\bibitem{Hsu} C. S. Hsu, Cell-to-Cell Mapping - A method of global analysis for nonlinear systems, Springer Verlag: 1987.
%\bibitem{Kemeny} J. Kemeny and J. Snell, Finite Markov Chains. New York: van Nostrand, 1960.
%%\bibitem{Lee} J.H. Lee, W. Wong, Approximate dynamic programming approach for process control, Int. J. Process Control, vol.  20 pp. 1038-1048, 2010.
%\bibitem{Poznyak} A.S. Poznyak, K. Najim and E. Gomez-Ramirez, Self-Learning Control of Finite Markov Chains, New York: Marcel Dekker, 2000.
%\bibitem{Li19} Z. Li, J. Jiang, L. Hong and J.-Q. Sun, On the data-driven generalized cell mapping method, Int. J. Bifurcation and Chaos, vol. 29 pp. 2019.
\end{thebibliography}


    % in the appendices.
\end{document}