DOW Experimental Design

DOW Experimental Design
State dimension:	1
Differential states:	11
Discrete control functions:	2
Path constraints:	4
Interior point equalities:	11

The DOW Experimental Design problem models the OED problem for the parameter estimation problem formulated by the DOW Chemical Co. in 1981. The following formulation is taken from [1].

Chemical background

The chemical species are disguised for proprietary reasons and the desired reaction is given by $HA+2BM\to AB+HBMH$, where $AB$ is the desired product. The reactions are described as follows:

Slow Kinetic Reactions:

$\begin{array}{c}{M}^{-}+BM\underset{k_{-1}}{\leftarrow }\stackrel{k_1}{\to }MB{M}^{-}\\ A^{-}+BM\stackrel{k_2}{\to }AB{M}^{-}\\ M^{-}+AB\underset{k_{-3}}{\leftarrow }\stackrel{k_3}{\to }AB{M}^{-}\end{array}$

Acid-Base Reactions:

$\begin{array}{c}MBMH\leftarrow K_1\to MB{M}^{-}+H^{+}\\ HA\leftarrow K_2\to A^{-}+H^{+}\\ HABM\leftarrow K_3\to AB{M}^{-}+H^{+}\end{array}$

In order to devise a model to account for these reactions, it is first necessary to distinguish between the overall concentration of a species and the concentration of its neutral form. Overall concentrations are defined for three components based on neutral and ionic species

$\begin{array}{rl}\left[HBMH\right]& =[(MBMH{)}_N]+\left[MB{M}^{-}\right]\\ \left[HA\right]& =[(HA{)}_N]+\left[A^{-}\right]\\ \left[HABM\right]& =[(HABM{)}_N]+\left[AB{M}^{-}\right]\end{array}$

Here $[]$ denotes the concentration of the species in $mol/kg$. By assuming the rapid acid-base reactions are at equilibrium, the equilibrium constants $K_1,K_2,K_3$ can be defined as

$\begin{array}{rl}{K}_1& =\frac{[MB{M}^{-}][H^{+}]}{[(MBMH{)}_N]}\\ K_2& =\frac{[A^{-}][H^{+}]}{[(HA{)}_N]}\\ K_3& =\frac{[AB{M}^{-}][H^{+}]}{[(HABM{)}_N]}\end{array}$

The anionic species may then be represented by

$\begin{array}{rlrl}\left[MB{M}^{-}\right]& =\frac{K_1[MBMH]}{K_1+[H^{+}]}& & (a)\\ \left[A^{-}\right]& =\frac{K_2[HA]}{K_2+[H^{+}]}& & (b)\\ \left[AB{M}^{-}\right]& =\frac{K_3[HABM]}{K_3+[H^{+}]}& & (c)\end{array}$

Material balance equations for the three reactants in the slow kinetic reactions yield:

$\begin{array}{rlrl}\frac{d[M^{-}]}{dt}& =-k_1\left[M^{-}\right]\left[BM\right]+k_{-1}\left[MB{M}^{-}\right]-k_3\left[M^{-}\right]\left[AB\right]+k_{-1}\left[AB{M}^{-}\right]& & (d)\\ \frac{d[BM]}{dt}& =-k_1\left[M^{-}\right]\left[BM\right]+k_{-1}\left[MB{M}^{-}\right]-k_2\left[A^{-}\right]\left[BM\right]& & (e)\\ \frac{d[AB]}{dt}& =-k_3\left[M^{-}\right]\left[AB\right]+k_{-3}\left[AB{M}^{-}\right]& & (f)\end{array}$

From stoichiometry, rate expressions can also be written for the total species:

$\begin{array}{rlrl}\frac{d[MBMH]}{dt}& =k_1\left[M^{-}\right]\left[BM\right]-k_{-1}\left[MB{M}^{-}\right]& & (g)\\ \frac{d[HA]}{dt}& =k_2\left[A^{-}\right]\left[BM\right]& & (h)\\ \frac{d[HABM]}{dt}& =k_2\left[A^{-}\right]\left[BM\right]+k_3\left[M^{-}\right]\left[AB\right]-k_{-3}\left[AB{M}^{-}\right]& & (i)\end{array}$

An electroneutrality constraint gives the hydrogen ion concentration $\left[H^{+}\right]$ as

$\left[H^{+}\right]+\left[Q^{+}\right]=\left[M^{-}\right]+\left[MB{M}^{-}\right]+\left[A^{-}\right]+\left[AB{M}^{-}\right](j)$

Based on similarities of reacting species, we assume

$k_3=k_1,k_{-3}=\frac{1}{2}{k}_{-1}$

With these assumptions, the three rate constants $k_1,k_2$ and $k_{-1}$ must be estimated. Each of these can be fitted with two adjustable model parameters, assuming an Arrhenius temperature dependence. That is

$\begin{array}{c}{k}_i={\alpha }_i\exp (-E_i/(RT)),i\in \{-1,1,2\}\end{array}$

Here $R\approx 1.98720425864083cal/(K\cdot mol)$ is the gas constant and $T$ is the reaction temperature in Kelvins. The parameters ${\alpha }_i$, given in $mol/(kg\cdot h)$, represent the pre-exponential factors and the $E_i$, with unit $cal/mol$, are the activation energies.

Mathematical formulation

The chemical processes $(a)-(j)$ can be expressed mathematically as six differential equations and four algebraic equations:

$\begin{array}{rlrl}{\dot{y}}_1& =-k_2{y}_8{y}_2& & (1),(h)\\ {\dot{y}}_2& =-k_1{y}_6{y}_2+k_{-1}{y}_{10}-k_2{y}_8{y}_2& & (2),(e)\\ {\dot{y}}_3& =k_2{y}_8{y}_2+k_1{y}_6{y}_4-\frac{1}{2}{k}_{-1}{y}_9& & (3),(i)\\ {\dot{y}}_4& =-k_1{y}_6{y}_4+\frac{1}{2}{k}_{-1}{y}_9& & (4),(f)\\ {\dot{y}}_5& =k_1{y}_6{y}_2-k_{-1}{y}_{10}& & (5),(g)\\ {\dot{y}}_6& =-k_1(y_6{y}_2+y_6{y}_4)+k_{-1}(y_{10}+\frac{1}{2}{y}_9)& & (6),(d)\\ y_7& =-\left[Q^{+}\right]+y_6+y_8+y_9+y_{10}& & (7),(j)\\ y_8& =\frac{{\theta }_8{y}_1}{{\theta }_8+y_7}& & (8),(b)\\ y_9& =\frac{{\theta }_9{y}_3}{{\theta }_9+y_7}& & (9),(c)\\ y_{10}& =\frac{{\theta }_7{y}_5}{{\theta }_7+y_7}& & (10),(a)\end{array}$

Here the letters in parentheses stand for the corresponding chemical process and the quantity $\left[Q^{+}\right]=0.0131$ is a constant during the reaction. The nine parameters form the vector

$\theta =({\alpha }_1,E_1,{\alpha }_2,E_2,{\alpha }_{-1},E_{-1},K_1,K_2,K_3)$

The predicted concentrations form the vector

$y=(HA,BM,HABM,AB,MBMH,M^{-},H^{+},A^{-},AB{M}^{-},MB{M}^{-})$

Let $f_k(\cdot )$ denote the right hand side of equation $(k)$ for $k=1,\dots ,6$. We reformulate the last four algebraic equations as differential ones:

$\begin{array}{rlrl}{\dot{y}}_7& =f_7=f_6+f_8+f_9+f_{10}& & (7^{\prime })\\ {\dot{y}}_8& =f_8=\frac{{\theta }_8{f}_1\cdot ({\theta }_8+y_7)-{\theta }_8{y}_1{f}_7}{({\theta }_8+y_7{)}^2}& & (8^{\prime })\\ {\dot{y}}_9& =f_9=\frac{{\theta }_9{f}_3\cdot ({\theta }_9+y_7)-{\theta }_9{y}_3{f}_7}{({\theta }_9+y_7{)}^2}& & (9^{\prime })\\ {\dot{y}}_{10}& =f_{10}=\frac{{\theta }_7{f}_5\cdot ({\theta }_7+y_7)-{\theta }_7{y}_5{f}_7}{({\theta }_7+y_7{)}^2}& & (10^{\prime })\end{array}$

The right hand sides of $(1)-(6)$ and $(7^{\prime })-(9^{\prime })$ are summarized as the vector-valued function $f$. Moreover, let

$\begin{array}{rl}{f}_{y,m,n}(\cdot )& =\frac{\partial f_m(\cdot )}{\partial y_n},m,n\in \{1,\dots ,10\}\\ f_{\theta ,m,n}(\cdot )& =\frac{\partial f_m(\cdot )}{\partial {\theta }_n},m\in \{1,\dots ,10\};n\in \{1,\dots ,9\}\end{array}$

Parameters

The initial parameter estimates are:

${\alpha }_1$	$2.0\cdot 10^{13}mol/(kg\cdot h)$
${\alpha }_2$	$2.0\cdot 10^{13}mol/(kg\cdot h)$
${\alpha }_{-1}$	$4.3\cdot 10^{15}mol/(kg\cdot h)$
$E_1$	$2.0\cdot 10^{4}cal/mol$
$E_2$	$2.0\cdot 10^{4}cal/mol$
$E_{-1}$	$2.0\cdot 10^{4}cal/mol$
$K_1$	$1.0\cdot 10^{-17}mol/kg$
$K_2$	$1.0\cdot 10^{-11}mol/kg$
$K_3$	$1.0\cdot 10^{-17}mol/kg$

Note that for the calculations all temperatures given in ${\circ }^{}C$ have to be rescaled to $K$ by adding $273.15$.

There are three datasets for different temperatures $T$, with corresponding starting values

	$40^{\circ }C$	$67^{\circ }C$	$100^{\circ }C$
$y_1(0)$	$1.7066$	$1.6749$	$1.5608$
$y_2(0)$	$8.32$	$8.2262$	$8.3546$
$y_3(0)$	$0.01$	$0.0104$	$0.0082$
$y_4(0)$	$0$	$0.0017$	$0.0086$

The initial model conditions in addition to those given in the data sets are:

$\begin{array}{l}{y}_5=0\\ y_6=[Q^{+}]\\ y_7=\frac{1}{2}\cdot (-K_2+\sqrt{K_2^2+4K_2{y}_1(0)})\\ y_8=y_7\\ y_9=0\\ y_{10}=0\end{array}$

To reduce the intercorrelation between the parameters in the rate constants, we apply the following reparametrization (cf. [4].):

$\begin{array}{rl}{k}_i& ={\alpha }_i\cdot \exp (-\frac{E_i}{RT})\\ & =k_{0,i}\cdot \exp (-\frac{E_i}{R}\cdot (\frac{1}{T}-\frac{1}{T_0})),i=1,2,-1\end{array}$

in which $k_{0,i}={\alpha }_i\cdot \exp (-\frac{E_i}{R{T}_0})$. The reference temperature in $T_0$ is chosen as the average over all performed experiments, i.e., $T_0=69^{\circ }C$. Additionally, we add a logarithmic transformation, which gives rise to the following transformed starting values:

$\ln k_{0,1}$	$1.194$
$\ln k_{0,2}$	$1.194$
$\ln k_{0,-1}$	$6.565$
$E_1$	$2.0\cdot 10^4$
$E_2$	$2.0\cdot 10^4$
$E_{-1}$	$2.0\cdot 10^4$
$\ln K_1$	$-34.54$
$\ln K_2$	$-25.33$
$\ln K_3$	$-39.14$

Optimal Experimental Design Problem

To be specified.

We are interested in when to measure (with an upper bound $M_i$ on the measuring time for each observable). We define

$\begin{array}{rlrl}{f}_y(\cdot )& \in {\mathbb{ℝ}}^{10\times 10}& & \text{ with }(f_y{)}_{i,j}=f_{y,i,j},\\ f_{\theta }(\cdot )& \in {\mathbb{ℝ}}^{10\times 9}& & \text{ with }(f_{\theta }{)}_{i,j}=f_{\theta ,i,j}\end{array}$

In this approach, we add the so-called sensitivities $G=dy/d\theta$. For the differential equations this means

$\dot{G}(t)=f_y(y(t),\theta )G(t)+f_{\theta }(y(t),\theta ),G(0)=\frac{\partial y(0)}{\partial \theta }$

Now we formulate the OED problem as described in [2].

$\begin{array}{lll}\underset{y,G,F,z,w}{\min }& & \text{trace}(F^{-1}(t_f))\\ \text{subject to}\\ \dot{y}(t)& =& f(y(t),\theta )\\ \dot{G}(t)& =& f_y(y(t),\theta )G(t)+f_{\theta }(y(t),\theta )\\ \dot{F}(t)& =& \sum _{i=1}^{n_o}{w}_i(t)(h_y^i(y(t))G(t){)}^T(h_y^i(y(t))G(t))\\ \dot{z}(t)& =& w(t),\\ y(0)& =& y_0\\ G(0)& =& \frac{\partial y(0)}{\partial \theta }\\ F(0)& =& 0,\\ z(0)& =& 0\\ w(t)& \in & {𝒲}\\ z_i(t_f)& \le & M_i\end{array}$

Here $h:{\mathbb{ℝ}}^{10}\to {\mathbb{ℝ}}^{n_o}$ is the observed function. The evolution of the symmetric matrix $F:[0,t_f]\to {\mathbb{ℝ}}^{9\times 9}$ is given by the weighted sum of observability Gramians $h_y^i(y(t))G(t),i=1,\dots ,n_o$ for each observed function of states. The weights $w_i(t)\in \{0,1\},i=1,\dots ,n_o$ are the (binary) sampling decisions, where $w_i(t)=1$ denotes the decision to perform a measurement at time $t$.

Miscellaneous and Further Reading

To be specified.

References

[1] "Nonlinear Parameter Estimation: a Case Study Comparison" by L. T. Biegler and J. J. Damiano
[2] "Optimal Experimental Design for Universal Differential Equations" by C. Plate, C.J. Martensen and S. Sager
[3] "Parameter estimation in nonlinear systems" by W.J.H. Stortelder
[4] "Parameter Estimation in Nonlinear Dynamical Systems" by Morten Rode Kristensen