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 "Nonlinear Parameter Estimation: a Case Study Comparison" by L. T. Biegler and J. J. Damiano add quote.

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_2\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}{c}\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 $\mathrm{gmol}/\mathrm{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}{c}{K}_1=\frac{[MB{M}^{-}][H^{+}]}{[(MBMH{)}_N]}\\ K_1=\frac{[A^{-}][H^{+}]}{[(HA{)}_N]}\\ K_1=\frac{[AB{M}^{-}][H^{+}]}{[(HABM{)}_N]}\end{array}$

The anionic species may then be represented by

$\begin{array}{cr}\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}{cr}\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}{cr}\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_3$ 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$ is the gas constant and $T$ is reaction temperature in Kelvins. The parameter ${\alpha }_i$, given in $\mathrm{mol}/(\mathrm{kg}\cdot \mathrm{h})$, represent the pre-exponential factor and $E_i$, with unit $\mathrm{cal}/\mathrm{mol}$, is the activation energy.

Mathematical formulation

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

$\begin{array}{lr}{\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]$ 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}{lr}{\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+f_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+f_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+f_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}{l}{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:

$\begin{array}{c}{\alpha }_1=2.0\times 10^{13}\\ E_1=2.0\times 10^4\\ {\alpha }_2=2.0\times 10^{13}\\ E_2=2.0\times 10^4\\ {\alpha }_{-1}=4.3\times 10^{15}\\ E_{-1}=2.0\times 10^4\\ K_1=1.0\times 10^{-17}\\ K_2=1.0\times 10^{-11}\\ K_3=1.0\times 10^{-17}\\ \end{array}$

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

$\begin{array}{lccc}& 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\end{array}$

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

$\begin{array}{l}{y}_5=0\\ y_6=0.0131\\ 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}$

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}{l}{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 (Optimal Experimental Design for Universal Differential Equations add quote)

$\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

To be added.