2016-05-26 2 views
0

Je suis nouveau en utilisant Pyomo, donc je m'excuse d'avance si c'est une question de base. Eh bien, je travaille avec des modèles cinétiques et mon but est d'estimer les paramètres cinétiques. J'ai commencé avec un 'modèle de jouet' pour mieux comprendre Pyomo avant d'essayer mon compliqué.Python Pyomo: estimation de paramètres dans un système ODE

Alors, mon modèle de jouet est un système simple de ODE de 3 équations:

dX1/dt = -k1*X1 
dX2/dt = k1*X1 - k2*X2 
dX3/dt = k2*X2 

Mon but est d'estimer les paramètres k1 et k2. Je légèrement changé le code de ce tutorial et est la suivante:

from pyomo.environ import * 
from pyomo.dae import * 

model = AbstractModel() 
model.t = ContinuousSet() 
model.MEAS_t = Set(within=model.t) 
model.x1_meas = Param(model.MEAS_t) 
model.x2_meas = Param(model.MEAS_t) 
model.x3_meas = Param(model.MEAS_t) 

model.x1 = Var(model.t) 
model.x2 = Var(model.t) 
model.x3 = Var(model.t) 

model.k1 = Var(bounds=(0,3)) 
model.k2 = Var(bounds=(0,3)) 

model.x1dot = DerivativeVar(model.x1,wrt=model.t) 
model.x2dot = DerivativeVar(model.x2,wrt=model.t) 
model.x3dot = DerivativeVar(model.x3,wrt=model.t) 

def _x1dot(model,i): 
    return model.x1dot[i] == -model.k1*model.x1[i] 
model.x1dotcon = Constraint(model.t, rule=_x1dot) 

def _x2dot(model,i): 
    return model.x2dot[i] == model.k1*model.x1[i]-model.k2*model.x2[i] 
model.x2dotcon = Constraint(model.t, rule=_x2dot) 

def _x3dot(model,i): 
    return model.x3dot[i] == model.k2*model.x2[i] 
model.x3dotcon = Constraint(model.t, rule=_x3dot) 


def _obj(model): 
    return sum((model.x1[i]-model.x1_meas[i])**2+(model.x2[i]-model.x2_meas[i])**2+(model.x3[i]-model.x3_meas)**2 for i in model.MEAS_t) 
model.obj = Objective(rule=_obj) 

model.pprint() 

instance = model.create_instance('data2.dat') 
instance.t.pprint() 

discretizer = TransformationFactory('dae.collocation') 
discretizer.apply_to(instance,nfe=8,ncp=5) 

solver=SolverFactory('ipopt') 

results = solver.solve(instance,tee=True) 

instance.k1.pprint() 
instance.k2.pprint() 

Une fois que je lance ce code, je reçu le message suivant:

TypeError: Cannot convert object of type 'IndexedParam' (value = x3_meas) to a numeric value. 

Cependant, quand j'efface toutes les lignes correspondant à x3_meas dans mon code, ainsi que les données dans mon fichier .dat, cela fonctionne parfaitement.

Est-ce que quelqu'un sait quel est le problème?

Mes données ressemble:

set t := 0.00 0.66 1.33 2.00 2.66 3.33 4.00 4.66 5.33 6.00 ; 
set MEAS_t := 0.00 0.66 1.33 2.00 2.66 3.33 4.00 4.66 5.33 6.00 ; 
param x1_meas := 
0.00 1.000000 
0.66 0.263597 
1.33 0.069483 
2.00 0.018316 
2.66 0.004828 
3.33 0.001273 
4.00 0.000335 
4.66 0.000088 
5.33 0.000023 
6.00 0.000006 
; 
param x2_meas := 
0.00 0.000000 
0.66 0.499640 
1.33 0.388227 
2.00 0.234039 
2.66 0.129311 
3.33 0.068803 
4.00 0.035960 
4.66 0.018630 
5.33 0.009609 
6.00 0.004945 
; 
param x3_meas := 
0.00 0.000000 
0.66 0.236763 
1.33 0.542289 
2.00 0.747645 
2.66 0.865861 
3.33 0.929925 
4.00 0.963704 
4.66 0.981281 
5.33 0.990367 
6.00 0.995049 
; 

Répondre

0

Après avoir demandé de l'aide dans le groupe Google de Pyomo, je fait remarquer que je commis une erreur stupide:

Dans votre fonction objectif vous avez oublié de x3_meas d'index. Votre fonction d'objectif doit être:

def _obj(model): 
    return sum((model.x1[i]-model.x1_meas[i])**2+(model.x2[i]-model.x2_meas[i])**2+(model.x3[i]-model.x3_meas[i])**2 for i in model.MEAS_t)