Difference between revisions of "DOW Experimental Design"

From mintOC
Jump to: navigation, search
(→‎References)
(→‎Optimal Experimental Design Problem)
Line 297: Line 297:
 
</p>
 
</p>
  
−
Now we formulate the OED problem as described in [[#OEDUDE | [2]]] (Optimal Experimental Design for Universal Differential Equations <span style="color:red">add quote</span>)
+
Now we formulate the OED problem as described in [[#OEDUDE | [2]]].
 
<p>
 
<p>
 
<math>
 
<math>

Revision as of 13:26, 12 November 2024

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_2 \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{array}{c} \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{array}

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{array}{c} K_1 = \frac{[MBM^-][H^+]}{[(MBMH)_N]} \\ K_1 = \frac{[A^-][H^+]}{[(HA)_N]} \\ K_1 = \frac{[ABM^-][H^+]}{[(HABM)_N]} \end{array}

The anionic species may then be represented by

\begin{array}{cr} \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{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[ 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{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[ 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{array}

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 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{array}{lr} \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{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^-,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{array}{l r} \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 + f_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 + f_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 + f_7) - \theta_7 y_5 f_7}{(\theta_7 + y_7)^2} & \quad (10') \end{array}

The right hand sides of (1)-(6) and (7')-(9') 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}, \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{array}

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

\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 \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{array}{l} 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{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), \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

Retrieved from "https://mintoc.de/index.php?title=DOW_Experimental_Design&oldid=2418"