Difference between revisions of "DOW Experimental Design"

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

Revision as of 13:45, 23 September 2024

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.

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 con- centrations 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 $gmol/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 con- centration $\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, $T$ is reaction temperature in Kelvins and the parameters $\alpha_i,E_i$ represent the pre-exponential factor and activation energy, respectively, for the appropriate rate constant.

Mathematical formulation

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

$\begin{array}{lr} \frac{d y_1}{dt} = -k_2 y_8 y_2 & \quad (1),(h) \\ \frac{d y_2}{dt} = -k_1 y_6 y_2 + k_{-1} y_{10} - k_2 y_8 y_2 & \quad (2),(e) \\ \frac{d y_3}{dt} = -k_2 y_8 y_2 + k_1 y_6 y_4 - \frac{1}{2} k_{-1} y_9 & \quad (3),(i) \\ \frac{d y_4}{dt} = -k_1 y_6 y_4 + \frac{1}{2} k_{-1} y_9 & \quad (4),(f) \\ \frac{d y_5}{dt} = -k_1 y_6 y_2 + k_{-1} y_{10} & \quad (5),(g) \\ \frac{d y_6}{dt} = -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 letter stands for the corresponding chemical process. 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} \frac{d y_7}{dt} = f_7 = f_6 + f_8 + f_9 + f_{10} & \quad (7') \\ \frac{d y_8}{dt} = 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') \\ \frac{d y_9}{dt} = 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')\\ \frac{d y_{10}}{dt} = 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,6\} \\ f_{\theta,m,n}(\cdot) = \frac{\partial f_m(\cdot)}{\partial y_n} \end{array}$

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} = \frac{\partial f_i}{\partial \theta_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 (Optimal Experimental Design for Universal Differential Equations add quote)

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

Parameters

We use $t_f=12$ , $p_1 = p_2 = p_3 = p_4 = 1$ , and $p_5 = 0.4$ , $p_6 = 0.2$ . The upper bound on the measurement time intervals is chosen as $M=4$ .

Miscellaneous and Further Reading

The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper [Sager2006]Address: Heidelberg
Author: S. Sager; H.G. Bock; M. Diehl; G. Reinelt; J.P. Schl\"oder
Booktitle: Recent Advances in Optimization
Editor: A. Seeger
Note: ISBN 978-3-5402-8257-0
Pages: 269--289
Publisher: Springer
Series: Lectures Notes in Economics and Mathematical Systems
Title: Numerical methods for optimal control with binary control functions applied to a Lotka-Volterra type fishing problem
Volume: 563
Year: 2009
and revisited in his PhD thesis [Sager2005]Address: Tönning, Lübeck, Marburg
Author: S. Sager
Editor: ISBN 3-89959-416-9
Publisher: Der andere Verlag
Title: Numerical methods for mixed--integer optimal control problems
Url: http://mathopt.de/PUBLICATIONS/Sager2005.pdf
Year: 2005
. These are also the references to look for more details. The experimental design problem has been described in the habilitation thesis of Sager, [Sager2011d]Author: S. Sager
How published: University of Heidelberg
Month: August
Note: Habilitation
Title: On the Integration of Optimization Approaches for Mixed-Integer Nonlinear Optimal Control
Url: http://mathopt.de/PUBLICATIONS/Sager2011d.pdf
Year: 2011
.

References

[Sager2005]	S. Sager (2005): Numerical methods for mixed--integer optimal control problems. (%edition%). Der andere Verlag, Tönning, Lübeck, Marburg, %pages%
[Sager2006]	S. Sager; H.G. Bock; M. Diehl; G. Reinelt; J.P. Schl\"oder (2009): Numerical methods for optimal control with binary control functions applied to a Lotka-Volterra type fishing problem. Springer, Recent Advances in Optimization
[Sager2011d]	S. Sager: On the Integration of Optimization Approaches for Mixed-Integer Nonlinear Optimal Control, 2011

@@ Line 146: / Line 146: @@
 </p>
-Let <math>f_k(\cdot)</math> denote the right hand side of equation <math>(k)</math> for <math>k=1,\ldots,6</math>. Moreover, let
+Let <math>f_k(\cdot)</math> denote the right hand side of equation <math>(k)</math> for <math>k=1,\ldots,6</math>.
+We reformulate the last four algebraic equations as differential ones:
+<p>
+ <math>
+  \begin{array}{l r}
+   \frac{d y_7}{dt} = f_7 = f_6 + f_8 + f_9 + f_{10} & \quad (7') \\
+   \frac{d y_8}{dt} = 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') \\
+   \frac{d y_9}{dt} = 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')\\
+   \frac{d y_{10}}{dt} = 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}
+ </math>
+</p>
+The right hand sides of <math>(1)-(6)</math> and <math>(7')-(9')</math> are summarized as the vector-valued function <math>f</math>.
+Moreover, let
 <p>
   <math>
    \begin{array}{l}
-   f_{ymn}(\cdot) = \frac{\partial f_m(\cdot)}{\partial y_n}, \quad m,n \in \{1,\ldots,6\} \\
+   f_{y,m,n}(\cdot) = \frac{\partial f_m(\cdot)}{\partial y_n}, \quad m,n \in \{1,\ldots,6\} \\
-   f_{\theta mn}(\cdot) = \frac{\partial f_m(\cdot)}{\partial y_n}
+   f_{\theta,m,n}(\cdot) = \frac{\partial f_m(\cdot)}{\partial y_n}
    \end{array}
   </math>
 </p>
-The non-zero derivatives are given by
+<!--
+To simplify the notation, the derivatives of the <math>k_i</math> with respect to the <math>\alpha_i</math> and <math>E_i</math> are written as
+<p>
+ <math>
+  \begin{array}{c}
+  k_{i,\alpha} := \frac{\partial k_i}{\partial \alpha_i} = \exp(-E_i/(RT)) \\
+  k_{i,E} := \frac{\partial k_i}{\partial E_i} = -\frac{\alpha_i}{RT} \exp(-E_i/(RT))
+  \end{array}
+ </math>
+</p>
+The non-zero derivatives for the first six components are given by
 <p>
 <math>
-  \begin{array}{l}
+  \begin{array}{lcl}
-  f_{y12} = -k_2 y_8 \\
+  f_{y,1,2} = -k_2 y_8 & & \\
-  f_{y18} = -k_2 y_2 \\
+  f_{y,1,8} = -k_2 y_2 & & \\
-  f_{y22} = -k_1 y_6 - k_2 y_8 \\
+  f_{y,2,2} = -k_1 y_6 - k_2 y_8 & & \\
-  f_{y26} = -k_1 y_2 \\
+  f_{y,2,6} = -k_1 y_2 & & \\
-  f_{y28} = -k_2 y_2 \\
+  f_{y,2,8} = -k_2 y_2 & & \\
-  f_{y32} = -k_2 y_8 \\
+  f_{y,2,10} = k_{-1} & & \\
-  f_{y34} = k_1 y_6 \\
+ f_{y,3,2} = -k_2 y_8 & & \\
-  f_{y36} = k_1 y_4 \\
+  f_{y,3,4} = k_1 y_6 & & \\
-  f_{y44} = -k_1 y_6 \\
+  f_{y,3,6} = k_1 y_4 & & \\
-  f_{y46} = -k_1 y_4 \\
+  f_{y,3,8} = -k_2 y_2 & & \\
-  f_{52} = -k_1 y_6 \\
+ f_{y,3,9} = -\frac{1}{2} k_{-1} & & \\
-  f_{y56} = -k_1 y_2 \\
+ f_{y,4,4} = -k_1 y_6 & & \\
-  f_{y62} = -k_1 y_6 \\
+  f_{y,4,6} = -k_1 y_4 & & \\
-  f_{y64} = -k_1 y_6 \\
+  f_{y,4,9} = \frac{1}{2} k_{-1} & & \\
-  f_{y66} = -k_1 (y_2 + y_4)
+ f_{y,5,2} = -k_1 y_6 & & \\
+  f_{y,5,6} = -k_1 y_2 & & \\
+  f_{y,5,10} = k_{-1} & & \\
+ f_{y,6,2} = -k_1 y_6 & & \\
+  f_{y,6,4} = -k_1 y_6 & & \\
+  f_{y,6,6} = -k_1 (y_2 + y_4) & & \\
+ f_{y,6,9} = \frac{1}{2} k_{-1} & & \\
+ f_{y,6,10} = k_{-1} & &
   \end{array}
 </math>
+</p>
+-->
+We are interested in when to measure (with an upper bound <math>M_i</math> on the measuring time for each observable).
+We define
+<p>
+ <math>
+\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} = \frac{\partial f_i}{\partial \theta_j}
+ \end{array}
+ </math>
+</p>
+In this approach, we add the so-called sensitivities <math>G=dy/d\theta</math>. For the differential equations this means
+<p>
+ <math>
+  \dot{G}(t) = f_y(y(t),\theta) G(t) + f_\theta(y(t),\theta), \quad G(0) = \frac{\partial y(0)}{\partial \theta}
+ </math>
 </p>
+Now we formulate the OED problem as described in (Optimal Experimental Design for Universal Differential Equations <span style="color:red">add quote</span>)
+<p>
+<math>
+ \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}
+</math>
+</p>
-We are interested in when to measure, with an upper bound <math>M</math> on the measuring time. We can measure the states directly, <math>h^1(x(t)) = x_1(t)</math> and <math>h^2(x(t)) = x_2(t)</math>. We use two different sampling functions, <math>w^1(\cdot)</math> and <math>w^2(\cdot)</math> in the same experimental setting. This can be seen either as a two-dimensional measurement function <math>h(x(t))</math>, or as a special case of a multiple experiment, in which <math>u(\cdot), x(\cdot)</math>, and <math>G(\cdot)</math> are identical. The experimental design problem then reads
+Here <math>h:\mathbb{R}^{10} \rightarrow \mathbb{R}^{n_o}</math> is the observed function. The
+evolution of the symmetric matrix <math>F: \left[0,t_f \right] \rightarrow \mathbb{R}^{9 \times 9}</math> is given by the weighted sum of observability Gramians
+<math>h^i_y (y(t)) G(t), \ i = 1,\ldots,n_o</math> for each observed function of states. The weights <math>w_i (t) \in \{0, 1\}, \ i = 1,\ldots,n_o</math>
+are the (binary) sampling decisions, where <math> w_i (t) = 1 </math> denotes the decision to perform a measurement at
+time <math>t</math>.
 == Parameters ==

Difference between revisions of "DOW Experimental Design"

Revision as of 13:45, 23 September 2024

Contents

Mathematical formulation

Parameters

Miscellaneous and Further Reading

References

Navigation menu

Personal tools

Namespaces

Variants

Views

More

Search

Navigation

Support

Tools