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

Abstract

Reliable estimation of kinetic parameters in chemical systems comprising both slow and rapid reaction steps and rapidly reacting intermediate species is a difficult differential-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 differential equations and nonlinear algebraic equations is converted to ordinary differential 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.

1. 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.

2. 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.

2.1. 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,

$$\frac{d{C}_A}{dt}=-r_1$$

(1)

$$\frac{d{C}_R}{dt}=r_1-r_2$$

(2)

$$\frac{d{C}_S}{dt}=r_2$$

(3)

$$r_1=k_1\left(C_A-\frac{C_R}{K_1}\right)$$

(4)

$$r_2=k_2\left(C_R-\frac{C_S}{K_2}\right)$$

(5)

For the special case, where reaction step 1 is slow and step 2 is rapid, a reduced model can be derived:

$$\frac{C_S}{C_R}=K_2$$

(6)

$$\frac{d{C}_S}{dt}=K_2\frac{d{C}_R}{dt}$$

(7)

$$\frac{d{C}_A}{dt}+\frac{d{C}_R}{dt}+\frac{d{C}_S}{dt}=0$$

(8)

Substituting $\frac{d{C}_S}{dt}$ into Equation (8) yields

$$\frac{d{C}_A}{dt}+\left(1+K_2\right)\frac{d{C}_R}{dt}=0$$

(9)

$$\frac{d{C}_R}{dt}=-\frac{1}{1+K_2} \frac{d{C}_A}{dt}=+\frac{1}{1+K_2} r_1$$

(10)

$$\frac{d{C}_S}{dt}=+\frac{K_2}{1+K_2} r_1$$

(11)

As a summary we obtain for the special case of the consecutive reaction system, where step 1 is slow and step 2 is rapid:

$$\begin{gathered} \frac{d{C}_A}{dt}=-r_1 \\ \frac{d{C}_A}{dt}=-r_1 \\ \frac{d{C}_R}{dt}=\frac{1}{1+K_2} r_1 \\ \frac{d{C}_S}{dt}=\frac{K_2}{1+K_2} r_1 \\ {r}_1=k_1\left(C_A-\frac{C_R}{K_1}\right) \end{gathered}$$

The general case was compared with the special case by increasing the value of the kinetic constant k₂ of the rapid reaction step 2 to find the value of k₂ by which the solution of the general case approaches that of the special case.

By setting k₁ = 1 and assuming arbitrary values for the equilibrium constants K₁ = 10 and K₂ = 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₂ of the rapid reaction step 2. Other parameters k₁ = 1, K₁ = 10 and K₂ = 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₂ 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, k₂ = 1, 5, 10, 100 and 300, K₁ = 10 and K₂ = 10 compared to the special case are presented in

Table 1. The sum of squared residuals of the general case simulations as compared to the special case (k₁ = 1, K₁ = 10 and K₂ = 10).

k2	S
1	3.6113
5	0.2373
10	0.0644
100	0.006522
300	0.00071645

.

The behavior of the sum of squared residuals as a function k₂ is displayed in Figure 4. The general case approaches the special case as the rate constant k₂ of the rapid reaction step 2 exceeds 100 i.e., 100 times the value of the rate parameter k₁ of the slow reaction step 1.

2.2. 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.

The mass balances and rate equations of this reaction system can be written in the general case as follows,

$$\frac{d{C}_A}{dt}=-r_1-r_2$$

(12)

$$\frac{d{C}_R}{dt}=r_1$$

(13)

$$\frac{d{C}_S}{dt}=r_2$$

(14)

$$r_1=k_1\left(C_A-\frac{C_R}{K_1}\right)$$

(15)

$$r_2=k_2\left(C_A-\frac{C_S}{K_2}\right)$$

(16)

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:

$$\frac{d{C}_A}{dt}+\frac{d{C}_R}{dt}+\frac{d{C}_S}{dt}=0$$

(17)

$$\frac{C_S}{C_A}=K_2$$

(18)

$$\frac{d{C}_S}{dt}=K_2\frac{d{C}_A}{dt}$$

(19)

Substituting $\frac{d{C}_S}{dt}$ into Equation (17) yields:

$$\left(1+K_2\right)\frac{d{C}_A}{dt}+\frac{d{C}_R}{dt} =0$$

(20)

$$\frac{d{C}_A}{dt}=-\frac{1}{1+K_2} \frac{d{C}_R}{dt}=-\frac{1}{1+K_2} r_1$$

(21)

$$\frac{d{C}_S}{dt}=-\frac{K_2}{1+K_2} r_1$$

(22)

At time t = 0,

$$C_{R\left(0\right)}=0$$

(23)

$$\frac{C_{S\left(0\right)}}{C_{A\left(0\right)}-C_{S\left(0\right)}}=K_2$$

(24)

$$C_{S\left(0\right)}=\frac{K_2{C}_{0A}}{1+K_2}$$

(25)

$$C_{A\left(0\right)}=\frac{C_{0A}}{1+K_2}$$

(26)

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:

$$\begin{gathered} \frac{d{C}_A}{dt}=-\frac{1}{1+K_2} r_1 \\ \frac{d{C}_R}{dt}=r_1 \\ \frac{d{C}_S}{dt}=-\frac{K_2}{1+K_2} r_1 \\ {r}_1=k_1\left(C_A-\frac{C_R}{K_1}\right) \\ {C}_{A\left(0\right)}=\frac{C_{0A}}{1+K_2} \\ {C}_{R\left(0\right)}=0 \\ {C}_{S\left(0\right)}=\frac{K_2{C}_{0A}}{1+K_2} \end{gathered}$$

Initial values C_A(0) and C_S(0), C_R(0) = 0 as a function of equilibrium constant K₂ 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 k₁ = 1, K₁ = 10 and K₂ = 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 k₂ of the rapid reaction step 2 to find the value of k₂ 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₂ 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 k₂ while k₁ = 1, K₁ = 10 and K₂ = 2, (the number of simulated data points was 301 in the time interval of 0–15) are shown in

Table 2. The sum of squared residuals of the general case of parallel reaction system compared to the special case with different values of k₂ (k₁ = 1, K₁ = 10 and K₂ = 2).

k2	S
2	3.31901
5	1.68380
10	1.14954
50	0.89085
100	0.88929
1000	0.888893
10,000	0.888889

. 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₂ 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.

2.3. 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₂-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.

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₂ was constantly added to the system by keeping the pressure constant. Thus, the saturation concentration of CO₂ was presumed, and Henry’s law was applied to relate the partial pressure of CO₂ and the concentration of dissolved CO₂. 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.

2.3.1. Basic Mass Balances

The mass balances for bulk phase components in a batch reactor (assuming constant density) at a constant CO₂ concentration due to controlled pressure can be written as follows (ρ_B = mass of catalyst-to-liquid volume, i.e., the catalyst bulk density)

$$\frac{d\left[DMC\right]}{dt}=\left(r_3+r_5\right){\rho }_B$$

(27)

$$\frac{d\left[BG\right]}{dt}=r_3{\rho }_B$$

(28)

$$\frac{d\left[H_{2}O\right]}{dt}=r_5{\rho }_B$$

(29)

$$\frac{d\left[MeOH\right]}{dt}=-r_1{\rho }_B-r_4$$

(30)

$$\frac{d\left[BO\right]}{dt}=-r_4$$

(31)

$$\frac{d\left[MB\right]}{dt}=-r_3{\rho }_B+r_4$$

(32)

Applying a quasi-steady state to MC* and MeOH* gives relations (33) and (34),

$$r_{MC*}=r_2-r_3-r_5\approx 0$$

(33)

$$r_{MeOH*}=r_1-r_2-r_5\approx 0$$

(34)

Substituting Equations (33) and (34) into Equations (27)–(32) results in the following relationships:

$$\frac{d\left[DMC\right]}{dt}=r_2{\rho }_B$$

(35)

$$\frac{d\left[BG\right]}{dt}=r_3{\rho }_B$$

(36)

$$\frac{d\left[H_{2}O\right]}{dt}=r_5{\rho }_B$$

(37)

$$\frac{d\left[MeOH\right]}{dt}=-(r_2+r_5){\rho }_B-r_4$$

(38)

$$\frac{d\left[BO\right]}{dt}=-r_4$$

(39)

$$\frac{d\left[MB\right]}{dt}=-r_3{\rho }_B+r_4$$

(40)

2.3.2. Differential-Algebraic Problem

The addition of the mass balances further gives

$$\frac{d\left[MB\right]}{dt}+\frac{d\left[BG\right]}{dt}=r_4$$

(41)

$$\frac{d\left[MB\right]}{dt}+\frac{d\left[MeOH\right]}{dt}=-2r_2{\rho }_B$$

(42)

$$\frac{d\left[H_{2}O\right]}{dt}+\frac{d\left[BG\right]}{dt}=r_2{\rho }_B$$

(43)

The ODEs (Equations (33) and (35)–(37)) give the stoichiometric relations

$$\begin{gathered} {\displaystyle \underset{0}{\overset{\left[DMC\right]}{\int }}d\left[DMC\right]}={\displaystyle \underset{0}{\overset{\left[H_{2}O\right]}{\int }}d\left[H_{2}O\right]}+{\displaystyle \underset{0}{\overset{\left[BG\right]}{\int }}d\left[BG\right]} \\ \left[DMC\right]=\left[H_{2}O\right]+\left[BG\right] \end{gathered}$$

(44)

Equations (36), (39) and (40) give

$$\begin{gathered} -{\displaystyle \underset{\left[BO\right]o}{\overset{\left[BO\right]}{\int }}d\left[BO\right]}={\displaystyle \underset{0}{\overset{\left[BG\right]}{\int }}d\left[BG\right]}+{\displaystyle \underset{0}{\overset{\left[MB\right]}{\int }}d\left[MB\right]} \\ {\left[BO\right]}_O=\left[BO\right]+\left[BG\right]+\left[MB\right] \end{gathered}$$

(45)

From Equations (27), (30) and (32)–(34) is obtained

$$\begin{gathered} -2{\displaystyle \underset{0}{\overset{\left[DMC\right]}{\int }}d\left[DMC\right]}={\displaystyle \underset{\left[MeOH\right]o}{\overset{\left[MeOH\right]}{\int }}d\left[MeOH\right]}+{\displaystyle \underset{0}{\overset{\left[MB\right]}{\int }}d\left[MB\right]} \\ {\left[MeOH\right]}_O=\left[MeOH\right]+2\left[DMC\right]+\left[MB\right] \end{gathered}$$

(46)

The following relationship between the concentrations was obtained by consideration of the reaction mechanism [6] (α = K₁K₅/K₃),

$$\frac{\left[DMC\right]}{\left[BG\right]}=1+\alpha \frac{\left[MeOH\right]}{\left[MB\right]}$$

(47)

Equations (35), (39), (41)–(43) and (47) form a set of six differential-algebraic equations, with six unknowns: ([DMC], [BO], [BG], [H₂O], [MB], [MeOH]).

From Equation (47), the following expression can be obtained:

$$\frac{\left[BG\right]}{\left[DMC\right]}=\frac{1}{1+\alpha \frac{\left[MeOH\right]}{\left[MB\right]}}=\frac{\left[MB\right]}{\left[MB\right]+\alpha \left[MeOH\right]}$$

(48)

From Equations (44) and (47), the following relationship is obtained:

$$\frac{\left[H_{2}O\right]}{\left[DMC\right]}=1-\frac{\left[BG\right]}{\left[DMC\right]}$$

(49)

$$\frac{\left[H_{2}O\right]}{\left[DMC\right]}=\frac{\alpha \left[MeOH\right]}{\left[MB\right]+\alpha \left[MeOH\right]}$$

(50)

Equations (45) and (48) give the following relation:

$$\frac{\left[MB\right]\left[DMC\right]}{\left[MB\right]+\alpha \left[MeOH\right]}+\left[MB\right]={\left[BO\right]}_O-\left[BO\right]=\omega$$

(51)

i.e.,

$$\left[MB\right]\left[DMC\right]+{\left[MB\right]}^2+\alpha \left[MeOH\right]\left[MB\right]=\omega \left[MB\right]+\alpha \left[MeOH\right]\omega$$

(52)

i.e.,

$${\left[MB\right]}^2+\left(\left[DMC\right]-\omega \right)\left[MB\right]+\left(\alpha \left[MB\right]-\alpha \omega \right)({\left[MeOH\right]}_O-2\left[DMC\right]-\left[MB\right])=0$$

(53)

Rearranging of Equation (53) gives

$$\left(1-\alpha \right){\left[MB\right]}^2+\left(\left[DMC\right]-\omega +\alpha \gamma +\alpha \omega \right)\left[MB\right]-\alpha \gamma \omega =0$$

(54)

where $\gamma ={\left[MeOH\right]}_O-2\left[DMC\right]$. The solution of the second-degree Equation (54) becomes

$$\left[MB\right]=\frac{\beta \pm \sqrt{{\beta }^2+4\left(1-\alpha \right)\alpha \gamma \omega }}{2\left(1-\alpha \right)}$$

(55)

where $\beta =\omega -\left[DMC\right]-\alpha \left(\gamma +\omega \right)$. The sign + in the nominator gives the physically meaningful solution of Equation (55).

Equations (46), (48) and (50) give

$$\left[MeOH\right]=\gamma -\left[MB\right]$$

(56)

$$\left[BG\right]=\frac{\left[MB\right]\left[DMC\right]}{\left[MB\right]+\alpha \left[MeOH\right]}$$

(57)

$$\left[H_{2}O\right]=\frac{\alpha \left[MeOH\right]\left[DMC\right]}{\left[MB\right]+\alpha \left[MeOH\right]}$$

(58)

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.

2.3.3. 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 (48),

$$\frac{\left[BG\right]}{\left[DMC\right]}=\frac{\frac{MB}{MeOH}}{\frac{MB}{MeOH}+\alpha }$$

(59)

$$\left[BG\right]=\frac{y}{\alpha +y}\left[DMC\right]=f(y)\left[DMC\right]$$

(60)

where $y=\frac{\left[MB\right]}{\left[MeOH\right]}$.

The procedure is continued by differentiation as shown below,

$$\frac{df}{dt}=\frac{\frac{dy}{dt}\left(\alpha +y\right)-\frac{dy}{dt}\left(y\right)}{{\left(\alpha +y\right)}^2}=\frac{\alpha \frac{dy}{dt}}{{\left(\alpha +y\right)}^2}$$

(61)

$$\frac{d\left[BG\right]}{dt}=\frac{df}{dt}\left[DMC\right]+f\frac{d\left[DMC\right]}{dt}$$

(62)

$$\frac{d\left[BG\right]}{dt}=\frac{\alpha \left[DMC\right]}{{\left(\alpha +y\right)}^2}\frac{dy}{dt}+\frac{y}{\alpha +y}\frac{d\left[DMC\right]}{dt}$$

(63)

The definition y gives the relation

$$\frac{d\left[MB\right]}{dt}=\frac{dy}{dt}\left[MeOH\right]+y\frac{d\left[MeOH\right]}{dt}$$

(64)

$$\frac{dy}{dt}=\frac{\frac{d\left[MB\right]}{dt}}{\left[MeOH\right]}-y\frac{\frac{d\left[MeOH\right]}{dt}}{\left[MeOH\right]}$$

(65)

Equation (16) gives

$$\frac{d\left[MeOH\right]}{dt}=-\frac{d\left[MB\right]}{dt}-2r_2{\rho }_B$$

(66)

A combination of Equations (65) and (66) gives

$$\frac{dy}{dt}=\frac{\frac{d\left[MB\right]}{dt}}{\left[MeOH\right]}+y\frac{\left(\frac{d\left[MB\right]}{dt}+2r_2{\rho }_B\right)}{\left[MeOH\right]}$$

(67)

Equations (41) and (63) give

$$\begin{array}{l}\frac{d\left[BG\right]}{dt}=r_4-\frac{d\left[MB\right]}{dt}\\ =\frac{\alpha \left[DMC\right]}{{\left(\alpha +y\right)}^2}\frac{dy}{dt}+\frac{y}{\alpha +y}{r}_2{\rho }_B\end{array}$$

(68)

From Equation (68) dy/dt is solved:

$$\begin{array}{l}\frac{dy}{dt}=\frac{{\left(\alpha +y\right)}^2}{\alpha \left[DMC\right]}\\ \left(r_4-\frac{y}{\alpha +y}{r}_2{\rho }_B-\frac{d\left[MB\right]}{dt}\right)\end{array}$$

(69)

Equations (67) and (69) are equal, therefore,

$$\begin{array}{l}\frac{\left(1+y\right)}{\left[MeOH\right]}\frac{d\left[MB\right]}{dt}+2r_2{\rho }_B\frac{y}{\left[MeOH\right]}\\ =\frac{{\left(\alpha +y\right)}^2}{\alpha \left[DMC\right]}\left(r_4-\frac{y}{\alpha +y}{r}_2{\rho }_B\right)-\frac{{\left(\alpha +y\right)}^2}{\alpha \left[DMC\right]}\frac{d\left[MB\right]}{dt}\end{array}$$

(70)

from which the time derivative is solved explicitly:

$$\begin{array}{l}\left[\frac{\left(1+y\right)\left[DMC\right]}{\left[MeOH\right]}+\frac{{\left(\alpha +y\right)}^2}{\alpha }\right]\frac{d\left[MB\right]}{dt}\\ =-\left[\frac{\left(\alpha +y\right)y}{\alpha }-\frac{2y\left[DMC\right]}{\left[MeOH\right]}\right]r_2{\rho }_B+\frac{{\left(\alpha +y\right)}^2}{\alpha }{r}_4\end{array}$$

(71)

$$\begin{array}{l}\left[\frac{\alpha \left(1+y\right)\left[DMC\right]}{{\left(\alpha +y\right)}^2\left[MeOH\right]}+1\right]\frac{d\left[MB\right]}{dt}\\ =r_4-\frac{y}{\alpha +y}\left[1-\frac{2\alpha }{\alpha +y}\frac{\left[DMC\right]}{\left[MeOH\right]}\right]r_2{\rho }_B\end{array}$$

(72)

The final forms of the mass balances for the liquid-phase components in the batch reactor thus become (α = constant):

$$\frac{d\left[DMC\right]}{dt}=r_2{\rho }_B$$

(73)

$$\frac{d\left[BO\right]}{dt}=-r_4$$

(74)

$$\begin{gathered} \frac{d\left[MB\right]}{dt}={\omega }^{-1}\left\{r_4-\left[\frac{y}{\alpha +y}\right]\left[1-\frac{2\alpha }{\alpha +y}\frac{\left[DMC\right]}{\left[MeOH\right]}\right]r_2{\rho }_B\right\} \\ \omega =\left[\frac{\alpha \left(1+y\right)\left[DMC\right]}{{\left(\alpha +y\right)}^2\left[MeOH\right]}+1\right]y=\frac{\left[MB\right]}{\left[MeOH\right]} \end{gathered}$$

(75)

$$\frac{d\left[BG\right]}{dt}=r_4-\frac{d\left[MB\right]}{dt}$$

(76)

$$\frac{d\left[H_{2}O\right]}{dt}=r_2{\rho }_B-r_4+\frac{d\left[MB\right]}{dt}$$

(77)

$$\frac{d\left[MeOH\right]}{dt}=-2r_2{\rho }_B-\frac{d\left[MB\right]}{dt}$$

(78)

The rates of the rate-determining steps (r₂ and r₄) were obtained from the mechanism—these rates include only the concentrations of CO₂, MeOH, DMC, MB, BO, BG and H₂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:

$$r_2=\frac{k_2\left(\left[MeOH\right]P_{C{O}_2}-\frac{1}{K_{}}\frac{\left[DMC\right]\left[BG\right]}{\left[MB\right]}\right)}{1+K_1\left[MeOH\right]+\frac{1}{K_3}\frac{\left[DMC\right]\left[BG\right]}{\left[MB\right]}}$$

(79)

$$r_4=k_4\left[BO\right]\left[MeOH\right]$$

(80)

2.3.4. Parameter Estimation Results

Rate parameters k₂ and k₄ 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₂ initially. Parameters K₁ and 1/K₃ (

Table 3. Parameter estimation results.

Parameter	Value		Error/%
k2		1.60 × 10−5	6.8
k4		1.12 × 10−2	6.7
K1 = 0, 1/K3 = 0, α = 1.3 × 10−3, K = 0.08 × 10−5, at 150 °C

, Equations (47) and (48)) turned out not to be significant and were approximated to zero in the rate expression r₂. Parameter α was determined separately from the plot according to Equation (47); α = 1.3 × 10⁻³. The thermodynamic equilibrium constant was estimated from theoretical calculations: K = 0.08 × 10⁻⁵, 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%.

3. 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.