Difference between revisions of "DOW Experimental Design"

From mintOC
Jump to: navigation, search
(Parameters)
 
(13 intermediate revisions by the same user not shown)
Line 7: Line 7:
 
}}
 
}}
  
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 [[#Biegler1986]] "Nonlinear Parameter Estimation:
+
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 [[#Biegler1986|[1]]].
a Case Study Comparison" by L. T. Biegler and J. J. Damiano <span style="color:red">add quote</span>.
+
  
 
== Chemical background ==
 
== Chemical background ==
Line 41: Line 40:
 
<p>
 
<p>
 
<math>
 
<math>
  \begin{array}{c}  
+
  \begin{align}
   \left[HBMH\right] = \left[ (MBMH)_N \right] + \left[ MBM^- \right] \\
+
   \left[HBMH\right] &= \left[ (MBMH)_N \right] + \left[ MBM^- \right] \\
   \left[HA\right] = \left[ (HA)_N \right] + \left[ A^- \right] \\
+
   \left[HA\right] &= \left[ (HA)_N \right] + \left[ A^- \right] \\
   \left[HABM\right] = \left[ (HABM)_N \right] + \left[ ABM^- \right]
+
   \left[HABM\right] &= \left[ (HABM)_N \right] + \left[ ABM^- \right]
  \end{array}  
+
  \end{align}  
 
</math>
 
</math>
 
</p>
 
</p>
Line 52: Line 51:
 
<p>
 
<p>
 
<math>
 
<math>
\begin{array}{c}
+
\begin{align}
  K_1 = \frac{[MBM^-][H^+]}{[(MBMH)_N]} \\
+
  K_1 &= \frac{[MBM^-][H^+]}{[(MBMH)_N]} \\
  K_1 = \frac{[A^-][H^+]}{[(HA)_N]} \\
+
  K_1 &= \frac{[A^-][H^+]}{[(HA)_N]} \\
  K_1 = \frac{[ABM^-][H^+]}{[(HABM)_N]}  
+
  K_1 &= \frac{[ABM^-][H^+]}{[(HABM)_N]}  
\end{array}
+
\end{align}
 
</math>
 
</math>
 
</p>
 
</p>
Line 63: Line 62:
 
<p>
 
<p>
 
<math>
 
<math>
\begin{array}{cr}
+
\begin{align}
  \left[MBM^-\right] = \frac{K_1[MBMH]}{K_1 + [H^+]} & \quad (a) \\
+
  \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[A^-\right] &= \frac{K_2[HA]}{K_2 + [H^+]} && \quad (b)\\
  \left[ABM^-\right] = \frac{K_3[HABM]}{K_3 + [H^+]}& \quad (c)
+
  \left[ABM^-\right] &= \frac{K_3[HABM]}{K_3 + [H^+]} && \quad (c)
\end{array}
+
\end{align}
 
</math>
 
</math>
 
</p>
 
</p>
Line 74: Line 73:
 
<p>
 
<p>
 
<math>
 
<math>
\begin{array}{cr}
+
\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[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[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)  
+
  \frac{d[AB]}{dt} &= -k_3 \left[ M^- \right] \left[ AB \right] + k_{-3} \left[ ABM^- \right]  &&  \quad (f)  
\end{array}
+
\end{align}
 
</math>
 
</math>
 
</p>
 
</p>
Line 85: Line 84:
 
<p>
 
<p>
 
<math>
 
<math>
\begin{array}{cr}
+
\begin{align}
  \frac{d[MBMH]}{dt} = k_1 \left[ M^- \right] \left[ BM \right] + k_{-1} \left[ MBM^- \right]  &  \quad (g)\\
+
  \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[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)
+
  \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}
+
\end{align}
 
</math>
 
</math>
 
</p>
 
</p>
Line 118: Line 117:
 
<p>
 
<p>
 
  <math>
 
  <math>
   \begin{array}{lr}
+
   \begin{align}
   \dot{y}_1 = -k_2 y_8 y_2 & \quad (1),(h) \\
+
   \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}_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}_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}_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}_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) \\
+
   \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_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_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_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)
+
   y_{10} &= \frac{\theta_7 y_5}{\theta_7 + y_7} && \quad (10),(a)
   \end{array}
+
   \end{align}
 
  </math>
 
  </math>
 
</p>
 
</p>
Line 150: Line 149:
 
<p>
 
<p>
 
  <math>
 
  <math>
   \begin{array}{l r}
+
   \begin{align}
   \dot{y}_7 = f_7 = f_6 + f_8 + f_9 + f_{10} & \quad (7') \\
+
   \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}_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}_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')  
+
   \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}
+
   \end{align}
 
  </math>
 
  </math>
 
</p>
 
</p>
Line 162: Line 161:
 
<p>
 
<p>
 
  <math>
 
  <math>
   \begin{array}{l}
+
   \begin{align}
   f_{y,m,n}(\cdot) = \frac{\partial f_m(\cdot)}{\partial y_n}, \quad m,n \in \{1,\ldots,10\} \\
+
   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\}
+
   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}
+
   \end{align}
 
  </math>
 
  </math>
</p>  
+
</p>
  
 
== Parameters ==
 
== Parameters ==
  
 
The initial parameter estimates are:
 
The initial parameter estimates are:
<table border="1.5">
+
<table style="border-collapse: collapse;" border="1.">
 
             <tr>
 
             <tr>
 
                 <td style="text-align: center; padding:5pt"> <math> \alpha_1 </math> </td>
 
                 <td style="text-align: center; padding:5pt"> <math> \alpha_1 </math> </td>
Line 210: Line 209:
 
             </tr>
 
             </tr>
 
     </table>
 
     </table>
 +
  
 
There are three datasets for different temperatures <math>T</math>, with corresponding starting values
 
There are three datasets for different temperatures <math>T</math>, with corresponding starting values
 +
<table style="border-collapse: collapse;" border="1.">
 +
            <tr style="border-bottom: 2pt solid black;">
 +
                <th style="border-right: 2pt solid black;"></th>
 +
                <th style="text-align: center; padding:5pt"> <math>40^\circ C </math></th>
 +
                <th style="text-align: center; padding:5pt"> <math>67^\circ C </math></th>
 +
                <th style="text-align: center; padding:5pt"> <math>100^\circ C </math></th>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt; border-right: 2pt solid black;"> <math> y_1(0) </math> </td>
 +
                <td style="text-align: center; padding:5pt"> <math> 1.7066 </math> </td>
 +
                <td style="text-align: center; padding:5pt"> <math> 1.6749 </math> </td>
 +
                <td style="text-align: center; padding:5pt"> <math> 1.5608 </math> </td>           
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt; border-right: 2pt solid black"><math>y_2(0)</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>8.32</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>8.2262</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>8.3546</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt; border-right: 2pt solid black"><math>y_3(0)</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>0.01</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>0.0104</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>0.0082</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt; border-right: 2pt solid black"><math>y_4(0)</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>0</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>0.0017</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>0.0086</math></td>
 +
            </tr>
 +
    </table>
 +
<!--
 
<p>
 
<p>
 
  <math>
 
  <math>
Line 223: Line 256:
 
  </math>
 
  </math>
 
</p>
 
</p>
 +
-->
 
The initial model conditions in addition to those given in the
 
The initial model conditions in addition to those given in the
 
data sets are:
 
data sets are:
Line 245: Line 279:
 
<p>
 
<p>
 
  <math>
 
  <math>
   \begin{array}{c}
+
   \begin{align}
   k_{i,\alpha} := \frac{\partial k_i}{\partial \alpha_i} = \exp(-E_i/(RT)) \\
+
   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))
+
   k_{i,E} &:= \frac{\partial k_i}{\partial E_i} = -\frac{\alpha_i}{RT} \exp(-E_i/(RT))
   \end{array}
+
   \end{align}
 
  </math>
 
  </math>
 
</p>
 
</p>
Line 285: Line 319:
 
<p>
 
<p>
 
  <math>
 
  <math>
\begin{array}{l}
+
\begin{align}
   f_y(\cdot) \in \mathbb{R}^{10 \times 10} \quad \text{ with } (f_y)_{i,j} = f_{y,i,j}, \\
+
   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}
+
   f_\theta(\cdot) &\in \mathbb{R}^{10 \times 9} \quad &&\text{ with } (f_\theta)_{i,j} = f_{\theta,i,j}
  \end{array}
+
  \end{align}
 
  </math>
 
  </math>
 
</p>  
 
</p>  
Line 298: Line 332:
 
</p>
 
</p>
  
Now we formulate the OED problem as described in (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>
Line 328: Line 362:
  
 
== References ==
 
== References ==
To be added. <br>
+
<span id="Biegler1986">[1]</span> "Nonlinear Parameter Estimation: a Case Study Comparison" by L. T. Biegler and J. J. Damiano <br>
<!--
+
<span id="OEDUDE">[2]</span> "Optimal Experimental Design for Universal Differential Equations" by C. Plate, C.J. Martensen and S. Sager <br>
<biblist />
+
<span id="Stortelder1998">[3]</span> "Parameter estimation in nonlinear systems" by W.J.H. Stortelder <br>
-->
+
<span id="Biegler1986">[1]</span> "Nonlinear Parameter Estimation: a Case Study Comparison" by L. T. Biegler and J. J. Damiano  
+
  
 
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here -->
 
<!--List of all categories this page is part of. List characterization of solution behavior, model properties, ore presence of implementation details (e.g., AMPL for AMPL model) here -->

Latest revision as of 09:56, 13 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{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_1 &= \frac{[A^-][H^+]}{[(HA)_N]} \\
 K_1 &= \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 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] 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 + 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{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