Microsoft Word - ETASR_V11_N2_pp6919-6929 Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6919 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … Markov Chain Monte Carlo Analysis of the Variable- Volume Exothermic Model for a Continuously Stirred Tank Reactor Jean Pierre Muhirwa Department of Applied Mathematics and Computational Science (AMCS), Nelson Mandela African Institution of Science and Technology Tanzania (NM-AIST) Arusha, Tanzania and Department of Mathematics, University of Rwanda, College of Science and Technology (UR-CST), Kigali, Rwanda muhirwaj@nm-aist.ac.tz Isambi Sailon Mbalawata Research Department African Institute for Mathematical Sciences (AIMS) Kigali, Rwanda imbalawata@nexteinstein.org Verdiana Grace Masanja Department of Applied Mathematics and Computational Science (AMCS) Nelson Mandela African Institution of Science and Technology Tanzania (NM-AIST) Arusha, Tanzania verdiana.masanja@nm-aist.ac.tz Abstract-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-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 waste- water 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 Corresponding author: Jean Pierre Muhirwa Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6920 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … 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]. First- order 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]. Non- parametric 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-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. Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6921 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … • 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): 1 1 ( ) 0 1 1 ( ) * 0 , ( ) , ( ) ( ) , ( ) ( ) , mean mean in c in out E R T T in E R T T c in p p c c c c c c c c p dV F F dt dC F C C k e C dt V dT F H k e C UA T T T T dt V C V C dT F UA T T T T dt V V C ρ ρ ρ − − − −  = −   = − −   − − = − + −   − = − +  (1) where Vis the tank volume, Fin = F is the feeding velocity, Fout is the outlet velocity, Cin is the feeding concentration, C is the mixture concentration, E is the activation energy, T is the temperature of the tank, Tin is the feeding temperature, Tmean 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, Tc is the cooling temperature, Fc is the cooling velocity, Vc is the jacket volume, inc T is the feeding coolant temperature, ρc is the coolant density, cp C is the jacket's heat capacity, and � � � � � � �� � �� � ����� � is the Arrhenius equation [41], with � being the gas law constant and � the pre-Arrhenius frequency factor. Fig. 1. Schematic diagram that describes the dynamics of CSTR. Fig. 2. Schematic illustration of the RTT. III. MATERIALS AND METHODS 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- 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), � � #$%&',$, � ,(,)*+,-,. ∗ ,�,��,0,�,$1, 21,�1,��34 and �� � �2, �,),)5�. 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, �6,�7,…,�" ) can be used to approximate the intractable integral (distribution), as: (#8���/ ���,�� ,…, �"�4 9 6 " ∑ 8��; " ;<6 � where (#8���/ ���,�� ,…,�"�4is the expectation and 8��;� is the density function. The Metropolis algorithm, Metropolis Hastings, Hamiltonian Monte Carlo, Gibbs Sampler, Reversible Jump Monte Carlo, Metropolis Adjusted Langevin, Slice Sampling, Multiple Try Metropolis, and Delayed Rejection Metropolis are among the most used MCMC algorithms [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: =>? =' @��, �A,B;C � 0 where � is the number of state variable and �� is the vector of Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6922 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … state variables. For the model (1), � � 4, B; is the discrete time sampling, and �� � �2, �, ), )5 ), �A is the set of parameters of the model to be identified, where F � 1,2, . . . ,14. So, in this paper, the 14 parameters to be identified are: �A = [$%&' , $ , � , (, )*+,- , .∗ , �, ��, 0, �, $1 , 21, �1 , ��3 ] 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: J�; = K@��;, �A , B; C + M�; where J�; represents the observations, K@��; , �A , B; C are the solutions of the model (1) at time B; starting from a fixed set of parameters � , and M�; is the residual term. The sum of squared residuals function is obtained by taking the sum of squares of M�; : N.N. � � ∑ �M�;� 7 � ∑ �J�; − * ;<6 * ;<6 K@��;, �A, B; C)7 Hence, the least square method searches the best fitting parameters that minimize the N. N. � taken as the likelihood function. Thus, the fitting set of parameters �A can be obtained after solving (2): P�Q.Q.R� PST � 0 ≡ P(∑ (V?W�WX� �Y@>?W ,ST,'W C) PST = 0 (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, Σ � 6 [\]T ^A , has been used, where j is the length of model’s parameters to be estimated, and _A � _A×A which is an F × F identity matrix. 1) Algorithm 1: DRAM Draw the initial point � , from initial distribution � ���. Set an initial non-adaptive period a and initial covariance matrix Σ . For b � 1, 2, … perform the following: (i) Sample a current point �c from the proposal distribution d(�c/�;�6 ) (ii) Compute the acceptance probability using: e@�;�6 , �cC = min �1, �(� c/��� , �� , … , ��" )d(�;�6 / �c� ���;�6/��� , �� , … , ��" )d(�c/�;�6 )) (iii) Acceptance/rejection rule by setting: �; = h �c, b8 i < e@�;�6, �cC kℎ�M� i~n�0, 1) �;�6, oBℎ�Mkbp�. (iv) If q ≪ b (or after every q 's iterations), then update the covariance matrix by using: Σ; � �ot�� , �6, … , �; ) + u_= where _= is an v × v identity matrix and u is a small positive number that makes the matrix Σ; to be non-singular [46, 52, 53]. (v) i ← i+1 TABLE I. VARIABLES, CONSTANTS, AND PARAMETERS VALUES Parameters and variables Symbol Unit Physical meaning Values Source �;- wxyz x{| x} Feeding concentration 316.8 [50] � wxyz x{| x} Initial concentration 316.8 [50] � wxyz x{| x} Mixture concentration State variable Simulated );- ~ Feeding temperature 298.35 [50] ) ~ Initial temperature 298.35 [50] ) ~ Mixture temperature State variable Simulated .∗ w��z wxyz Enthalpy -1004.3×10 3 [50] )1� K Initial cool temperature 288.15 [50] )1W� ~ Feeding cool temperature 293 [38] )1 ~ Jacket temperature State variable Simulated � w� ~ wxyz Gas constant 8.314 [50] 2 x } Initial tank volume 100 [51] 2 x} Tank volume State variable Simulated ( w� wxyz Activation energy 0.5 [50] )*+,- ~ Reference temp 298.15 [50] �� w��z ~ w� Heat capacity 4186 [50] 0 w� ~ x{| x� Heat transfer coefficient 100000 Assumed $1 x} x{| Cooling velocity 46.5×10 -6 [50] 21 x} Jacket volume 50×10 -6 [50] �1 w� x} Density of the coolant 1000 [51] ��3 w��z ~ w� Heat capacity of jacket 4186 [50] $%&' m� min Outlet velocity 130×10 -6 [50] $ m� min Feeding velocity 130×10 -2 Assumed � min -1 Pre-Arrhenius factor 0.9 [50] � kg m� Density 1000 [51] � m7 Area 0.015 [50] Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6923 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … IV. NUMERICAL ANALYSIS OF THE MODEL AND DISCUSSION The formulated model (1) is solved numerically by the 4's 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 J�; 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. Fig. 3. 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. Fig. 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 Σ � 6 [\]�\`�\ √6� . 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- 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. Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6924 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … TABLE II. THE STATISTICAL INFERENCE ON THE ESTIMATED PARAMETERS BY USING THE LEAST SQUARE AND MCMC METHODS Param LSQ and MCMC estimates Initial values LSQ estim MCMC posterior mean MCMC posterior median Cred interv Std MCerr tau geweke Kurtosis Skewness $%&' 130×10 -6 130.016×10 -6 130.015×10 -6 130.016×10 -6 [130.015×10 -6 , 130.016×10 -6 ] 7.32 ×10 -8 1.0633×10 -9 41.003 0.99994 2.9626 -0.0630 $ 130×10-2 130.09×10-2 129.971139×10-2 129.9704569×10-2 [129.9705396×10 -2 , 129.9716882×10 -2 ] 0.00093 3.9735×10 -5 62.713 0.99997 3.0915 -0.0225 � 0.9 0.900553 0.900592376 0.900596750 [0.900586482, 0.900598270] 0.000951 2.2856×10 -5 50.329 0.99994 3.0930 -0.0016 ( 0.5 0.493183 0.493486288 0.493470050 [0.493480384, 0.493492191] 0.0009525 1.7754×10 -5 46.313 0.9985 3.1854 0.0934 )*+,- 298.15 298.147 298.151888471 298.151879315 [298.151882554, 298.151894388] 0.00095466 1.9926×10 -5 58.812 1 3.0988 -0.0220 .∗ -1004.3×103 1026.03×103 -1028032.32 9702401 -1028032.329 677610 [-1028032.32970, -1028032.3296 966452] 0.00095973 3.2217×10 -5 57.28 1 3.1340 0.0125 � 1000 1017.16 1026.121280804 1026.121304498 [1026.121274963, 1026.121286644] 0.00094236 2.9081×10 -5 59.644 1 3.0811 -0.0299 �� 4186 4231.14 4217.368624072 4217.368601824 [4217.368617960, 4217.368630185] 0.00098619 2.1133×10 -5 48 1 3.3398 0.0578 0 100,000 100681 100098.50255273 8 100098.502551528 [100098.502547056, 100098.502558419] 0.00091663 4.2388×10 -5 69.113 1 3.2362 0.0790 � 0.015 0.0151979 0.014959065 0.014962113 [0.014953136, 0.014964995] 0.00095668 3.0706×10 -5 53.978 0.99458 3.0965 -0.0279 $1 46.5×10 -6 46.7705×10 -6 74.7936×10 -5 62.4641×10 -5 [74.4298×10-5, 75.1574×10-5] 0.00058698 3.7635×10-5 85.568 0.21093 3.9661 1.0355 21 50×10 -6 49.6872×10 -6 76.2673×10 -5 63.7067×10 -5 [75.8967×10-5, 76.6379××10-5] 0.00059791 3.2879×10-5 69.882 0.41264 3.8227 1.0014 �1 1000 987.751 994.859555082 994.859566589 [994.859549347, 994.859560817] 0.00092534 2.464×10 -5 63.819 1 3.2461 0.0145 ��3 4186 4134.6 4163.918586616 4163.918608754 [4163.918580689, 4163.918592542] 0.00095621 1.7679×10 -5 45.604 1 3.0466 -0.0150 Note: Param=Parameters, LSQ=least squares, estim=estimates, Cred interv=Credible interval, Std=Standard deviation, MCerr=MCMC error Fig. 5. The trace plots 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. 4) Autocorrelation Plots 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. Figure 8 depicts how density distributions of all estimated parameters follow Gaussian distribution and their values are taken as means of distributions. Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6925 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … Fig. 6. The scatter plots of pairs of sampled parameters. We can see that none of the parameters #$%&', $, � ,(, )*+,-,. ∗ , �,��, 0, �4 is correlated to the others because there are no correlation patterns observed. Fig. 7. Autocorrelation plots of the sampled parameters. Fig. 8. The marginal density distribution of the sampled parameters. Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6926 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … Fig. 9. Histograms of the 14 sampled parameters. Fig. 10. Parameter PRCCs' plot for the model outputs. 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 $1 and 21. 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 $1 and 0.41264 for 21 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 $1 and 21 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 Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6927 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … of the temperature profiles positively depends on the increase in the 6 th , 7 th and 8 th parameter values whereas the increase of the 8 th and 10 th parameter values decreases and increases the cooling temperature profiles respectively. So, in total there are 7 identified parameters that may have a great impact on the model response and their values are represented by (*) in Table III. TABLE III. PRCC RESULTS Parameter Partial rank correlation coefficient values for each variable � � � �� $%&' −0.9847 �∗) 0.2545 −0.0479 −0.0429 $ −0.0269 0.4127 0.1482 −0.0119 � −0.1060 −0.7418 �∗) 0.0302 −0.0761 ( 0.0070 0.3054 −0.0743 −0.0973 )*+,- −0.1713 −0.0515 −0.2141 0.0032 .∗ 0.0813 −0.1444 0.6114�∗) −0.3979 � 0.0819 −0.3102 0.5000(∗) −0.2924 �� −0.0345 −0.2137 0.6391(∗) −0.5186(∗) 0 0.0232 −0.0547 −0.1470 0.4348 � −0.0091 0.0298 −0.2086 0.5232(∗) $1 0.1761 −0.2802 0.0415 −0.3810 21 0.0038 −0.0401 −0.1139 0.1187 �1 0.0928 −0.0246 0.1285 −0.2502 ��3 −0.0655 −0.2189 0.0294 −0.3894 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. ACKNOWLEDGEMENT The authors acknowledge the financial support received from the German Academic Exchange Service (DAAD) and the research department at African Institute for Mathematical Sciences (AIMS)-Rwanda Center. REFERENCES [1] A. S. Al-Araji, "Modeling of Continuous Stirred Tank Reactor based on Artificial Neural Network," Al-Nahrain Journal for Engineering Sciences, vol. 18, no. 2, pp. 202–207, 2015. [2] A. O. Ahmed, G. A. Gasmelseed, A. B. Karama, and A. E. Musa, "Cascade Control of a Continuous Stirred Tank Reactor (CSTR)," Journal of Applied and Industrial Sciences, vol. 1, no. 4, pp. 16–23, 2013. [3] H. A. Maddah, "Numerical Analysis for the Oxidation of Phenol with TiO2 in Wastewater Photocatalytic Reactors," Engineering, Technology & Applied Science Research, vol. 8, no. 5, pp. 3463–3469, Oct. 2018, https://doi.org/10.48084/etasr.2304. [4] A. Simorgh, A. Razminia, and V. I. Shiryaev, "System identification and control design of a nonlinear continuously stirred tank reactor," Mathematics and Computers in Simulation, vol. 173, pp. 16–31, Jul. 2020, https://doi.org/10.1016/j.matcom.2020.01.010. [5] A. Uppal, W. H. Ray, and A. B. Poore, "On the dynamic behavior of continuous stirred tank reactors," Chemical Engineering Science, vol. 29, no. 4, pp. 967–985, Apr. 1974, https://doi.org/10.1016/0009- 2509(74)80089-8. [6] A. S. Ibrehem, "Modified Mathematical Model For Neutralization System In Stirred Tank Reactor," Bulletin of Chemical Reaction Engineering & Catalysis, vol. 6, no. 1, pp. 47–52, May 2011, https://doi.org/10.9767/bcrec.6.1.825.47-52. [7] A. Z. Al-Khazaal, F. Ahmad, and N. Ahmad, "Study on the Removal of Thiosulfate from Wastewater by Catalytic Oxidation," Engineering, Technology & Applied Science Research, vol. 9, no. 2, pp. 4053–4056, Apr. 2019, https://doi.org/10.48084/etasr.2553. [8] B. G. Osorio, H. B. Castro, and J. D. S. Torres, "State and unknown input estimation in a CSTR using higher-order sliding mode observer," in IEEE IX Latin American Robotics Symposium and IEEE Colombian Conference on Automatic Control, Bogota, Colombia, Oct. 2011, pp. 1– 5, https://doi.org/10.1109/LARC.2011.6086829. [9] C. G. Hill, An Introduction to Chemical Engineering Kinetics and Reactor Design. New York, USA: Wiley, 1977. [10] J. Du, C. Song, and P. Li, "Modeling and Control of a Continuous Stirred Tank Reactor Based on a Mixed Logical Dynamical Model," Chinese Journal of Chemical Engineering, vol. 15, no. 4, pp. 533–538, Aug. 2007, https://doi.org/10.1016/S1004-9541(07)60120-7. [11] D. P. Karadimou, P. A. Papadopoulos, and N. C. Markatos, "Mathematical modelling and numerical simulation of two-phase gas- liquid flows in stirred-tank reactors," Journal of King Saud University - Science, vol. 31, no. 1, pp. 33–41, Jan. 2019, https://doi.org/10.1016/ j.jksus.2017.05.015. [12] D. Ndanguza, J. P. Muhirwa, and A. Uwimana, "Modeling and parameters estimation of a Spatial Predator-Prey distribution," Rwanda Journal of Engineering, Science, Technology and Environment, vol. 2, no. 1, pp. 1–17, Jul. 2019, https://doi.org/10.4314/rjeste.v2i1.5. [13] D. Tamboli and R. Chile, "Multi-model approach for 2-DOF control of nonlinear CSTR process," International Journal of Modelling, Identification and Control, vol. 30, no. 2, pp. 143–161, Jan. 2018, https://doi.org/10.1504/IJMIC.2018.094208. [14] E. Shakeri, G. Latif-Shabgahi, and A. E. Abharian, "Design of an intelligent stochastic model predictive controller for a continuous stirred tank reactor through a Fokker-Planck observer," Transactions of the Institute of Measurement and Control, vol. 40, no. 10, pp. 3010–3022, Jun. 2018, https://doi.org/10.1177/0142331217712583. [15] E. Vlahakis and G. Halikias, "Temperature and concentration control of exothermic chemical processes in continuous stirred tank reactors," Transactions of the Institute of Measurement and Control, vol. 41, no. 15, pp. 4274–4284, Nov. 2019, https://doi.org/10.1177/0142331219 855591. [16] E. A. Buehler, J. A. Paulson, and A. Mesbah, "Lyapunov-based stochastic nonlinear model predictive control: Shaping the state probability distribution functions," in American Control Conference, Boston, MA, USA, Jul. 2016, pp. 5389–5394, https://doi.org/10.1109/ ACC.2016.7526514. Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6928 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … [17] E. H. Karimi and K. B. McAuley, "A Bayesian Method for Estimating Parameters in Stochastic Differential," IFAC-PapersOnLine, vol. 48, no. 8, pp. 147–152, Jan. 2015, https://doi.org/10.1016/j.ifacol.2015.08.172. [18] F. Remo, L. S. Luboobi, I. S. Mabalawata, and B. K. Nannyonga, "A mathematical model for the dynamics and MCMC analysis of tomato bacterial wilt disease," International Journal of Biomathematics, vol. 11, no. 1, Aug. 2017, Art. no. 1850001, https://doi.org/10.1142/ S1793524518500018. [19] G. I. Valderrama-Bahamondez and H. Frohlich, "MCMC Techniques for Parameter Estimation of ODE Based Models in Systems Biology," Frontiers in Applied Mathematics and Statistics, vol. 5, 2019, Art. no. 55, https://doi.org/10.3389/fams.2019.00055. [20] G. L. Foutch and A. H. Johannes, "Reactors in Process Engineering," in Encyclopedia of Physical Science and Technology, 3rd ed., R. A. Meyers, Ed. New York, NY, USA: Academic Press, 2003, pp. 23–43. [21] H. Ballesteros-Moncada, E. J. Herrera-Lopez, and J. Anzurez-Marin, "Fuzzy model-based observers for fault detection in CSTR," ISA Transactions, vol. 59, pp. 325–333, Nov. 2015, https://doi.org/10.1016/ j.isatra.2015.10.006. [22] H. Haario, E. Saksman, and J. Tamminen, "Adaptive proposal distribution for random walk Metropolis algorithm," Computational Statistics, vol. 14, no. 3, pp. 375–395, Sep. 1999, https://doi.org/ 10.1007/s001800050022. [23] H. Haario, E. Saksman, and J. Tamminen, "An adaptive Metropolis algorithm," Bernoulli, vol. 7, no. 2, pp. 223–242, Apr. 2001, https://doi.org/bj/1080222083. [24] I. A. Stepanov, "Exact Calculation of the Internal Energy of the Ideal Gas in Statistical Mechanics," Physical Science International Journal, vol. 14, no. 5, pp. 1–5, Apr. 2017, https://doi.org/10.9734/PSIJ/ 2017/32822. [25] I. S. Mbalawata, S. Sarkka, and H. Haario, "Parameter estimation in stochastic differential equations with Markov chain Monte Carlo and non-linear Kalman filtering," Computational Statistics, vol. 28, no. 3, pp. 1195–1223, Jun. 2013, https://doi.org/10.1007/s00180-012-0352-y. [26] I. S. Mbalawata, S. Sarkka, M. Vihola, and H. Haario, "Adaptive Metropolis algorithm using variational Bayesian adaptive Kalman filter," Computational Statistics & Data Analysis, vol. 83, pp. 101–115, Mar. 2015, https://doi.org/10.1016/j.csda.2014.10.006. [27] J. C. Etchells, "Process Intensification: Safety Pros and Cons," Process Safety and Environmental Protection, vol. 83, no. 2, pp. 85–89, Mar. 2005, https://doi.org/10.1205/psep.04241. [28] J. Jiang, J. Wu, S. Poncin, and H. Z. Li, "Effect of hydrodynamic shear on biogas production and granule characteristics in a continuous stirred tank reactor," Process Biochemistry, vol. 51, no. 3, pp. 345–351, Mar. 2016, https://doi.org/10.1016/j.procbio.2015.12.014. [29] J. Vojtesek and P. Dostal, "Simulation analysis of continuous stirred tank reactor," in 22nd European Conference on Modelling and Simulation, Nicosia, Cyprus, Jun. 2008, pp. 1–6. [30] J. Vojtesek and P. Dostal, "Simulation of adaptive control of continuous stirred tank reactor," International Journal of Simulation Modelling, vol. 8, no. 3, pp. 133–144, 2009. [31] J. P. Muhirwa and D. Ndanguza, "Effect of random noise, quasi random noise and systematic random noise on unknown continuous stirred tank reactor (cstr)," Applied Mathematical Sciences, vol. 11, no. 62, pp. 3051–3071, 2017. [32] K. Lopez Buritica, S. Casanova Trujillo, C. D. Acosta, and H. A. Granada Diaz, "Dynamical Analysis of a Continuous Stirred-Tank Reactor with the Formation of Biofilms for Wastewater Treatment," Mathematical Problems in Engineering, vol. 2015, Jun. 2015, Art. no. e512404, https://doi.org/10.1155/2015/512404. [33] K. Cahyari, Sarto, S. Syamsiah, and A. Prasetya, "Performance of continuous stirred tank reactor (CSTR) on fermentative biohydrogen production from melon waste," IOP Conference Series: Materials Science and Engineering, vol. 162, no. 1, Nov. 2016, Art. no. 012013, https://doi.org/10.1088/1757-899X/162/1/012013. [34] K. J. Keesman et al., "Aquaponics Systems Modelling," in Aquaponics Food Production Systems: Combined Aquaculture and Hydroponic Production Technologies for the Future, S. Goddek, A. Joyce, B. Kotzen, and G. M. Burnell, Eds. Cambridge, England: Springer International Publishing, 2019, pp. 267–299. [35] L. P. Russo and B. W. Bequette, "Cstr performance limitations due to cooling jacket dynamics," in Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, J. G. Balchen, Ed. Oxford, UK: Pergamon, 1993, pp. 149–154. [36] M. Danish, M. K. Al Mesfer, and M. M. Rashid, "Effect of operating conditions on cstr performance: an experimental study," International journal of engineering research and applications, vol. 5, no. 2, pp. 74– 78, 2015. [37] M. Laine, Adaptive MCMC methods with applications in environmental and geophysical models. Helsinki, Finland: Finnish Meteorological Institute, 2008. [38] D. F. A. M. N. Esmaeel, "Fuzzy logic Control of Continuous Stirred Tank Reactor," Tikrit Journal of Engineering Sciences, vol. 20, no. 2, pp. 70–80, 2013. [39] M. A. S. Aboelela and R. H. M. Hennas, "Development of a fractional order pid controller using adaptive weighted pso and genetic algorithms with applications," in Fractional Order Systems, A. T. Azar, A. G. Radwan, and S. Vaidyanathan, Eds. Cambridge, Massachusetts, MA, USA: Academic Press, 2018, pp. 511–551. [40] M.-H. Cui, D. Cui, L. Gao, H.-Y. Cheng, and A.-J. Wang, "Efficient azo dye decolorization in a continuous stirred tank reactor (CSTR) with built-in bioelectrochemical system," Bioresource Technology, vol. 218, pp. 1307–1311, Oct. 2016, https://doi.org/10.1016/j.biortech.2016. 07.135. [41] N. M. Ramli and M. S. Mohamad, "Modelling for Temperature Non- Isothermal Continuous Stirred Tank Reactor Using Fuzzy Logic," International Journal of Environmental and Ecological Engineering, vol. 11, no. 2, pp. 169–175, Jan. 2017. [42] R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport Phenomena, 2nd edition. New York, NY, USA: Wiley, 2009. [43] R. Kandiyoti, Fundamentals of Reaction Engineering. London, UK: Bookboon, 2009. [44] R. C. S. Dias and M. R. P. F. N. Costa, "Transient Behavior and Gelation of Free Radical Polymerizations in Continuous Stirred Tank Reactors," Macromolecular Theory and Simulations, vol. 14, no. 4, pp. 243–255, 2005, https://doi.org/10.1002/mats.200400086. [45] R. M. Sudhanan and D. P. Poongodi, "Comparative Analysis of Various Control Strategies for a Nonlinear CSTR System," International Journal of Nonlinear Sciences and Numerical Simulation, vol. 18, no. 3–4, pp. 269–276, Jun. 2017, https://doi.org/10.1515/ijnsns-2015-0125. [46] S. Brooks, "Markov chain Monte Carlo method and its application," Journal of the Royal Statistical Society: Series D (The Statistician), vol. 47, no. 1, pp. 69–100, 1998, https://doi.org/10.1111/1467-9884.00117. [47] S. S. Jang and R. B. Gopaluni, "Parameter estimation in nonlinear chemical and biological processes with unmeasured variables from small data sets," Chemical Engineering Science, vol. 66, no. 12, pp. 2774– 2787, Jun. 2011, https://doi.org/10.1016/j.ces.2011.03.029. [48] S. Masoumi, T. A. Duever, and P. M. Reilly, "Sequential Markov Chain Monte Carlo (MCMC) model discrimination," The Canadian Journal of Chemical Engineering, vol. 91, no. 5, pp. 862–869, 2013, https://doi.org/10.1002/cjce.21711. [49] S. Nanda, "Reactors and Fundamentals of Reactors Design for Chemical Reaction," Ph.D. dissertation, Maharshi Dayanand University, Haryana, India, 2008. [50] S. Sharma, "Markov Chain Monte Carlo Methods for Bayesian Data Analysis in Astronomy," Annual Review of Astronomy and Astrophysics, vol. 55, pp. 213–259, Aug. 2017, https://doi.org/10.1146/annurev-astro- 082214-122339. [51] S. Sinharay, "Assessing Convergence of the Markov Chain Monte Carlo Algorithms: A Review," ETS Research Report Series, vol. 2003, no. 1, pp. i–52, 2003, https://doi.org/10.1002/j.2333-8504.2003.tb01899.x. [52] S. Tronci, M. Grosso, J. Alvarez, and R. Baratti, "Stochastic dynamical nonlinear behavior analysis of a class of single-state CSTRs," IFAC Proceedings Volumes, vol. 42, no. 11, pp. 697–702, Jan. 2009, https://doi.org/10.3182/20090712-4-TR-2008.00113. Engineering, Technology & Applied Science Research Vol. 11, No. 2, 2021, 6919-6929 6929 www.etasr.com Muhirwa et al.: Markov Chain Monte Carlo Analysis of the Variable-Volume Exothermic Model for … [53] S. Zhang, D. Muller, H. Arellano-Garcia, and G. Wozny, "CFD simulation of the fluid hydrodynamics in a continuous stirred-tank reactor," Chemical Engineering Transactions, vol. 32, pp. 1441–1446, Jun. 2013. [54] S. N. Naikwad and S. V. Dudul, "Identification of a Typical CSTR Using Optimal Focused Time Lagged Recurrent Neural Network Model with Gamma Memory Filter," Applied Computational Intelligence and Soft Computing, vol. 2009, Jan. 2010, Art. no. 385757, https://doi.org/10.1155/2009/385757. [55] S. R. Tofighi, F. Bayat, and F. Merrikh-Bayat, "Robust feedback linearization of an isothermal continuous stirred tank reactor: H∞ mixed- sensitivity synthesis and DK-iteration approaches," Transactions of the Institute of Measurement and Control, vol. 39, no. 3, pp. 344–351, Mar. 2017, https://doi.org/10.1177/0142331215603446. [56] T. Niederberger, "Markov Chain Monte Carlo Methods for Parameter Identification in Systems Biology Models," Ph.D. dissertation, Ludwig Maximilian University of Munich, Bad Reichenhall, Germany, 2012. [57] T. Rajagopalan and V. Seshadri, "Analysis of continuous stirred tank reactor as a multivariable process and algorithms for computer determination of the equilibrium states," International Journal of Control, vol. 15, no. 3, pp. 497–507, Mar. 1972, https://doi.org/10.1080/ 00207177208932165. [58] V. Nicoulaud-Gouin, L. Garcia-Sanchez, M. Giacalone, J. C. Attard, A. Martin-Garin, and F. Y. Bois, "Identifiability of sorption parameters in stirred flow-through reactor experiments and their identification with a Bayesian approach," Journal of Environmental Radioactivity, vol. 162– 163, pp. 328–339, Oct. 2016, https://doi.org/10.1016/j.jenvrad. 2016.06.008. [59] V. Roy, "Convergence Diagnostics for Markov Chain Monte Carlo," Annual Review of Statistics and Its Application, vol. 7, no. 1, pp. 387– 412, 2020, https://doi.org/10.1146/annurev-statistics-031219-041300. [60] Y. Lu, Z. Fang, and C. Gao, "Stabilization of (state, input)-disturbed CSTRs through the port-Hamiltonian systems approach," arXiv:1707.01560 [math], Jun. 2017, Accessed: Feb. 20, 2021. [Online]. Available: http://arxiv.org/abs/1707.01560. [61] Z. Prokopova and R. Prokop, "Modelling and Simulation of Chemical Industrial Reactors," in 23rd European Conference on Modelling and Simulation, Madrid, Spain, Jun. 2009, pp. 378–383.