中具有缺失值的ODE系統Pyomo的參數估計在時間序列

Question

我有7個方程的ODE系統用於說明一組特定的形式的微生物動力學：

凡是涉及不同的化學和微生物物種（對於化學化合物偶數子索引），則是產量係數和是僞反應：

我使用Pyomo來估計所有我未知的參數，它們基本上是所有的屈服係數和動力學常數（共15個）。

下面的代碼完美地工作時與每一個動態變量的完整的實驗時間序列使用：

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) 

不過，我試圖運行在另一個實驗數據相同的推斷程序有遺漏值在一個或最多兩個時間序列的一些動態變量的末尾。

換句話說，這些完整的實驗數據看起來像（在.dat文件）：

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 
; 

雖然我不完全數據集一個可以將所有時間序列完整，但有這樣的：

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 . 
; 

我知道可以指定Pyomo對不同時間系列採取某些變量的導數。然而，經過嘗試後，它沒有奏效，我想這是因爲這些是耦合ODE。所以基本上我的問題是在Pyomo中是否有辦法解決這個問題。

在此先感謝。

Answer 1

我認爲，所有你需要做的是稍微修改你的目標函數是這樣的：

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) 

這種雙重檢查是我是每個組測量有效的索引，並稱指數總和前。如果您知道哪些測量數據缺失數據，那麼您可以通過僅對這些數據集進行檢查並像以前一樣對其他數據進行求和來簡化此功能。