I am trying to solve DAE using pyomo. However, I am getting Keyerrors.The material balance equation in chemical engineering, which includes differential equations.
import pyomo.environ as pyo
from pyomo import dae
from pyomo import opt
import numpy as np
gasname = ['Gas1', 'Gas2']
x_init_gas1 = 0.2
x_init_gas2 = 0.8
x_init_list = [x_init_gas1, x_init_gas2]
q_list = [1.2e-6, 5.0e-5]
x_init_dict = dict(zip(gasname, x_init_list))
q_dict = dict(zip(gasname, q_list))
area1_start = 0
area1_end = 100
ph_1, pl_1 = [3, 1]
U = 20
model = pyo.ConcreteModel(name='membrane sepalation')
model.gas_name = pyo.Set(initialize=gasname)
model.area1 = dae.ContinuousSet(bounds=(area1_start, area1_end))
discretizer_area1 = pyo.TransformationFactory('dae.finite_difference')
discretizer_area1.apply(model, wrt=model.area1)
area1_delta = list(model.area1)
model.area1_len = pyo.Set(initialize=range(0, len(area1_delta)))
model.ode_F1 = pyo.Var(model.gas_name, model.area1, domain=pyo.PositiveReals, initialize=(10))
model.F1_x = pyo.Var(model.gas_name, bounds=(0, 30), initialize=(10))
model.V1_x = pyo.Var(model.gas_name, bounds=(0, 30), initialize=(10))
model.G1_x = pyo.Var(model.gas_name, bounds=(0, 30), initialize=(10))
model.F1 = sum(model.F1_x[s] for s in gasname)
model.V1 = sum(model.V1_x[s] for s in gasname)
model.G1 = sum(model.G1_x[s] for s in gasname)
ph_1, pl_1 = [2.98, 1]
# 1
def const_F1_x(m, gas):
return m.F1_x[gas] == x_init_dict[gas]*U
model.const_F1_x = pyo.Constraint(model.gas_name, rule=const_F1_x)
def const_F1_flow(m, gas):
return m.F1_x[gas] == m.V1_x[gas] + m.G1_x[gas]
model.const_F1_flow = pyo.Constraint(model.gas_name, rule=const_F1_flow)
def const_module1_flow(m):
return m.F1 == m.V1 + m.G1
model.const_module1_flow = pyo.Constraint(rule=const_module1_flow)
model.dF1_da = dae.DerivativeVar(model.ode_F1, wrt=model.area1, domain=pyo.PositiveReals)
def _ode_rule_module1(m, a, gas):
if a == 0:
return pyo.Constraint.Skip
i = area1_delta.index(a)
return m.dF1_da[gas, a] == -q_dict[gas]*((ph_1*m.ode_F1[gas, a]/sum(m.ode_F1[:, a]))-pl_1*((m.ode_F1[gas, a] - m.ode_F1[gas, area1_delta[i-1]])/(sum(m.ode_F1[:, a])-sum(m.ode_F1[:, area1_delta[i-1]])+1e-10)))
model.ode_rule_module1 = pyo.Constraint(model.area1, model.gas_name, rule=_ode_rule_module1)
def _ode_rule_module1_entry(m, gas):
return m.F1_x[gas] == m.ode_F1[gas, area1_start]
model.ode_rule_module1_entry = pyo.Constraint(model.gas_name, rule=_ode_rule_module1_entry)
def _ode_rule_module1_exit(m, gas):
return m.V1_x[gas] == m.ode_F1[gas, area1_end]
model.ode_rule_module1_exit = pyo.Constraint(model.gas_name, rule=_ode_rule_module1_exit)
def calc_objective(m):
return m.V1_x[gasname[0]] / (U*x_init_gas1)
model.OBJ = pyo.Objective(rule=calc_objective, sense=pyo.maximize)
solver = opt.SolverFactory('ipopt')
solver.options["print_level"] = 1
res = solver.solve(model, tee=True)
print(model.display())
and, the error I got was
File ~.py:154 in <module>
res = solver.solve(model, tee=True)
File ~\anaconda3\lib\site-packages\pyomo\opt\base\solvers.py:570 in solve
self._presolve(*args, **kwds)
~~~
File pyomo\repn\plugins\ampl\ampl_.pyx:1039 in pyomo.repn.plugins.ampl.ampl_.ProblemWriter_nl._print_model_NL
File pyomo\repn\plugins\ampl\ampl_.pyx:1033 in genexpr
File pyomo\repn\plugins\ampl\ampl_.pyx:1033 in genexpr
KeyError: 2447832121968
The model.display() showed like that.
4 Set Declarations
area1_len : Size=1, Index=None, Ordered=Insertion
Key : Dimen : Domain : Size : Members
None : 1 : Any : 11 : {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
gas_name : Size=1, Index=None, Ordered=Insertion
Key : Dimen : Domain : Size : Members
None : 1 : Any : 2 : {'Gas1', 'Gas2'}
ode_F1_index : Size=1, Index=None, Ordered=True
Key : Dimen : Domain : Size : Members
None : 2 : gas_name*area1 : 22 : {('Gas1', 0), ('Gas1', 10.0), ('Gas1', 20.0), ('Gas1', 30.0), ('Gas1', 40.0), ('Gas1', 50.0), ('Gas1', 60.0), ('Gas1', 70.0), ('Gas1', 80.0), ('Gas1', 90.0), ('Gas1', 100), ('Gas2', 0), ('Gas2', 10.0), ('Gas2', 20.0), ('Gas2', 30.0), ('Gas2', 40.0), ('Gas2', 50.0), ('Gas2', 60.0), ('Gas2', 70.0), ('Gas2', 80.0), ('Gas2', 90.0), ('Gas2', 100)}
ode_rule_module1_index : Size=1, Index=None, Ordered=True
Key : Dimen : Domain : Size : Members
None : 2 : area1*gas_name : 22 : {(0, 'Gas1'), (0, 'Gas2'), (10.0, 'Gas1'), (10.0, 'Gas2'), (20.0, 'Gas1'), (20.0, 'Gas2'), (30.0, 'Gas1'), (30.0, 'Gas2'), (40.0, 'Gas1'), (40.0, 'Gas2'), (50.0, 'Gas1'), (50.0, 'Gas2'), (60.0, 'Gas1'), (60.0, 'Gas2'), (70.0, 'Gas1'), (70.0, 'Gas2'), (80.0, 'Gas1'), (80.0, 'Gas2'), (90.0, 'Gas1'), (90.0, 'Gas2'), (100, 'Gas1'), (100, 'Gas2')}
1 RangeSet Declarations
area1_domain : Dimen=1, Size=Inf, Bounds=(0, 100)
Key : Finite : Members
None : False : [0..100]
4 Var Declarations
F1_x : Size=2, Index=gas_name
Key : Lower : Value : Upper : Fixed : Stale : Domain
Gas1 : 0 : 10 : 30 : False : False : Reals
Gas2 : 0 : 10 : 30 : False : False : Reals
G1_x : Size=2, Index=gas_name
Key : Lower : Value : Upper : Fixed : Stale : Domain
Gas1 : 0 : 10 : 30 : False : False : Reals
Gas2 : 0 : 10 : 30 : False : False : Reals
V1_x : Size=2, Index=gas_name
Key : Lower : Value : Upper : Fixed : Stale : Domain
Gas1 : 0 : 10 : 30 : False : False : Reals
Gas2 : 0 : 10 : 30 : False : False : Reals
ode_F1 : Size=22, Index=ode_F1_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
('Gas1', 0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 10.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 20.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 30.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 40.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 50.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 60.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 70.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 80.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 90.0) : 0 : 10 : None : False : False : PositiveReals
('Gas1', 100) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 10.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 20.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 30.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 40.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 50.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 60.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 70.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 80.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 90.0) : 0 : 10 : None : False : False : PositiveReals
('Gas2', 100) : 0 : 10 : None : False : False : PositiveReals
1 Objective Declarations
OBJ : Size=1, Index=None, Active=True
Key : Active : Sense : Expression
None : True : maximize : 0.2557544757033248*V1_x[Gas1]
6 Constraint Declarations
const_F1_flow : Size=2, Index=gas_name, Active=True
Key : Lower : Body : Upper : Active
Gas1 : 0.0 : F1_x[Gas1] - (V1_x[Gas1] + G1_x[Gas1]) : 0.0 : True
Gas2 : 0.0 : F1_x[Gas2] - (V1_x[Gas2] + G1_x[Gas2]) : 0.0 : True
const_F1_x : Size=2, Index=gas_name, Active=True
Key : Lower : Body : Upper : Active
Gas1 : 3.91 : F1_x[Gas1] : 3.91 : True
Gas2 : 15.64 : F1_x[Gas2] : 15.64 : True
const_module1_flow : Size=1, Index=None, Active=True
Key : Lower : Body : Upper : Active
None : 0.0 : F1_x[Gas1] + F1_x[Gas2] - (V1_x[Gas1] + V1_x[Gas2] + G1_x[Gas1] + G1_x[Gas2]) : 0.0 : True
ode_rule_module1 : Size=20, Index=ode_rule_module1_index, Active=True
Key : Lower : Body : Upper : Active
(10.0, 'Gas1') : 0.0 : dF1_da[Gas1,10.0] + 1.2e-06*(2.98*ode_F1[Gas1,10.0]/(ode_F1[Gas1,10.0] + ode_F1[Gas2,10.0]) - (ode_F1[Gas1,10.0] - ode_F1[Gas1,0])/(ode_F1[Gas1,10.0] + ode_F1[Gas2,10.0] - (ode_F1[Gas1,0] + ode_F1[Gas2,0]) + 1e-10)) : 0.0 : True
(10.0, 'Gas2') : 0.0 : dF1_da[Gas2,10.0] + 5e-05*(2.98*ode_F1[Gas2,10.0]/(ode_F1[Gas1,10.0] + ode_F1[Gas2,10.0]) - (ode_F1[Gas2,10.0] - ode_F1[Gas2,0])/(ode_F1[Gas1,10.0] + ode_F1[Gas2,10.0] - (ode_F1[Gas1,0] + ode_F1[Gas2,0]) + 1e-10)) : 0.0 : True
(20.0, 'Gas1') : 0.0 : dF1_da[Gas1,20.0] + 1.2e-06*(2.98*ode_F1[Gas1,20.0]/(ode_F1[Gas1,20.0] + ode_F1[Gas2,20.0]) - (ode_F1[Gas1,20.0] - ode_F1[Gas1,10.0])/(ode_F1[Gas1,20.0] + ode_F1[Gas2,20.0] - (ode_F1[Gas1,10.0] + ode_F1[Gas2,10.0]) + 1e-10)) : 0.0 : True
(20.0, 'Gas2') : 0.0 : dF1_da[Gas2,20.0] + 5e-05*(2.98*ode_F1[Gas2,20.0]/(ode_F1[Gas1,20.0] + ode_F1[Gas2,20.0]) - (ode_F1[Gas2,20.0] - ode_F1[Gas2,10.0])/(ode_F1[Gas1,20.0] + ode_F1[Gas2,20.0] - (ode_F1[Gas1,10.0] + ode_F1[Gas2,10.0]) + 1e-10)) : 0.0 : True
(30.0, 'Gas1') : 0.0 : dF1_da[Gas1,30.0] + 1.2e-06*(2.98*ode_F1[Gas1,30.0]/(ode_F1[Gas1,30.0] + ode_F1[Gas2,30.0]) - (ode_F1[Gas1,30.0] - ode_F1[Gas1,20.0])/(ode_F1[Gas1,30.0] + ode_F1[Gas2,30.0] - (ode_F1[Gas1,20.0] + ode_F1[Gas2,20.0]) + 1e-10)) : 0.0 : True
(30.0, 'Gas2') : 0.0 : dF1_da[Gas2,30.0] + 5e-05*(2.98*ode_F1[Gas2,30.0]/(ode_F1[Gas1,30.0] + ode_F1[Gas2,30.0]) - (ode_F1[Gas2,30.0] - ode_F1[Gas2,20.0])/(ode_F1[Gas1,30.0] + ode_F1[Gas2,30.0] - (ode_F1[Gas1,20.0] + ode_F1[Gas2,20.0]) + 1e-10)) : 0.0 : True
(40.0, 'Gas1') : 0.0 : dF1_da[Gas1,40.0] + 1.2e-06*(2.98*ode_F1[Gas1,40.0]/(ode_F1[Gas1,40.0] + ode_F1[Gas2,40.0]) - (ode_F1[Gas1,40.0] - ode_F1[Gas1,30.0])/(ode_F1[Gas1,40.0] + ode_F1[Gas2,40.0] - (ode_F1[Gas1,30.0] + ode_F1[Gas2,30.0]) + 1e-10)) : 0.0 : True
(40.0, 'Gas2') : 0.0 : dF1_da[Gas2,40.0] + 5e-05*(2.98*ode_F1[Gas2,40.0]/(ode_F1[Gas1,40.0] + ode_F1[Gas2,40.0]) - (ode_F1[Gas2,40.0] - ode_F1[Gas2,30.0])/(ode_F1[Gas1,40.0] + ode_F1[Gas2,40.0] - (ode_F1[Gas1,30.0] + ode_F1[Gas2,30.0]) + 1e-10)) : 0.0 : True
(50.0, 'Gas1') : 0.0 : dF1_da[Gas1,50.0] + 1.2e-06*(2.98*ode_F1[Gas1,50.0]/(ode_F1[Gas1,50.0] + ode_F1[Gas2,50.0]) - (ode_F1[Gas1,50.0] - ode_F1[Gas1,40.0])/(ode_F1[Gas1,50.0] + ode_F1[Gas2,50.0] - (ode_F1[Gas1,40.0] + ode_F1[Gas2,40.0]) + 1e-10)) : 0.0 : True
(50.0, 'Gas2') : 0.0 : dF1_da[Gas2,50.0] + 5e-05*(2.98*ode_F1[Gas2,50.0]/(ode_F1[Gas1,50.0] + ode_F1[Gas2,50.0]) - (ode_F1[Gas2,50.0] - ode_F1[Gas2,40.0])/(ode_F1[Gas1,50.0] + ode_F1[Gas2,50.0] - (ode_F1[Gas1,40.0] + ode_F1[Gas2,40.0]) + 1e-10)) : 0.0 : True
(60.0, 'Gas1') : 0.0 : dF1_da[Gas1,60.0] + 1.2e-06*(2.98*ode_F1[Gas1,60.0]/(ode_F1[Gas1,60.0] + ode_F1[Gas2,60.0]) - (ode_F1[Gas1,60.0] - ode_F1[Gas1,50.0])/(ode_F1[Gas1,60.0] + ode_F1[Gas2,60.0] - (ode_F1[Gas1,50.0] + ode_F1[Gas2,50.0]) + 1e-10)) : 0.0 : True
(60.0, 'Gas2') : 0.0 : dF1_da[Gas2,60.0] + 5e-05*(2.98*ode_F1[Gas2,60.0]/(ode_F1[Gas1,60.0] + ode_F1[Gas2,60.0]) - (ode_F1[Gas2,60.0] - ode_F1[Gas2,50.0])/(ode_F1[Gas1,60.0] + ode_F1[Gas2,60.0] - (ode_F1[Gas1,50.0] + ode_F1[Gas2,50.0]) + 1e-10)) : 0.0 : True
(70.0, 'Gas1') : 0.0 : dF1_da[Gas1,70.0] + 1.2e-06*(2.98*ode_F1[Gas1,70.0]/(ode_F1[Gas1,70.0] + ode_F1[Gas2,70.0]) - (ode_F1[Gas1,70.0] - ode_F1[Gas1,60.0])/(ode_F1[Gas1,70.0] + ode_F1[Gas2,70.0] - (ode_F1[Gas1,60.0] + ode_F1[Gas2,60.0]) + 1e-10)) : 0.0 : True
(70.0, 'Gas2') : 0.0 : dF1_da[Gas2,70.0] + 5e-05*(2.98*ode_F1[Gas2,70.0]/(ode_F1[Gas1,70.0] + ode_F1[Gas2,70.0]) - (ode_F1[Gas2,70.0] - ode_F1[Gas2,60.0])/(ode_F1[Gas1,70.0] + ode_F1[Gas2,70.0] - (ode_F1[Gas1,60.0] + ode_F1[Gas2,60.0]) + 1e-10)) : 0.0 : True
(80.0, 'Gas1') : 0.0 : dF1_da[Gas1,80.0] + 1.2e-06*(2.98*ode_F1[Gas1,80.0]/(ode_F1[Gas1,80.0] + ode_F1[Gas2,80.0]) - (ode_F1[Gas1,80.0] - ode_F1[Gas1,70.0])/(ode_F1[Gas1,80.0] + ode_F1[Gas2,80.0] - (ode_F1[Gas1,70.0] + ode_F1[Gas2,70.0]) + 1e-10)) : 0.0 : True
(80.0, 'Gas2') : 0.0 : dF1_da[Gas2,80.0] + 5e-05*(2.98*ode_F1[Gas2,80.0]/(ode_F1[Gas1,80.0] + ode_F1[Gas2,80.0]) - (ode_F1[Gas2,80.0] - ode_F1[Gas2,70.0])/(ode_F1[Gas1,80.0] + ode_F1[Gas2,80.0] - (ode_F1[Gas1,70.0] + ode_F1[Gas2,70.0]) + 1e-10)) : 0.0 : True
(90.0, 'Gas1') : 0.0 : dF1_da[Gas1,90.0] + 1.2e-06*(2.98*ode_F1[Gas1,90.0]/(ode_F1[Gas1,90.0] + ode_F1[Gas2,90.0]) - (ode_F1[Gas1,90.0] - ode_F1[Gas1,80.0])/(ode_F1[Gas1,90.0] + ode_F1[Gas2,90.0] - (ode_F1[Gas1,80.0] + ode_F1[Gas2,80.0]) + 1e-10)) : 0.0 : True
(90.0, 'Gas2') : 0.0 : dF1_da[Gas2,90.0] + 5e-05*(2.98*ode_F1[Gas2,90.0]/(ode_F1[Gas1,90.0] + ode_F1[Gas2,90.0]) - (ode_F1[Gas2,90.0] - ode_F1[Gas2,80.0])/(ode_F1[Gas1,90.0] + ode_F1[Gas2,90.0] - (ode_F1[Gas1,80.0] + ode_F1[Gas2,80.0]) + 1e-10)) : 0.0 : True
(100, 'Gas1') : 0.0 : dF1_da[Gas1,100] + 1.2e-06*(2.98*ode_F1[Gas1,100]/(ode_F1[Gas1,100] + ode_F1[Gas2,100]) - (ode_F1[Gas1,100] - ode_F1[Gas1,90.0])/(ode_F1[Gas1,100] + ode_F1[Gas2,100] - (ode_F1[Gas1,90.0] + ode_F1[Gas2,90.0]) + 1e-10)) : 0.0 : True
(100, 'Gas2') : 0.0 : dF1_da[Gas2,100] + 5e-05*(2.98*ode_F1[Gas2,100]/(ode_F1[Gas1,100] + ode_F1[Gas2,100]) - (ode_F1[Gas2,100] - ode_F1[Gas2,90.0])/(ode_F1[Gas1,100] + ode_F1[Gas2,100] - (ode_F1[Gas1,90.0] + ode_F1[Gas2,90.0]) + 1e-10)) : 0.0 : True
ode_rule_module1_entry : Size=2, Index=gas_name, Active=True
Key : Lower : Body : Upper : Active
Gas1 : 0.0 : F1_x[Gas1] - ode_F1[Gas1,0] : 0.0 : True
Gas2 : 0.0 : F1_x[Gas2] - ode_F1[Gas2,0] : 0.0 : True
ode_rule_module1_exit : Size=2, Index=gas_name, Active=True
Key : Lower : Body : Upper : Active
Gas1 : 0.0 : V1_x[Gas1] - ode_F1[Gas1,100] : 0.0 : True
Gas2 : 0.0 : V1_x[Gas2] - ode_F1[Gas2,100] : 0.0 : True
1 ContinuousSet Declarations
area1 : Size=1, Index=None, Ordered=Sorted
Key : Dimen : Domain : Size : Members
None : 1 : [0..100] : 11 : {0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100}
1 DerivativeVar Declarations
dF1_da : Size=22, Index=ode_F1_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
('Gas1', 0) : 0 : None : None : False : True : PositiveReals
('Gas1', 10.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 20.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 30.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 40.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 50.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 60.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 70.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 80.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 90.0) : 0 : None : None : False : True : PositiveReals
('Gas1', 100) : 0 : None : None : False : True : PositiveReals
('Gas2', 0) : 0 : None : None : False : True : PositiveReals
('Gas2', 10.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 20.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 30.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 40.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 50.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 60.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 70.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 80.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 90.0) : 0 : None : None : False : True : PositiveReals
('Gas2', 100) : 0 : None : None : False : True : PositiveReals
18 Declarations: gas_name area1_domain area1 area1_len ode_F1_index ode_F1 F1_x V1_x G1_x const_F1_x const_F1_flow const_module1_flow dF1_da ode_rule_module1_index ode_rule_module1 ode_rule_module1_entry ode_rule_module1_exit OBJ
Similar questions existed before, but I could not understand them.
Could you please help me solve it? Thank you very much!
You're applying the discretization transformation before declaring all the DerivativeVar components in your model. The discretization transformation discretizes the ContinuousSets in the model and also adds discretization equations to the model for every DerivativeVar that has been declared so it is crucial that the transformation not be applied until after all DerivativeVar components have been declared.