Difference between revisions of "DOW Experimental Design"

From mintOC
Jump to: navigation, search
(→‎Optimal Experimental Design Problem)
(→‎Parameters)
 
(53 intermediate revisions by 2 users 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 "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 31: Line 30:
 
  MBMH \leftarrow K_1 \rightarrow MBM^- + H^+ \\
 
  MBMH \leftarrow K_1 \rightarrow MBM^- + H^+ \\
 
  HA \leftarrow K_2 \rightarrow A^- + H^+ \\
 
  HA \leftarrow K_2 \rightarrow A^- + H^+ \\
−
  HABM \leftarrow K_2 \rightarrow ABM^- + H^+  
+
  HABM \leftarrow K_3 \rightarrow ABM^- + H^+  
 
  \end{array}
 
  \end{array}
 
</math>
 
</math>
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>
−
Here <math>[ \ ]</math> denotes the concentration of the species in <math>\operatorname{gmol}/\operatorname{kg}</math>.
+
Here <math>[ \ ]</math> denotes the concentration of the species in <math>mol/kg</math>.
 
By assuming the rapid acid-base reactions are at equilibrium, the equilibrium constants <math>K_1,K_2,K_3</math> can be defined as
 
By assuming the rapid acid-base reactions are at equilibrium, the equilibrium constants <math>K_1,K_2,K_3</math> can be defined as
 
<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_2 &= \frac{[A^-][H^+]}{[(HA)_N]} \\
−
  K_1 = \frac{[ABM^-][H^+]}{[(HABM)_N]}  
+
  K_3 &= \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>
−
An electroneutrality constraint gives the hydrogen ion con-
+
An electroneutrality constraint gives the hydrogen ion concentration <math>\left[ H^+ \right]</math> as
−
centration <math>\left[ H^+ \right]</math> as
+
 
<p>
 
<p>
 
  <math>
 
  <math>
Line 104: Line 102:
 
  <math> k_3 = k_1, \quad k_{-3} = \frac{1}{2} k_{-1} </math>
 
  <math> k_3 = k_1, \quad k_{-3} = \frac{1}{2} k_{-1} </math>
 
</p>
 
</p>
−
With these assumptions, the three rate constants <math>k_1,k_2</math> and <math>k_3</math> must be estimated. Each of these can be fitted with two adjustable model parameters, assuming an Arrhenius temperature dependence. That is
+
With these assumptions, the three rate constants <math>k_1,k_2</math> and <math>k_{-1}</math> must be estimated. Each of these can be fitted with two adjustable model parameters, assuming an Arrhenius temperature dependence. That is
 
<p>
 
<p>
 
  <math>  
 
  <math>  
Line 112: Line 110:
 
  </math>
 
  </math>
 
</p>
 
</p>
−
Here <math>R</math> is the gas constant and <math>T</math> is reaction temperature in Kelvins. The parameter <math>\alpha_i</math>, given in <math>\operatorname{mol}/( \operatorname{kg} \cdot \operatorname{h})</math>, represent the pre-exponential factor and <math>E_i</math>, with unit <math>\operatorname{cal}/{\operatorname{mol}}</math>, is the activation energy.  
+
Here <math>R \approx 1.98720425864083 \ cal/(K \cdot mol) </math> is the gas constant and <math>T</math> is the reaction temperature in Kelvins. The parameters <math>\alpha_i</math>, given in <math>mol/( kg \cdot h)</math>, represent the pre-exponential factors and the <math>E_i</math>, with unit <math>cal/mol</math>, are the activation energies.
  
 
== Mathematical formulation ==
 
== Mathematical formulation ==
Line 119: 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>
−
Here the letters in parentheses stand for the corresponding chemical process and the quantity <math>\left[Q^+\right]</math> is a constant during the reaction.
+
Here the letters in parentheses stand for the corresponding chemical process and the quantity <math>\left[Q^+\right]=0.0131</math> is a constant during the reaction.
 
The nine parameters form the vector  
 
The nine parameters form the vector  
 
<p>
 
<p>
Line 151: 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 + 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 + 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 + 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 + 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 + y_7) - \theta_7 y_5 f_7}{(\theta_7 + y_7)^2} && \quad (10')  
−
   \end{array}
+
   \end{align}
 
  </math>
 
  </math>
 
</p>
 
</p>
Line 163: 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:
−
<p>
+
<table style="border-collapse: collapse;" border="1.">
−
<math>
+
            <tr>
−
  \begin{array}{c}
+
                <td style="text-align: center; padding:5pt"> <math> \alpha_1 </math> </td>
−
  \alpha_1 = 2.0 \times 10^{13} \\
+
                <td style="text-align: center; padding:5pt"> <math>2.0 \cdot 10^{13} \ mol / (kg \cdot h)</math> </td>
−
  E_1 = 2.0 \times 10^{4} \\
+
            </tr>
−
  \alpha_2 = 2.0 \times 10^{13} \\
+
            <tr>
−
  E_2 = 2.0 \times 10^{4} \\
+
                <td style="text-align: center; padding:5pt"><math>\alpha_2</math></td>
−
  \alpha_{-1} = 4.3 \times 10^{15} \\
+
                <td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{13} \ mol / (kg \cdot h)</math></td>
−
  E_{-1} = 2.0 \times 10^{4} \\
+
            </tr>
−
  K_1 = 1.0 \times 10^{-17} \\
+
            <tr>
−
  K_2 = 1.0 \times 10^{-11} \\
+
                <td style="text-align: center; padding:5pt"><math>\alpha_{-1}</math></td>
−
  K_3 = 1.0 \times 10^{-17} \\
+
                <td style="text-align: center; padding:5pt"><math>4.3 \cdot 10^{15} \ mol / (kg \cdot h)</math></td>
−
  \end{array}
+
            </tr>
−
</math>
+
            <tr>
−
</p>
+
                <td style="text-align: center; padding:5pt"><math>E_1</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4} \ cal / mol</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>E_2</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4} \ cal / mol </math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>E_{-1}</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4} \ cal / mol</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>K_1</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>1.0\cdot 10^{-17} \ mol / kg</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>K_2</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>1.0\cdot 10^{-11} \ mol / kg</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>K_3</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>1.0\cdot 10^{-17} \ mol / kg</math></td>
 +
            </tr>
 +
    </table>
 +
Note that for the calculations all temperatures given in <math>^{\circ}C</math> have to be rescaled to <math>K</math> by adding <math>273.15</math>.
 +
 
 
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 200: 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 206: Line 263:
 
   \begin{array}{l}
 
   \begin{array}{l}
 
   y_5 = 0 \\
 
   y_5 = 0 \\
−
   y_6 = 0.0131 \\
+
   y_6 = [Q^+] \\
 
   y_7 = \frac{1}{2} \cdot \left( -K_2 + \sqrt{K_2^2 + 4K_2 y_1(0)} \right) \\
 
   y_7 = \frac{1}{2} \cdot \left( -K_2 + \sqrt{K_2^2 + 4K_2 y_1(0)} \right) \\
 
   y_8 = y_7 \\
 
   y_8 = y_7 \\
Line 214: Line 271:
 
  </math>
 
  </math>
 
</p>
 
</p>
 +
 +
To reduce the intercorrelation between the parameters in the rate constants, we apply the following reparametrization (cf. [[#Kristensen2004|[4]]].):
 +
<p>
 +
<math>
 +
  \begin{align}
 +
  k_i &= \alpha_i \cdot \exp \left( - \frac{E_i}{RT} \right) \\
 +
      &= k_{0,i} \cdot \exp \left( - \frac{E_i}{R} \cdot \left( \frac{1}{T} - \frac{1}{T_0} \right) \right), \quad i=1,2,-1
 +
\end{align}
 +
</math>
 +
</p>
 +
in which <math>k_{0,i} = \alpha_i \cdot \exp\left( - \frac{E_i}{RT_0} \right)</math>. The reference temperature in <math>T_0</math> is chosen as the average over all performed experiments, i.e., <math>T_0 = 69^\circ C</math>. Additionally, we add a logarithmic transformation, which gives rise to the following transformed starting values:
 +
 +
<table style="border-collapse: collapse;" border="1.">
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"> <math> \ln k_{0,1} </math> </td>
 +
                <td style="text-align: center; padding:5pt"> <math>1.194</math> </td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>\ln k_{0,2}</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>1.194</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>\ln k_{0,-1}</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>6.565</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>E_1</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4}</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>E_2</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4} </math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>E_{-1}</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>2.0 \cdot 10^{4}</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>\ln K_1</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>-34.54</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>\ln K_2</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>-25.33</math></td>
 +
            </tr>
 +
            <tr>
 +
                <td style="text-align: center; padding:5pt"><math>\ln K_3</math></td>
 +
                <td style="text-align: center; padding:5pt"><math>-39.14</math></td>
 +
            </tr>
 +
    </table>
  
 
== Optimal Experimental Design Problem ==
 
== Optimal Experimental Design Problem ==
Line 222: Line 329:
 
<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 262: Line 369:
 
<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_i_{\theta_j,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 275: Line 382:
 
</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 305: Line 412:
  
 
== References ==
 
== References ==
−
To be added.
+
<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="Kristensen2004">[4]</span> "Parameter Estimation in Nonlinear Dynamical Systems" by Morten Rode Kristensen <br>
  
 
<!--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 11:31, 18 February 2025

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 mol/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_{-1} 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 1.98720425864083 \ cal/(K \cdot mol) is the gas constant and T is the reaction temperature in Kelvins. The parameters \alpha_i, given in mol/( kg \cdot h), represent the pre-exponential factors and the E_i, with unit cal/mol, are the activation energies.

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.0131 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 \cdot 10^{13} \ mol / (kg \cdot h)
\alpha_2 2.0 \cdot 10^{13} \ mol / (kg \cdot h)
\alpha_{-1} 4.3 \cdot 10^{15} \ mol / (kg \cdot h)
E_1 2.0 \cdot 10^{4} \ cal / mol
E_2 2.0 \cdot 10^{4} \ cal / mol
E_{-1} 2.0 \cdot 10^{4} \ cal / mol
K_1 1.0\cdot 10^{-17} \ mol / kg
K_2 1.0\cdot 10^{-11} \ mol / kg
K_3 1.0\cdot 10^{-17} \ mol / kg

Note that for the calculations all temperatures given in ^{\circ}C have to be rescaled to K by adding 273.15.

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 = [Q^+] \\ 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}

To reduce the intercorrelation between the parameters in the rate constants, we apply the following reparametrization (cf. [4].):

\begin{align} k_i &= \alpha_i \cdot \exp \left( - \frac{E_i}{RT} \right) \\ &= k_{0,i} \cdot \exp \left( - \frac{E_i}{R} \cdot \left( \frac{1}{T} - \frac{1}{T_0} \right) \right), \quad i=1,2,-1 \end{align}

in which k_{0,i} = \alpha_i \cdot \exp\left( - \frac{E_i}{RT_0} \right). The reference temperature in T_0 is chosen as the average over all performed experiments, i.e., T_0 = 69^\circ C. Additionally, we add a logarithmic transformation, which gives rise to the following transformed starting values:

\ln k_{0,1} 1.194
\ln k_{0,2} 1.194
\ln k_{0,-1} 6.565
E_1 2.0 \cdot 10^{4}
E_2 2.0 \cdot 10^{4}
E_{-1} 2.0 \cdot 10^{4}
\ln K_1 -34.54
\ln K_2 -25.33
\ln K_3 -39.14

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
[4] "Parameter Estimation in Nonlinear Dynamical Systems" by Morten Rode Kristensen

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