Markov Chain Monte Carlo Analysis of the Variable- Volume Exothermic Model for a Continuously Stirred Tank Reactor

In this paper, a variable-volume Continuously Stirred Tank Reactor (CSTR) deterministic exothermic model has been formulated based on the Reynold Transport Theorem. The numerical analysis of the formulated model and the identifiability of its physical parameters are done by using the least squares and the Delayed-Rejection Adaptive Metropolis (DRAM) method. The least square estimates provide the prior information for the DRAM method. The overall numerical results show that the model gives an insight in describing the dynamics of CSTR processes, and 14 parameters of the CSTR are well identified through DRAM convergence diagnostic tests, such as trace, scatter, autocorrelation, histograms, and marginal density plots. Global sensitivity analysis was further performed, by using the partial rank correlation coefficients obtained from the Latin hypercube sampling method, in order to study and quantify the impact of estimated parameters, uncertainties on the model outputs. The results showed that 7 among the 14 estimated model parameters are very sensitive to the model outcomes and so those parameters need to be handled and treated carefully. Keywords-parameter identifiability; variable volume; exothermic; CSTR; RTT; MCMC; DRAM


I. INTRODUCTION
During the past decades, Continuously Stirred Tank Reactors (CSTRs) have gained research momentum as important industrial and chemical production tools. For controlling the reactors, various methods have been proposed to tackle the complexity and the non-linearity operational behaviors that are present in the tank reactor during the production processes [1][2][3][4][5][6]. Discussions about CSTRs seem to be broad and range from general to specific purposes. For example, mathematical modeling and numerical simulations of two-phases which are gas-liquid flow in the CSTR can be found in [7] and the Fokker-Plank Equation was applied for a two-state stochastic CSTR system in [8]. In [9], the robust feedback linearization of an isothermal CSTR was conducted by using the mixed sensitivity synthesis and iteration approaches in the presence of uncertainties. A one state variable, temperature, of a non-isothermal CSTR was analyzed by using Proportional Integral Derivative (PID) and fuzzy logic controllers, and the results from simulations and temperature control show that fuzzy logic can be adopted as a good controller of the process compared with the PID controller [10,17]. The effects of hydrodynamic shear on biogas production in the CSTR using the Metzner-Otto method were analyzed and discussed in [11]. The Bayesian approach was used in [12] as the sorption parameter identifiability tool. The research outputs showed that the Bayesian inference is a more preferable method for the analysis of CSTR experiments as per numerical identifiability as well as the sorption parameter identifiability. The efficient Azo Dye color identification in the CSTR with the built-in bio-electrochemical system was developed for Azo dye alizarin Yellow R (AYR) which in turn can help in wastewater treatment [13]. Authors in [14] investigated the performance of CSTR as bioreactor for producing biohydrogen from water melon waste in the anaerobic digester. The Lyapunov-based stochastic non-linear model predictive control was used in [15] to shape the state probability density functions in the CSTR with the exothermic reaction k A B → , where ݇ is the reaction rate, ‫ܣ‬ is the reactant and ‫ܤ‬ is the product. The Luenberger fuzzy observer with and without sliding modes, the Walcott-Zak fuzzy observer, and the Utkin fuzzy observer were adopted and used as fault detection sensors of the CSTR in [16]. A two state CSTR stochastic model was studied and analyzed in [18] by using the Approximate Expectation Maximization (AEM) and Bayesian algorithms. It was revealed that the Bayesian algorithm is an effective method to be applied on the CSTR's stochastic model since it provided more accurate parameter estimates compared with AEM, and it is even more applicable for an unknown system with a small number of data sets. Authors in [19] used Monod and Haldane kinetics methods to perform the stability analysis of a system that models the formation of biofilms inside the CSTR during the waste-water treatment process. Even though both methods performed well, still the Monod kinetics provided biofilms formation in a shorter time compared to the opponent method.
The effect of operating conditions on the CSTR's performance with a saponification experiment was conducted in [20] and the results showed that the increase in conversion scale may depend on the increase in CSTR's volume. The neural network approach was used in [21] in order to identify the dynamics of two states, namely the temperature and the concentration of the CSTR's model, and reasonable and precise results were obtained. Three-dimensional Computational Fluids Dynamics (CFD) simulations were carried out [22] in order to identify the flow behaviors in the CSTR. Authors in [23] used the cascade control strategy to control the temperature of the exothermic CSTR with a cooling jacket. The stability analysis of the system was investigated with the Routh-Hurwitz and Argand diagram techniques. As a result, the cascade control was pointed out to be an efficient control method for the CSTR processes. A modified CSTR model for the neutralization process was studied and analyzed in [24]. This model has been used to assess the effects of strong acid (HCl) and strong base (NaOH) on the flow rates of ionic concentrations [24]. Firstorder and higher-order sliding mode observer methods were used in [25] in order to design and estimate states and unknown inputs of the CSTR, and it was shown that higher-order sliding modes may be adopted to reduce the noise into the system compared with the first-order sliding mode. Parameter estimation of non-linear chemical and biological processes with non-measured variables from a number of data sets was performed in [26] using the Bayesian approach. Furthermore, a mathematical model and simulation of reactors with the production experiment of hexane from benzene were conducted in [27] and the numerical investigation of phenol oxidation from waste-water inside a reactor was carried out in [28]. Nonparametric and non-linear stochastic dynamical model along with the behavior analysis of a class of the single state isothermal CSTR was also studied and analyzed in [29]. The adaptive method with recursive identification and the polynomial synthesis with placement of poles were applied on the CSTR system to determine its dynamics, however this method provided inappropriate control responses and overshoots [30]. The problem of characterizing the global dynamics of a single state non-linear stochastic CSTR system was addressed in [29] using the Fokker-Plank as the state probability distribution function, but the study of a several-state non-linear stochastic system is paramount as recommended in this article. In the same way, different types of reactors and types of reactions in the chemical engineering processes are widely defined, explained and described in [31]. The modeling and control of the CSTR were also conducted based on a mixed logical dynamical model which resulted in a satisfactory performance of the tank [32]. In [33], the general model of the CSTR was developed and the transient behavior for irreversible non-linear polymerization process in the CSTR was studied and analyzed. Chemical process hazards, causes and proposed measures of safety of batch and semi-batch processes were discussed in [34]. The dynamical behavior of the CSTR through a single first order reaction was researched and analyzed in [35]. The limitations of CSTRs' performance due to cooling jacket dynamics with both open and closed loops are well explained and discussed in [36]. Authors in [37] performed multivariate character and stability analysis of irreversible exothermic CSTRs. The signal flow diagram and the equilibrium states were determined by taking into consideration the first and the second-order reactions.
Often, the volume is treated as a constant [4,25,38,39], however, in reality the tanks may expand and lead to variable volume. Moreover, the parameters are randomly assumed and the estimation of physical parameters that are very influential in the system variation using MCMC methods is lagging behind. It is with these reasons this paper aims to treat the CSTR volume as a variable and identify 14 CSTR physical parameters.

II. MODEL FORMULATION
The assumptions on the parameters and physical properties inside and outside the tank have been considered during model formulation and are adopted from [25,[38][39][40], except for the volume which is considered to be a variable. To depict the problem, the CSTR dynamics are schematically illustrated using the CSTR and the Reynold Transport Theorem (RTT) diagrams.

A. Assumptions
The CSTRs under investigation assume perfect mixing to avoid spatial gradients of velocity, temperature, concentration and of other properties of the mixture.
• The CSTRs assume non-viscous liquid and gas phases which are frequently supplied in it and a static mixer which makes the shaft work produced by the stirring process to be negligible.
• No pressure drop is taking place in the CSTRs, i.e. CSTRs work at a constant pressure.
• Kinetic energy, potential energy, and other forms of external energy are infinitesimally small compared to the heat exchange and the heat from the chemical reactions.
• The wall is isolated and its temperature is negligible. Only the heat exchange is channelled through the designed area.
• The CSTRs assume the volume variation.
• The CSTRs assume density ߩ and specific heat capacitiy ‫ܥ‬ to be constants. • Since there is negligible external stress acting on the system, it is assumed that there is also a negligible momentum on the system.

B. Formulation
Based on schematic diagrams presented in Figures 1 and 2, considering the above listed assumptions, and according to the transport phenomena, mass and energy conservation, the variable-volume exothermic CSTR model formulated is given by the system of equations (1): where Vis the tank volume, F in = F is the feeding velocity, F out is the outlet velocity, C in is the feeding concentration, C is the mixture concentration, E is the activation energy, T is the temperature of the tank, T in is the feeding temperature, T mean is the reference temperature, H * is the reaction enthalpy, ρ is the density, Cp is the heat capacity, U is the heat transfer coefficient, A is the cross-sectional area, T c is the cooling temperature, F c is the cooling velocity, V c is the jacket volume, in c T is the feeding coolant temperature, ρ c is the coolant density, c p C is the jacket's heat capacity, and ݇ ൌ ݇ ݁ Arrhenius equation [41], with ܴ being the gas law constant and ݇ the pre-Arrhenius frequency factor.

A. Markov Chain Monte Carlo Method
The Markov Chain Monte Carlo (MCMC) is a numerical analysis method used as a statistical and Bayesian technique to identify the complex ordinary differential equations' parameters that fit the dynamics of chemical and biological models [42][43][44]. In this paper, the DRAM method is used to estimate 14 parameters of the model (1). The Bayesian inference, qualified to be a very powerful statistical technique, has been widely used to identify the model's parameters ߠ which are obtained after evaluating the parameters' posterior density ‫‬ሺߠ/ ܺ భ , ܺ మ , … , ܺ ே ሻ, where ܺ భ , ܺ మ , … , ܺ ே are the measurement points of the chemical process. For the model (1), Further, the MCMC method is a sampling technique that combines Monte Carlo integration and Markov chains [45]. The overall implementation process starts from proposing a suitable distribution, called proposal distribution, and drawing samples from it [46]. The proposal distribution sometimes depends on the present value to form the chain, which in turn is considered as a Markov chain. The acceptance or rejection mechanisms are employed in the simulation to rectify the trial proposal distribution which ends up with the target distribution. In the end, a simulated chain of parameters (drawn samples, ߠ ଵ , ߠ ଶ , … , ߠ ே ) can be used to approximate the intractable integral (distribution), as: where ‫ܧ‬ሾ݂ሺߠሻ/ ܺ భ , ܺ మ , … , ܺ ே ሻሿis the expectation and ݂ሺߠ ሻ is the density function.  [47,48]. However, the method considered in this paper is the DRAM, presented in Algorithm 1 since it has the desirable feature of tuning the reliable proposal distribution without defining it manually [45,46]. This method overcomes the tedious task of trial and error of tuning the suitable proposal distribution that may appear in the Metropolis-Hastings technique [49]. In this case, the Gaussian distribution is mostly used as a proposal distribution. For a startup, MCMC needs the initial parameter values and a proposal distribution. The DRAM method is used to propose a distribution and its initial parameter values are computed using the classical least square method. The details of the numerical results obtained from the least squares and the MCMC methods are well explained and discussed in Section IV.

B. Least Squares Method
The model in (1) is a system of time dependent differential equations of the form: The least squares method one of the classical optimization methods with the purpose of minimizing the sum of squared residuals. To obtain the residuals, a predictive model is described as: where ܻ represents the observations, ݃൫ܺ , ߠ , ‫ݐ‬ ൯ are the solutions of the model (1) at time ‫ݐ‬ starting from a fixed set of parameters ߠ , and ‫ݎ‬ is the residual term.
The sum of squared residuals function is obtained by taking the sum of squares of ‫ݎ‬ : Hence, the least square method searches the best fitting parameters that minimize the ܵ. ܵ. ܴ taken as the likelihood function. Thus, the fitting set of parameters ߠ can be obtained after solving (2): simultaneously. Almost all chemical processes are intractable and complex due to their non-linearity behaviors. In fact, stepping back to the model (1) with 14 parameters of interest, if we could manually solve (2), we would end up with solving 14 nonlinear equations simultaneously, which is a complicated task. As a simplification, the numerical simulations become a usual way of solving the problem.

C. Delayed Adaptive Markov Chain Monte Carlo
In this paper, the delayed-rejection adaptive metropolis (DRAM) algorithm is used for the parameter identifiability of the deterministic variable-volume exothermic CSTR model described in (1). The DRAM algorithm is chosen since it combines and utilizes both the Delayed-Rejection (DR) and the Adaptive Metropolis (AM) methods to efficiently improve and speed-up the computation process for a slow-start adaptation. The inputs of the algorithm are the initial values of the parameters and all of them are the least square optimal parameter values shown in Table II. The proposed distribution which is a Gaussian with mean 0 and standard deviation Σ taken as the initial covariance matrix, Σ ൌ ଵ షర ூ ೕ ඥ , has been used, where j is the length of model's parameters to be estimated, and ‫ܫ‬ ൌ ‫ܫ‬ × which is an ݆ × ݆ identity matrix.

1) Algorithm 1: DRAM
Draw the initial point ߠ , from initial distribution ‫‬ ሺߠሻ. Set an initial non-adaptive period ܰ and initial covariance matrix Σ .
For ݅ ൌ 1, 2, … perform the following: (iii) Acceptance/rejection rule by setting: (iv) If ݊ ≪ ݅ (or after every ݊ ௧ iterations), then update the covariance matrix by using: where ‫ܫ‬ ௗ is an ݀ × ݀ identity matrix and ߝ is a small positive number that makes the matrix Σ to be non-singular [46,52,53].

IV. NUMERICAL ANALYSIS OF THE MODEL AND DISCUSSION
The formulated model (1) is solved numerically by the 4 ௧ order Runge-Kutta method which is the ode45 solver available in Matlab R2016b software package. The parameter identifiability is done using the least squares and MCMC.

A. Numerical Simulations
Due to lack of actual information (real data) about the system, the model (1) is simulated by using literature and assumed values (Table I), the discrete sampling time points are 100, and 100×4 numerical solutions of the model ܻ presented in (1) are obtained. The numerical results from the subplot 1 of Figure 3 reveal that the volume of tank reactor increases from 100m 3 to approximately 126m 3 and this is an indication of having non-constant flow rates of reactants due change at both the inlets and outlets. Numerical solutions of the model. Figure 3 also shows that the reactants are consumed inside the reacting tank as its concentration approaches zero, which means the complete mixing and at the same time symbolizes a non-partial conversion of reactants into products that may lead to time residence distribution analysis as one of the inconveniences of CSTRs. Along this process, there is a covering cooling jacket that communicates with the reacting tank through a designed cross-sectional area of 0.015m 2 to cool down the rising temperature inside the reacting tank. The covering cooling jacket contains the substance whose temperature is initially lower than the starting temperature of the reacting tank to disable the explosion of the reaction. During this process, the temperature of the reacting tank rises from 298.35K to its operating working temperature that is373.48K. Meanwhile, the cooling temperature of the covering cooling jacket also rises from its initial temperature 288.15K until it reaches 363.14K. As a result, if this scenario is selected to be a piloting tank, then all the simulated information and the working conditions that are described above have to be taken into consideration quantitatively.

B. Least Square Results
The deterministic variable-volume exothermic CSTR model presented in (1) has been numerically analyzed by using the least square method and 14 parameters of the model were optimized from the noisy measurements of the system. The noisy measurements are obtained after corrupting the obtained empirical data with Gaussian noise of 0.05 standard deviation. The obtained results are presented in Figure 4, and the 14 estimated parameters and their corresponding initial values are also presented in Table II. According to the results, all parameters relatively converge to their corresponding initial values and the measurements fit the exothermic CSTR model as can be viewed from fitted state variables shown in Figure 4. Least square fitting model. We can see that the measurements fit the model well as parameters vary.

C. Markov Chain Monte Carlo Results
The initial values for the MCMC are estimated using the least squares method. Since we have fourteen parameters to be estimated, a sample consisting of 100,000 ×14 parameters has been generated during simulations and the method has been adapted 100 times. Also the initialized covariance matrix of MCMC is chosen to be Σ ൌ ଵ షర ூ భరൈభర √ଵସ . Details of the MCMC results are discussed by the following MCMC diagnostic tests.

1) MCMC Diagnostics
The identifiability of the model parameters is mainly based on the convergence of the MCMC. There are several statistical and graphical convergence tests for the MCMC methods [54][55][56][57]. In this paper, we use the common statistical inference and graphical analyses of the trace (time series plots), scatter plots, autocorrelation plots, histograms, and the marginal density distributions for each drawn sample of the parameters to identify the model and diagnose the convergence of the generated MCMC samples.

2) Trace Plots
One way to analyze the convergence of the MCMC is to check the mixing of the generated sample of posteriors through trace plots. If the generated chain of posteriors becomes stationary for several initial values, and there are no obvious spikes, it is an indication of having a good mixing which is a good sign of convergence. In Figure 5 we can see that the mixing of samples is very good so the chain converges, and the 14 sampled values of posteriors are the means (centers) of the samples.

3) Scatter Plots (pairs)
A poor convergence of MCMC method inventively leads to high correlation between estimated parameters. Since there are 14 parameters to be identified, then there are 91 scatter plots which we need to investigate if there are strong correlations between them. Figure 6 shows the scatter plots for the first 10 sampled parameters equivalent to 45 scatter plots and how these posterior samples correlate to each other. We observe that there is no high correlation between the pairs of the estimated parameters and so the parameters of the model (1) are adequately identified. Figure 7 determines and examines the correlation between consecutive samples during posteriors chain sampling. We can see that the coefficients of autocorrelation functions of all generated samples (x-axis) tend to zero as the number of lags (y-axis) increases and get stationary around zero after 100 lags. That is an indication of having a good mix.

5) Marginal Distribution (Density) of the Sampled Parameters
Another diagnostic test is to plot the posteriors density distributions. Normally for better mixing we expect the distributions of all density estimations to follow a Gaussian distribution.  Autocorrelation plots of the sampled parameters.

6) Histograms of the Sampled Parameters
From Figure 9, we can see that almost all parameter values follow the normal distribution. The MCMC's convergence tends to be bad for ‫ܨ‬ and ܸ . However, the other 12 parameters are fine. This could be caused by different and uncertain reasons. It may stem from either the complexity of the model, the little information about their mean values, or the method itself has failed to sample well these two parameters from the proposed distribution.
The results in Table II show that almost all parameters converge to their initial values and posterior means are within their credible intervals. We can see that the posterior means and posterior medians are nearly equal. We can in addition observe that the Markov chain errors and the standard deviation of the model parameters are significantly small which shows that the method performed well in identifying the best estimates. The computed autocorrelation time (tau) values are small with the significance of having independence in sampling the successive posteriors as well as the convergence of the method. All Geweke values are almost 1, except 0.21093 for ‫ܨ‬ and 0.41264 for ܸ with the indication of having non-stationary means values by default for the start (10%) and the end (50%) of sampling of these two parameters. In reality, the skewness and kurtosis values for the Gaussian distribution are approximately expected to be 0 and 3 respectively. So, the skewness and kurtosis values for the obtained samples show that only ‫ܨ‬ and ܸ values are somehow skewed in the right hand side and do not incorporate the exact features of the Gaussian distribution.

D. Global Sensitivity Analysis
To quantify the uncertainty effects of estimated parameters on the model outputs, we performed global sensitivity analysis by using Latin hypercube sampling method. This method computes the Partial Rank Correlation Coefficients (PRCCs) and determines which parameters are globally very sensitive to the model. PRCCs values vary between -1 and +1 with high sensitivity index for values approaching ±1 and low sensitivity index for values which are far from ± 1. The PRCCs values are graphically displayed in Figure 10 and their results can be seen in Table III.
From Figure 10, we observe that the 1 st parameter significantly and negatively affects the model volume outputs. The increase in that parameter decreases the volume. From subplot 2, the 3 rd parameter, which is ݇ , has high negative effect on the concentration outputs of the model. The increase  Table  III.

V. CONCLUSION
This paper focused on the mathematical formulation of the deterministic variable-volume exothermic CSTR model, and its numerical simulations were tested to supplement the theoretical results as applications using literature values. The identifiability of the parameters required real system measurements, but due to the absence of real data, the system measurements were simulated, i.e. the numerical solutions were empirical data corrupted with Gaussian noise with a standard deviation of 0.05.
The identifiability of the physical parameters of the formulated model was numerically carried out by using the least squares and DRAM. The least square estimates converged to the literature values and were treated as prior information for the DRAM method. The generated DRAM samples were graphically and statistically analyzed to test the convergence. The results show that the model was solved and its physical parameters were well identified. To study and quantify the associated parameter uncertainties to the model outcomes, global sensitivity analysis was performed by using the Latin hypercube sampling method. Seven parameters among the 14 estimated model parameters were found to be very sensitive to the model outputs and therefore, these parameters need to be controlled effectively.