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 \rightarrow 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} \overset{k_1}{\rightarrow} MBM^- \\ A^- + BM \overset{k_2}{\rightarrow} ABM^- \\ M^- + AB \underset{k_{-3}}{\leftarrow} \overset{k_3}{\rightarrow} ABM^- \end{array}$

Acid-Base Reactions:

$\begin{array}{c} MBMH \leftarrow K_1 \rightarrow MBM^- + H^+ \\ HA \leftarrow K_2 \rightarrow A^- + H^+ \\ HABM \leftarrow K_3 \rightarrow ABM^- + 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{align} \left[HBMH\right] &= \left[ (MBMH)_N \right] + \left[ MBM^- \right] \\ \left[HA\right] &= \left[ (HA)_N \right] + \left[ A^- \right] \\ \left[HABM\right] &= \left[ (HABM)_N \right] + \left[ ABM^- \right] \end{align}$

Here $[ \ ]$ denotes the concentration of the species in $\operatorname{gmol}/\operatorname{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{align} K_1 &= \frac{[MBM^-][H^+]}{[(MBMH)_N]} \\ K_2 &= \frac{[A^-][H^+]}{[(HA)_N]} \\ K_3 &= \frac{[ABM^-][H^+]}{[(HABM)_N]} \end{align}$

The anionic species may then be represented by

$\begin{align} \left[MBM^-\right] &= \frac{K_1[MBMH]}{K_1 + [H^+]} && \quad (a) \\ \left[A^-\right] &= \frac{K_2[HA]}{K_2 + [H^+]} && \quad (b)\\ \left[ABM^-\right] &= \frac{K_3[HABM]}{K_3 + [H^+]} && \quad (c) \end{align}$

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

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

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

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

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

$\left[ H^+ \right] + \left[ Q^+ \right] = \left[ M^- \right] + \left[ MBM^- \right] + \left[ A^- \right] + \left[ ABM^- \right] \quad \quad (j)$

Based on similarities of reacting species, we assume

$k_3 = k_1, \quad 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)), \quad i \in \{-1,1,2\} \end{array}$

Here $R \approx 8.314462 \ J/(\operatorname{mol} K)$ is the gas constant and $T$ is reaction temperature in Kelvins. The parameter $\alpha_i$, given in $\operatorname{mol}/( \operatorname{kg} \cdot \operatorname{h})$, represent the pre-exponential factor and $E_i$, with unit $\operatorname{cal}/{\operatorname{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{align} \dot{y}_1 &= -k_2 y_8 y_2 && \quad (1),(h) \\ \dot{y}_2 &= -k_1 y_6 y_2 + k_{-1} y_{10} - k_2 y_8 y_2 && \quad (2),(e) \\ \dot{y}_3 &= -k_2 y_8 y_2 + k_1 y_6 y_4 - \frac{1}{2} k_{-1} y_9 && \quad (3),(i) \\ \dot{y}_4 &= -k_1 y_6 y_4 + \frac{1}{2} k_{-1} y_9 && \quad (4),(f) \\ \dot{y}_5 &= -k_1 y_6 y_2 + k_{-1} y_{10} && \quad (5),(g) \\ \dot{y}_6 &= -k_1 (y_6 y_2 + y_6 y_4) + k_{-1} (y_{10} + \frac{1}{2} y_9) && \quad (6),(d) \\ y_7 &= -\left[ Q^+ \right] + y_6 + y_8 + y_9 + y_{10} && \quad (7),(j)\\ y_8 &= \frac{\theta_8 y_1}{\theta_8 + y_7} && \quad (8),(b)\\ y_9 &= \frac{\theta_9 y_3}{\theta_9 + y_7} && \quad (9),(c)\\ y_{10} &= \frac{\theta_7 y_5}{\theta_7 + y_7} && \quad (10),(a) \end{align}$

Here the letters in parentheses stand for the corresponding chemical process and the quantity $\left[Q^+\right]=0.00131$ 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^-,ABM^-,MBM^-)$

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

$\begin{align} \dot{y}_7 &= f_7 = f_6 + f_8 + f_9 + f_{10} && \quad (7') \\ \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} && \quad (8') \\ \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} && \quad (9')\\ \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} && \quad (10') \end{align}$

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

$\begin{align} f_{y,m,n}(\cdot) &= \frac{\partial f_m(\cdot)}{\partial y_n}, \quad m,n \in \{1,\ldots,10\} \\ f_{\theta,m,n}(\cdot) &= \frac{\partial f_m(\cdot)}{\partial \theta_n}, \quad m \in \{1,\ldots,10\}; \ n\in \{1,\ldots,9\} \end{align}$

Parameters

The initial parameter estimates are:

$\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\cdot 10^{-17}$
$K_2$	$1.0\cdot 10^{-11}$
$K_3$	$1.0\cdot 10^{-17}$

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 = 0.0131 \\ y_7 = \frac{1}{2} \cdot \left( -K_2 + \sqrt{K_2^2 + 4K_2 y_1(0)} \right) \\ 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{align} f_y(\cdot) &\in \mathbb{R}^{10 \times 10} \quad &&\text{ with } (f_y)_{i,j} = f_{y,i,j}, \\ f_\theta(\cdot) &\in \mathbb{R}^{10 \times 9} \quad &&\text{ with } (f_\theta)_{i,j} = f_{\theta,i,j} \end{align}$

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), \quad G(0) = \frac{\partial y(0)}{\partial \theta}$

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

$\begin{array}{lll} \displaystyle \min_{y,G,F,z,w} && \text{trace} \; \left( F^{-1}(t_f) \right) \\ \text{subject to} \\ \quad \dot{y}(t) & = & f(y(t),\theta) \\ \quad \dot{G}(t) & = & f_y(y(t),\theta) G(t) + f_\theta(y(t),\theta) \\ \quad \dot{F}(t) & = & \sum_{i=1}^{n_o} w_i(t)(h^i_y(y(t))G(t))^T(h^i_y(y(t))G(t)) \\ \quad \dot{z}(t) & = & w(t), \\ \quad y(0) & = & y_0 \\ \quad G(0) & = & \frac{\partial y(0)}{\partial \theta} \\ \quad F(0) & = & 0, \\ \quad z(0) & = & 0 \\ \quad w(t) & \in & \mathcal{W} \\ \quad z_i(t_f) & \leq & M_i \end{array}$

Here $h:\mathbb{R}^{10} \rightarrow \mathbb{R}^{n_o}$ is the observed function. The evolution of the symmetric matrix $F: \left[0,t_f \right] \rightarrow \mathbb{R}^{9 \times 9}$ is given by the weighted sum of observability Gramians $h^i_y (y(t)) G(t), \ i = 1,\ldots,n_o$ for each observed function of states. The weights $w_i (t) \in \{0, 1\}, \ i = 1,\ldots,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