Pyomo's parameter estimation in an ODE system with missing values in time series
<p>I have an ODE system of 7 equations for explaining a particular set of microorganisms dynamics of the form:</p>

<p><img src="http://latex.codecogs.com/gif.latex?%5Cbegin%7Balign*%7D&space;%5Cfrac%7BdX_1%7D%7Bdt%7D&amp;=-Y_1%5C,.%5C,v_1-Y_2%5C,.%5C,v_3%5C%5C[4pt]&space;%5Cfrac%7BdX_2%7D%7Bdt%7D&amp;=v_1-v_2%5C%5C[4pt]&space;%5Cfrac%7BdX_3%7D%7Bdt%7D&amp;=Y_3%5C,.%5C,v_1-Y_4%5C,.%5C,v_5%5C%5C[4pt]&space;%5Cfrac%7BdX_4%7D%7Bdt%7D&amp;=v_3-v_4%5C%5C[4pt]&space;%5Cfrac%7BdX_5%7D%7Bdt%7D&space;&amp;=Y_5%5C,.%5C,v_3%5C%5C[4pt]&space;%5Cfrac%7BdX_6%7D%7Bdt%7D&space;&amp;=&space;v_3-v_6%5C%5C[4pt]&space;%5Cfrac%7BdX_7%7D%7Bdt%7D&amp;=Y_6%5C,.%5C,v_5%5C%5C[4pt]&space;%5Cend%7Balign*%7D" title="\begin{align*} \frac{dX_1}{dt}&amp;=-Y_1\,.\,v_1-Y_2\,.\,v_3\\[4pt] \frac{dX_2}{dt}&amp;=v_1-v_2\\[4pt] \frac{dX_3}{dt}&amp;=Y_3\,.\,v_1-Y_4\,.\,v_5\\[4pt] \frac{dX_4}{dt}&amp;=v_3-v_4\\[4pt] \frac{dX_5}{dt} &amp;=Y_5\,.\,v_3\\[4pt] \frac{dX_6}{dt} &amp;= v_3-v_6\\[4pt] \frac{dX_7}{dt}&amp;=Y_6\,.\,v_5\\[4pt] \end{align*}" /></p>

<p>Where the <img src="http://latex.codecogs.com/gif.latex?X_n" title="X_n" /> are the different chemical and microorganisms species involved (even sub-indexes for chemical compounds), the <img src="http://latex.codecogs.com/gif.latex?Y_n" title="Y_n" /> are the yield coefficients and the <img src="http://latex.codecogs.com/gif.latex?v_n" title="v_n" /> are the pseudo-reactions:</p>

<p><img src="http://latex.codecogs.com/gif.latex?%5Cbegin%7Balign*%7D&space;v_1&amp;=%5Cfrac%7Bk_1%5C,.%5C,X_1%7D%7Bk_2&plus;X_1%7D%5C,.%5C,X_2%5C%5C[4pt]&space;v_2&amp;=k_7%5C,.%5C,X_3%5C,.%5C,X_2%5C%5C[4pt]&space;v_3&amp;=%5Cfrac%7Bk_3%5C,.%5C,X_1%7D%7Bk_4&plus;X_1%7D%5C,.%5C,X_4%5C%5C[4pt]&space;v_4&amp;=k_8%5C,.%5C,X_3%5C,.%5C,X_4%5C%5C[4pt]&space;v_5&amp;=%5Cfrac%7Bk_5%5C,.%5C,X_3%7D%7Bk_6&plus;X_3%7D%5C,.%5C,X_5%5C%5C[4pt]&space;v_6&amp;=k_9%5C,.%5C,X_5%5C%5C[4pt]&space;%5Cend%7Balign*%7D" title="\begin{align*} v_1&amp;=\frac{k_1\,.\,X_1}{k_2+X_1}\,.\,X_2\\[4pt] v_2&amp;=k_7\,.\,X_3\,.\,X_2\\[4pt] v_3&amp;=\frac{k_3\,.\,X_1}{k_4+X_1}\,.\,X_4\\[4pt] v_4&amp;=k_8\,.\,X_3\,.\,X_4\\[4pt] v_5&amp;=\frac{k_5\,.\,X_3}{k_6+X_3}\,.\,X_5\\[4pt] v_6&amp;=k_9\,.\,X_5\\[4pt] \end{align*}" /></p>

<p>I am using Pyomo for the estimation of all my unknown parameters, which are basically all the yield coefficients and kinetic constants (15 in total).</p>

<p>The following code works perfectly when is used with complete experimental time series for each of the dynamical variables:</p>

<pre><code>from pyomo.environ import *
from pyomo.dae import *

m = AbstractModel()
m.t = ContinuousSet()   
m.MEAS_t = Set(within=m.t)  # Measurement times, must be subset of t
m.x1_meas = Param(m.MEAS_t)
m.x2_meas = Param(m.MEAS_t)
m.x3_meas = Param(m.MEAS_t)
m.x4_meas = Param(m.MEAS_t)
m.x5_meas = Param(m.MEAS_t)
m.x6_meas = Param(m.MEAS_t)
m.x7_meas = Param(m.MEAS_t)

m.x1 = Var(m.t,within=PositiveReals)
m.x2 = Var(m.t,within=PositiveReals)
m.x3 = Var(m.t,within=PositiveReals)
m.x4 = Var(m.t,within=PositiveReals)
m.x5 = Var(m.t,within=PositiveReals)
m.x6 = Var(m.t,within=PositiveReals)
m.x7 = Var(m.t,within=PositiveReals)

m.k1 = Var(within=PositiveReals)
m.k2 = Var(within=PositiveReals)
m.k3 = Var(within=PositiveReals)
m.k4 = Var(within=PositiveReals)
m.k5 = Var(within=PositiveReals)
m.k6 = Var(within=PositiveReals)
m.k7 = Var(within=PositiveReals)
m.k8 = Var(within=PositiveReals)
m.k9 = Var(within=PositiveReals)
m.y1 = Var(within=PositiveReals)
m.y2 = Var(within=PositiveReals)
m.y3 = Var(within=PositiveReals)
m.y4 = Var(within=PositiveReals)
m.y5 = Var(within=PositiveReals)
m.y6 = Var(within=PositiveReals)

m.x1dot = DerivativeVar(m.x1,wrt=m.t)
m.x2dot = DerivativeVar(m.x2,wrt=m.t)
m.x3dot = DerivativeVar(m.x3,wrt=m.t)
m.x4dot = DerivativeVar(m.x4,wrt=m.t)
m.x5dot = DerivativeVar(m.x5,wrt=m.t)
m.x6dot = DerivativeVar(m.x6,wrt=m.t)
m.x7dot = DerivativeVar(m.x7,wrt=m.t)

def _init_conditions(m):
    yield m.x1[0] == 51.963
    yield m.x2[0] == 6.289
    yield m.x3[0] == 0
    yield m.x4[0] == 6.799
    yield m.x5[0] == 0
    yield m.x6[0] == 4.08
    yield m.x7[0] == 0
m.init_conditions=ConstraintList(rule=_init_conditions)


def _x1dot(m,i):
    if i==0:
        return Constraint.Skip
    return m.x1dot[i] == - m.y1*m.k1*m.x1[i]*m.x2[i]/(m.k2+m.x1[i]) - m.y2*m.k3*m.x1[i]*m.x4[i]/(m.k4+m.x1[i])
m.x1dotcon = Constraint(m.t, rule=_x1dot)

def _x2dot(m,i):
    if i==0:
        return Constraint.Skip
    return m.x2dot[i] ==  m.k1*m.x1[i]*m.x2[i]/(m.k2+m.x1[i]) - m.k7*m.x2[i]*m.x3[i]
m.x2dotcon = Constraint(m.t, rule=_x2dot)

def _x3dot(m,i):
    if i==0:
        return Constraint.Skip
    return m.x3dot[i] ==  m.y3*m.k1*m.x1[i]*m.x2[i]/(m.k2+m.x1[i]) - m.y4*m.k5*m.x3[i]*m.x6[i]/(m.k6+m.x3[i])
m.x3dotcon = Constraint(m.t, rule=_x3dot)

def _x4dot(m,i):
    if i==0:
        return Constraint.Skip
    return m.x4dot[i] == m.k3*m.x1[i]*m.x4[i]/(m.k4+m.x1[i]) - m.k8*m.x4[i]*m.x3[i]
m.x4dotcon = Constraint(m.t, rule=_x4dot)

def _x5dot(m,i):
    if i==0:
        return Constraint.Skip
    return m.x5dot[i] == m.y5*m.k3*m.x1[i]*m.x4[i]/(m.k4+m.x1[i])
m.x5dotcon = Constraint(m.t, rule=_x5dot)

def _x6dot(m,i):
    if i==0:
        return Constraint.Skip
    return m.x6dot[i] == m.k5*m.x3[i]*m.x6[i]/(m.k6+m.x3[i]) - m.k9*m.x6[i]*m.x7[i]
m.x6dotcon = Constraint(m.t, rule=_x6dot)

def _x7dot(m,i):
    if i==0:
        return Constraint.Skip
    return m.x7dot[i] == m.y6*m.k5*m.x3[i]*m.x6[i]/(m.k6+m.x3[i])
m.x7dotcon = Constraint(m.t, rule=_x7dot)

def _obj(m):
    return sum((m.x1[i]-m.x1_meas[i])**2+(m.x2[i]-m.x2_meas[i])**2+(m.x3[i]-m.x3_meas[i])**2+(m.x4[i]-m.x4_meas[i])**2+(m.x5[i]-m.x5_meas[i])**2+(m.x6[i]-m.x6_meas[i])**2+(m.x7[i]-m.x7_meas[i])**2 for i in m.MEAS_t)
m.obj = Objective(rule=_obj)

m.pprint()

instance = m.create_instance('exp.dat')
instance.t.pprint()

discretizer = TransformationFactory('dae.collocation')
discretizer.apply_to(instance,nfe=30)#,ncp=3)

solver=SolverFactory('ipopt')

results = solver.solve(instance,tee=True)
</code></pre>

<p>However, I am trying to run the same estimation routine in another experimental data that have missing values at the end of one or maximum two time series of some of the dynamical variables.</p>

<p>In other words, these complete experimental data looks like (in the .dat file):</p>

<pre><code>set t := 0  6   12  18  24  30  36  42  48  54  60  66  72  84  96  120 144;
set MEAS_t := 0 6   12  18  24  30  36  42  48  54  60  66  72  84  96  120 144;
param x1_meas :=
0   51.963
6   43.884
12  24.25
18  26.098
24  11.871
30  4.607
36  1.714
42  4.821
48  5.409
54  3.701
60  3.696
66  1.544
72  4.428
84  1.086
96  2.337
120 2.837
144 3.486
;
param x2_meas :=
0   6.289
6   6.242
12  7.804
18  7.202
24  6.48
30  5.833
36  6.644
42  5.741
48  4.568
54  4.252
60  5.603
66  5.167
72  4.399
84  4.773
96  4.801
120 3.866
144 3.847
;
param x3_meas :=
0   0
6   2.97
12  9.081
18  9.62
24  6.067
30  11.211
36  16.213
42  10.215
48  20.106
54  22.492
60  5.637
66  5.636
72  13.85
84  4.782
96  9.3
120 4.267
144 7.448
;
param x4_meas :=
0   6.799
6   7.73
12  7.804
18  8.299
24  8.208
30  8.523
36  8.507
42  8.656
48  8.49
54  8.474
60  8.203
66  8.127
72  8.111
84  8.064
96  6.845
120 6.721
144 6.162
;
param x5_meas :=
0   0
6   0.267
12  0.801
18  1.256
24  1.745
30  5.944
36  3.246
42  7.787
48  7.991
54  6.943
60  8.593
66  8.296
72  6.85
84  8.021
96  7.667
120 7.209
144 8.117
;
param x6_meas :=
0   4.08
6   4.545
12  4.784
18  4.888
24  5.293
30  5.577
36  5.802
42  5.967
48  6.386
54  6.115
60  6.625
66  6.835
72  6.383
84  6.605
96  5.928
120 5.354
144 4.975
;
param x7_meas :=
0   0
6   0.152
12  1.616
18  0.979
24  4.033
30  5.121
36  2.759
42  3.541
48  4.278
54  4.141
60  6.139
66  3.219
72  5.319
84  4.328
96  3.621
120 4.208
144 5.93
;
</code></pre>

<p>While one of my incomplete data sets could have all time series complete, but one like this:</p>

<pre><code>param x6_meas :=
0   4.08
6   4.545
12  4.784
18  4.888
24  5.293
30  5.577
36  5.802
42  5.967
48  6.386
54  6.115
60  6.625
66  6.835
72  6.383
84  6.605
96  5.928
120 5.354
144 .
;
</code></pre>

<p>I have knowledge that one can specify to Pyomo to take the derivative of certain variables with respect to a different time serie. However, after tried it, it hadn't worked, and I guess that is because that these are coupled ODE. So basically my question is if there is a way to overcome this issue in Pyomo.</p>

<p>Thanks in advance.</p>

<p>I think all you need to do is slightly modify your objective function like this:</p>

<pre><code>def _obj(m):
    sum1 = sum((m.x1[i]-m.x1_meas[i])**2 for i in m.MEAS_t if i in m.x1_meas.keys())
    sum2 = sum((m.x2[i]-m.x2_meas[i])**2 for i in m.MEAS_t if i in m.x2_meas.keys())
    sum3 = sum((m.x3[i]-m.x3_meas[i])**2 for i in m.MEAS_t if i in m.x3_meas.keys())
    sum4 = sum((m.x4[i]-m.x4_meas[i])**2 for i in m.MEAS_t if i in m.x4_meas.keys())
    sum5 = sum((m.x5[i]-m.x5_meas[i])**2 for i in m.MEAS_t if i in m.x5_meas.keys())
    sum6 = sum((m.x6[i]-m.x6_meas[i])**2 for i in m.MEAS_t if i in m.x6_meas.keys())
    sum7 = sum((m.x7[i]-m.x7_meas[i])**2 for i in m.MEAS_t if i in m.x7_meas.keys())
    return sum1+sum2+sum3+sum4+sum5+sum6+sum7
m.obj = Objective(rule=_obj)
</code></pre>

<p>This double checks that <em>i</em> is a valid index for each set of measurements before adding that index to the sum. If you knew apriori which measurement sets were missing data then you could simplify this function by only doing this check on those sets and summing over the others like you were before. </p>

python ode pyomo