A Robust Method for the Estimation of Kinetic Parameters for Systems Including Slow and Rapid Reactions—From Di ﬀ erential-Algebraic Model to Di ﬀ erential Model

: Reliable estimation of kinetic parameters in chemical systems comprising both slow and rapid reaction steps and rapidly reacting intermediate species is a di ﬃ cult di ﬀ erential-algebraic problem. Consequently, any conventional approach easily leads to serious convergence and stability problems during the parameter estimation. A robust method is proposed to surmount this dilemma: the system of ordinary di ﬀ erential equations and nonlinear algebraic equations is converted to ordinary di ﬀ erential equations, which are solved in-situ during the parameter estimation. The approach was illustrated with two generic examples and an example from green chemistry: synthesis of dimethyl carbonate from carbon dioxide and methanol.


Introduction
A reliable estimation of kinetic and thermodynamic parameters is one of the most important tasks in chemical reaction engineering. Only in few cases, like in the case of some homogeneous gas-phase reactions, is it possible to determine the kinetic constants a priori, exclusively from theoretical calculations. Consequently, in most cases, an extensive matrix of experimental work is needed. Particularly, in the presence of heterogeneous catalysts, small variations in the chemical composition and physical structure of the catalyst can change the reaction kinetics, and a new experimental program is inevitable. Typically, kinetic experiments are carried out in batch reactors or in continuous reactors with a well-established flow pattern, i.e., perfect back-mixing or plug flow. After that, the kinetic and thermodynamic parameters are determined by non-linear regression analysis. In the most recent 30 years, numerical methods and computing power have advanced tremendously. The development of solvers for stiff differential equations has enabled the treatment of difficult problems in chemical kinetics and simulation of large chemical systems appearing in combustion and atmospheric processes. However, even a superficially small problem can become a difficult one, as the system usually consists of a set of slow and rapid reaction steps. From a mathematical viewpoint, the system of ordinary differential equations describing the behavior of the components in a batch reactor becomes a system of differential-algebraic equations (DAE), for which solvers have been developed in recent years. In principle, the task could be solved by coupling a DAE solver to a parameter estimator (a nonlinear regression routine). However, if the parameters glide to an unrealistic regime in the course of the parameter estimation, significant problems arise: the solutions of the algebraic equations induce problems, such as negative concentrations, which usually leads to the collapse of the whole computational process. Our experience with different case studies [1] has taught that a more robust approach is to convert the DAE problem to a set of ordinary differential equations (ODEs) and then solve these ODEs as a sub-problem in the parameter estimation. Here we illustrate the method with two generic examples and a highly relevant issue in chemical technology, namely the utilization of carbon dioxide as a raw material for the synthesis of green chemicals.

Development of a Robust Method
The tasks in the method are briefly summarized as follows: T1. Mass balances of all components in the batch reactor are written down (ODEs) T2. A quasi-steady-state hypothesis is applied to the intermediates-the ODEs become a system of DAEs is created T3. The DAE system is converted to a set of implicit ODEs by differentiation T4. The system of implicit ODEs is converted to a set of explicit ODEs T5. The system of explicit ODEs is implemented in a combined stiff ODE solver-parameter estimation software Modest.

Example: Consecutive Reactions with Slow and Rapid Steps
The following consecutive reaction scheme is considered in Figure 1.
Step 1 is presumed to be slow whereas step 2 is rapid. parameter estimator (a nonlinear regression routine). However, if the parameters glide to an unrealistic regime in the course of the parameter estimation, significant problems arise: the solutions of the algebraic equations induce problems, such as negative concentrations, which usually leads to the collapse of the whole computational process. Our experience with different case studies [1] has taught that a more robust approach is to convert the DAE problem to a set of ordinary differential equations (ODEs) and then solve these ODEs as a sub-problem in the parameter estimation. Here we illustrate the method with two generic examples and a highly relevant issue in chemical technology, namely the utilization of carbon dioxide as a raw material for the synthesis of green chemicals.

Development of a Robust Method
The tasks in the method are briefly summarized as follows: T1. Mass balances of all components in the batch reactor are written down (ODEs) T2. A quasi-steady-state hypothesis is applied to the intermediates-the ODEs become a system of DAEs is created T3. The DAE system is converted to a set of implicit ODEs by differentiation T4. The system of implicit ODEs is converted to a set of explicit ODEs T5. The system of explicit ODEs is implemented in a combined stiff ODE solver-parameter estimation software Modest.

Example: Consecutive Reactions with Slow and Rapid Steps
The following consecutive reaction scheme is considered in Figure 1.
Step 1 is presumed to be slow whereas step 2 is rapid. The mass balances and rate equations of this consecutive reaction system can be written in the general case as follows, For the special case, where reaction step 1 is slow and step 2 is rapid, a reduced model can be derived: The mass balances and rate equations of this consecutive reaction system can be written in the general case as follows, For the special case, where reaction step 1 is slow and step 2 is rapid, a reduced model can be derived: As a summary we obtain for the special case of the consecutive reaction system, where step 1 is slow and step 2 is rapid: The general case was compared with the special case by increasing the value of the kinetic constant k 2 of the rapid reaction step 2 to find the value of k 2 by which the solution of the general case approaches that of the special case.
By setting k 1 = 1 and assuming arbitrary values for the equilibrium constants K 1 = 10 and K 2 = 10, the concentration profiles depicted in Figure 2 were obtained for the special case.
Substituting into Equation (8) yields As a summary we obtain for the special case of the consecutive reaction system, where step 1 is slow and step 2 is rapid: The general case was compared with the special case by increasing the value of the kinetic constant k2 of the rapid reaction step 2 to find the value of k2 by which the solution of the general case approaches that of the special case.
By setting k1 = 1 and assuming arbitrary values for the equilibrium constants K1 = 10 and K2 = 10, the concentration profiles depicted in Figure 2 were obtained for the special case.  As revealed by Figure 2, the system behaves like a parallel reaction system, which is caused by the fact that as soon as R is formed, some part of it is immediately transformed to S.
The general case was simulated by increasing the value of the rate parameter k 2 of the rapid reaction step 2. Other parameters k 1 = 1, K 1 = 10 and K 2 = 10 were kept the same as in the special case. The sum of squared residuals was calculated for the difference between simulated data points of the special case and the general case (number of simulated data points was 101 in the time frame of 0-10 in each case). Concentration profiles for the general case of consecutive reaction system corresponding to k 2 values of 1, 5, 10, 100 and 300 are depicted in Figure 3.
The sums of squared residuals of the general case simulations of consecutive reaction system with k 1 = 1, k 2 = 1, 5, 10, 100 and 300, K 1 = 10 and K 2 = 10 compared to the special case are presented in Table 1. As revealed by Figure2, the system behaves like a parallel reaction system, which is caused by the fact that as soon as R is formed, some part of it is immediately transformed to S.
The general case was simulated by increasing the value of the rate parameter k2 of the rapid reaction step 2. Other parameters k1 = 1, K1 = 10 and K2 = 10 were kept the same as in the special case. The sum of squared residuals was calculated for the difference between simulated data points of the special case and the general case (number of simulated data points was 101 in the time frame of 0-10 in each case). Concentration profiles for the general case of consecutive reaction system corresponding to k2 values of 1, 5, 10, 100 and 300 are depicted in Figure 3. The sums of squared residuals of the general case simulations of consecutive reaction system with k1 = 1, k2 = 1, 5, 10, 100 and 300, K1 = 10 and K2 = 10 compared to the special case are presented in Table 1. Table 1. The sum of squared residuals of the general case simulations as compared to the special case (k1 = 1, K1 = 10 and K2 = 10). The behavior of the sum of squared residuals as a function k2 is displayed in Figure 4. The general case approaches the special case as the rate constant k2 of the rapid reaction step 2 exceeds 100 i.e., 100 times the value of the rate parameter k1 of the slow reaction step 1. The behavior of the sum of squared residuals as a function k 2 is displayed in Figure 4. The general case approaches the special case as the rate constant k 2 of the rapid reaction step 2 exceeds 100 i.e., 100 times the value of the rate parameter k 1 of the slow reaction step 1.

Parallel Reactions with Slow and Rapid Steps
The second example is the classical parallel reaction scheme displayed in Figure 5.
Step 1 is presumed to be slow while step 2 is rapid.
A R 1 Figure 4. Sum of squared residuals as a function of the rate parameter k 2 of the rapid reaction step 2 in the general case of consecutive reaction system compared to the special case (k 1 = 1, K 1 = 10 and K 2 = 10).

Parallel Reactions with Slow and Rapid Steps
The second example is the classical parallel reaction scheme displayed in Figure 5.
Step 1 is presumed to be slow while step 2 is rapid. Figure 4. Sum of squared residuals as a function of the rate parameter k2 of the rapid reaction step 2 in the general case of consecutive reaction system compared to the special case (k1 = 1, K1 = 10 and K2 = 10).

Parallel Reactions with Slow and Rapid Steps
The second example is the classical parallel reaction scheme displayed in Figure 5.
Step 1 is presumed to be slow while step 2 is rapid.
For the special case, where reaction step 1 is slow and step 2 is rapid, a reduced model can be derived for the parallel reaction system: The mass balances and rate equations of this reaction system can be written in the general case as follows, For the special case, where reaction step 1 is slow and step 2 is rapid, a reduced model can be derived for the parallel reaction system: Substituting dC S dt into Equation (17) yields: The initial conditions (25) and (26) arise from the fact that some amounts of component S is formed immediately in the system, because step 2 progresses with an infinite rate. As a summary we obtain for the special case of the parallel reaction system, where step 1 is slow and step 2 is rapid: Initial values C A (0) and C S (0), C R (0) = 0 as a function of equilibrium constant K 2 in the special case of a parallel reaction system are depicted in Figure 6. Some concentration profiles as an example for the special case of a parallel reaction system with parameters k1 = 1, K1 = 10 and K2 = 2 are displayed in Figure 7. Some concentration profiles as an example for the special case of a parallel reaction system with parameters k 1 = 1, K 1 = 10 and K 2 = 2 are displayed in Figure 7. Figure 6. Initial values CA(0) and CS(0), CR(0) = 0 as a function of equilibrium constant K2 in the special case of a parallel reaction system, where step 1 is slow and step 2 is rapid.
Some concentration profiles as an example for the special case of a parallel reaction system with parameters k1 = 1, K1 = 10 and K2 = 2 are displayed in Figure 7. The general case of a parallel reaction system was compared with the special case as in the previous example 2.1 by increasing the value of the rate parameter k2 of the rapid reaction step 2 to find the value of k2 by which the solution of the general case approaches that of the special case. Concentration profiles for the general case of parallel reactions corresponding to k2 values of 2, 10 and 50 are depicted in Figure 8. The general case of a parallel reaction system was compared with the special case as in the previous example 2.1 by increasing the value of the rate parameter k 2 of the rapid reaction step 2 to find the value of k 2 by which the solution of the general case approaches that of the special case. Concentration profiles for the general case of parallel reactions corresponding to k 2 values of 2, 10 and 50 are depicted in Figure 8. The sum of squared residuals calculated as in the example 2.1 between simulated data points of the special case and the general case with different values of the rate parameter k2 while k1 = 1, K1 = 10 and K2 = 2, (the number of simulated data points was 301 in the time interval of 0-15) are shown in Table 2. The behavior of the sum of squared residuals as a function is displayed in Figure 9. The general case approaches the special case in this example when k2 is >50. The sum of squared residuals calculated as in the example 2.1 between simulated data points of the special case and the general case with different values of the rate parameter k 2 while k 1 = 1, K 1 = 10 and K 2 = 2, (the number of simulated data points was 301 in the time interval of 0-15) are shown in Table 2. The behavior of the sum of squared residuals as a function is displayed in Figure 9. The general case approaches the special case in this example when k 2 is >50. This kind of parallel reaction system has been studied by Branco et al. [2], who developed a concept for the determination of the switching point between kinetic and thermodynamic control for cases when one of the reactions is very rapid, whereas the second one is thermodynamically favored. The same phenomenon can be seen in Figure 6: S is the kinetically favored product, with a rapidly increasing concentration in the beginning, but with an increasing reaction time the concentration of R becomes higher because it is favored by thermodynamics, i.e., its formation has a higher equilibrium constant. The sum of squared residuals calculated as in the example 2.1 between simulated data points of the special case and the general case with different values of the rate parameter k2 while k1 = 1, K1 = 10 and K2 = 2, (the number of simulated data points was 301 in the time interval of 0-15) are shown in Table 2. The behavior of the sum of squared residuals as a function is displayed in Figure 9. The general case approaches the special case in this example when k2 is >50.  Figure 9. The sum of the squared residuals as a function of the rate parameter k 2 of the fast reaction step 2 in the general case of a parallel reaction system compared to the special case (k 1 = 1, K 1 = 10 and K 2 = 2).

Example: Synthesis of Dimethyl Carbonate from Methanol and Carbon Dioxide
Carbon dioxide reacts with methanol (MeOH) to yield dimethyl carbonate (DMC), which is a green alternative to methyl tert-butyl ether (MTBE) and can be used as a carbonylating and methylating agent [3,4]. The presence of a heterogeneous catalyst is required for the synthesis of DMC. We used zirconia-based catalysts (ZrO 2 -MgO) in an isothermal and isobaric laboratory-scale batch reactor [5]. The thermodynamics for this reaction is extremely unfavorable, so a way to shift the equilibrium is to include an additive to the reaction mixture; the role of the additive was to act as a chemical dehydration agent, i.e., to capture the water formed in the reaction. Butylene oxide (BO) was selected as the additive [5]. In this way, the process gains a more irreversible character and can be forced to the side of the products. Methylene butylate (MB) appears as an intermediate species, forming butylene glycol (BG) and thus preventing the water formation. The reaction scheme is displayed below in Figure 10, where * denotes a vacant site on the surface of the catalyst, and MeOH* denotes adsorbed methanol on the catalyst surface.
Processes 2020, 8, 1552 9 of 16 a chemical dehydration agent, i.e., to capture the water formed in the reaction. Butylene oxide (BO) was selected as the additive [5]. In this way, the process gains a more irreversible character and can be forced to the side of the products. Methylene butylate (MB) appears as an intermediate species, forming butylene glycol (BG) and thus preventing the water formation. The reaction scheme is displayed below in Figure 10, where * denotes a vacant site on the surface of the catalyst, and MeOH* denotes adsorbed methanol on the catalyst surface. Steps 1-3 and 5 take place on the catalyst surface, while step 4 proceeds as a homogeneous liquid-phase reaction. Reaction steps 2 and 4 in the scheme were assumed to be rate determining, whereas the adsorption step of methanol was presumed to be rapid. In addition, steps 3 and 5 were taken as rapid steps compared to steps 2 and 4. A constant-density system was assumed and the mass balance of carbon dioxide was not included, since CO2 was constantly added to the system by keeping the pressure constant. Thus, the saturation concentration of CO2 was presumed, and Henry's law was applied to relate the partial pressure of CO2 and the concentration of dissolved CO2. A large excess of methanol was used. Based on the reaction mechanism displayed above, the rate equations for the rate determining steps were derived. The details of the derivation of the rate equations are given as supplementary material in Appendix A: Derivation of the Rate Equations. Steps 1-3 and 5 take place on the catalyst surface, while step 4 proceeds as a homogeneous liquid-phase reaction. Reaction steps 2 and 4 in the scheme were assumed to be rate determining, whereas the adsorption step of methanol was presumed to be rapid. In addition, steps 3 and 5 were taken as rapid steps compared to steps 2 and 4. A constant-density system was assumed and the mass balance of carbon dioxide was not included, since CO 2 was constantly added to the system by keeping the pressure constant. Thus, the saturation concentration of CO 2 was presumed, and Henry's law was applied to relate the partial pressure of CO 2 and the concentration of dissolved CO 2 . A large excess of methanol was used. Based on the reaction mechanism displayed above, the rate equations for the rate determining steps were derived. The details of the derivation of the rate equations are given as supplementary material in Appendix A: Derivation of the Rate Equations.

Basic Mass Balances
The mass balances for bulk phase components in a batch reactor (assuming constant density) at a constant CO 2 concentration due to controlled pressure can be written as follows (ρ B = mass of catalyst-to-liquid volume, i.e., the catalyst bulk density) Applying a quasi-steady state to MC* and MeOH* gives relations (33) and (34), Substituting Equations (33) and (34) into Equations (27)-(32) results in the following relationships:

Differential-Algebraic Problem
The addition of the mass balances further gives The ODEs (Equations (33) and (35)-(37)) give the stoichiometric relations The following relationship between the concentrations was obtained by consideration of the reaction mechanism [6] (α = K 1 K 5 /K 3 ),

[DMC]
[BG] i.e., Rearranging of Equation (53) gives . The solution of the second-degree Equation (54) becomes A system of two ODEs Equations (35) and (39) can be solved by employing Equations (55)-(58). Thus, the problem can in principle be solved as an ODE problem coupled to algebraic equations.

Transformation to ODEs
In order to obtain a more robust algorithm for parameter estimation, the differential-algebraic problem is transformed to ODEs as follows.
From Equation ( The rates of the rate-determining steps (r 2 and r 4 ) were obtained from the mechanism-these rates include only the concentrations of CO 2 , MeOH, DMC, MB, BO, BG and H 2 O, respectively (Appendix A: Derivation of the Rate Equations). The rate equations are given below. The ODEs were solved with the backward difference method during the parameter estimation, which was performed with the Levenberg-Marquardt algorithm [7].
Rate equations:

Parameter Estimation Results
Rate parameters k 2 and k 4 were estimated from isothermal experiments shown as Equations (79) and (80). The experimental temperature was 150 • C, and the pressure was 45 bar of CO 2 initially. Parameters K 1 and 1/K 3 (Table 3, Equations (47) and (48)) turned out not to be significant and were approximated to zero in the rate expression r 2 . Parameter α was determined separately from the plot according to Equation (47); α = 1.3 × 10 −3 . The thermodynamic equilibrium constant was estimated from theoretical calculations: K = 0.08 × 10 −5 , at 150 • C. Table 3. Parameter estimation results.

Parameter
Value Error/% k 2 1.60 × 10 −5 6.8 k 4 1.12 × 10 −2 6.7 K 1 = 0, 1/K 3 = 0, α = 1.3 × 10 −3 , K = 0.08 × 10 −5 , at 150 • C The solution of ODEs and parameter estimation progressed very smoothly and led to an excellent description of the experimental data. The errors of the kinetic parameters were clearly less than 10%. The performance of the procedure is illustrated in Figure 11 The numerical values of the estimated parameters are enlisted in Table 3. The overall degree of explanation was 99.98%.

Conclusions
Solution of a differential-algebraic problem in connection to both fast and slow reaction steps is a demanding task when coupled to a parameter estimation task. In order to surmount the numerical problems often appearing for DAE systems, we propose a robust procedure, which implies the transformation of the DAE system to a set of explicit ODEs that can be easily solved in-situ during the estimation of kinetic parameters. The success of the methodology was illustrated with two generic examples and a case study, synthesis of dimethyl carbonate (DMC) from methanol and carbon dioxide.

Conclusions
Solution of a differential-algebraic problem in connection to both fast and slow reaction steps is a demanding task when coupled to a parameter estimation task. In order to surmount the numerical problems often appearing for DAE systems, we propose a robust procedure, which implies the transformation of the DAE system to a set of explicit ODEs that can be easily solved in-situ during the estimation of kinetic parameters. The success of the methodology was illustrated with two generic examples and a case study, synthesis of dimethyl carbonate (DMC) from methanol and carbon dioxide.