我有以下使用python pyomo的代码。我认为我的代码是正确的，但它不能运行。我不明白怎么了。False statement casadi/core/integrator.cpp:1417: Assertion "!de_in[DE_X].is_empty()" failed:Ill-posed ODE - no state。
这是我的代码，为了简单起见，我只放了一部分代码:
def Model_8(t_i=0,t_e=600,dt=1,CAT0=1e-6,CoC=1e-3,H2=0,M1=2,M2=0.9,Volume=Volume_135_temp,temp=temp_135,temp_i=temp_i_135,_M1concen=_M1concen_135_temp,k0a={1:10,2:100,3:1000},k0ij={1:{1:0.1,2:0},2:{1:10,2:0},3:{1:50,2:0}},k0pij={1:{(1,1):2500,(1,2):800,(2,1):1200,(2,2):300},2:{(1,1):5500,(1,2):800,(2,1):3200,(2,2):300},3:{(1,1):10500,(1,2):800,(2,1):5200,(2,2):300}},k0tCoC={1:1e-3,2:1e-3,3:1e-3},k0tij={1:{(1,1):1e-2,(1,2):1e-2,(2,1):1e-2,(2,2):1e-2},2:{(1,1):1e-2,(1,2):1e-2,(2,1):1e-2,(2,2):1e-2},3:{(1,1):1e-2,(1,2):1e-2,(2,1):1e-2,(2,2):1e-2}},k0tH={1:0,2:0,3:0},k0tbeta={1:1,2:5,3:5},k0d={1:0.004,2:0.008,3:0.01},Tref=405,E={1:10000,2:20000,3:1000},R=8.314,MW_i={1:28,2:112},sites=3):
model = pyo.ConcreteModel()
# Sets
model.t = pyodae.ContinuousSet(bounds=(t_i,t_e),initialize= temp_i)
model.j = pyo.RangeSet(1,2)
model.site = pyo.RangeSet(1,sites+1,1)
#Variables
model.CAT = pyo.Var(model.t,model.site, domain = pyo.NonNegativeReals,initialize=CAT0)
model.Ca = pyo.Var(model.t,model.site, domain = pyo.NonNegativeReals,initialize=0)
model.Cd = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=0)
model.M1 = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=M1)
model.M2 = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=M2)
model.P01 = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=0)
model.P02 = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=0)
model.P11 = pyo.Var(model.t,model.site, domain = pyo.NonNegativeReals,initialize=0)
model.P12 = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=0)
model.P21 = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=0)
model.P22 = pyo.Var(model.t,model.site, domain = pyo.NonNegativeReals,initialize=0)
model.D0 = pyo.Var(model.t,model.site, domain = pyo.NonNegativeReals,initialize=0)
model.D1 = pyo.Var(model.t,model.site, domain = pyo.NonNegativeReals,initialize=0)
model.D2 = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=0)
model.qM1 = pyo.Var(model.t, model.site,domain = pyo.NonNegativeReals,initialize=0)
# Parameters
model.CoC = pyo.Param(model.site,initialize=CoC)
model.H2 = pyo.Param(model.site,initialize=H2)
model.M1concen = pyo.Param(model.t,initialize=_M1concen[0])
# Create dictionaries to store the collections
_ka = {}
_kij1 = {}
_kij2 = {}
_kpij11 = {}
_kpij12 = {}
_kpij21 = {}
_kpij22 = {}
_ktCoC = {}
_ktij11 = {}
_ktij12 = {}
_ktij21 = {}
_ktij22 = {}
_ktH = {}
_ktbeta = {}
_kd = {}
# Loop through sites and populate the dictionaries
for site in range(1, sites+1):
_ka[site] = {}
_kij1[site] = {}
_kij2[site] = {}
_kpij11[site] = {}
_kpij12[site] = {}
_kpij21[site] = {}
_kpij22[site] = {}
_ktCoC[site] = {}
_ktij11[site] = {}
_ktij12[site] = {}
_ktij21[site] = {}
_ktij22[site] = {}
_ktH[site] = {}
_ktbeta[site] = {}
_kd[site] = {}
for i in temp_i:
_ka[site][i] = k0a[site]*exp((0/R)*((1/temp[i])-(1/Tref)))
_kij1[site][i] = k0ij[site][1]*exp((0/R)*((1/temp[i])-(1/Tref)))
_kij2[site][i] = k0ij[site][2]*exp((0/R)*((1/temp[i])-(1/Tref)))
_kpij11[site][i] = k0pij[site][1,1]*exp((-E[1]/R)*((1/temp[i])-(1/Tref)))
_kpij12[site][i] = k0pij[site][1,2]*exp((-E[1]/R)*((1/temp[i])-(1/Tref)))
_kpij21[site][i] = k0pij[site][2,1]*exp((-E[1]/R)*((1/temp[i])-(1/Tref)))
_kpij22[site][i] = k0pij[site][2,2]*exp((-E[1]/R)*((1/temp[i])-(1/Tref)))
_ktCoC[site][i] = k0tCoC[site] * exp((-E[2]/R)*((1/temp[i])-(1/Tref)))
_ktij11[site][i] = k0tij[site][1,1]*exp((-E[2]/R)*((1/temp[i])-(1/Tref)))
_ktij12[site][i] = k0tij[site][1,2]*exp((-E[2]/R)*((1/temp[i])-(1/Tref)))
_ktij21[site][i] = k0tij[site][2,1]*exp((-E[2]/R)*((1/temp[i])-(1/Tref)))
_ktij22[site][i] = k0tij[site][2,2]*exp((-E[2]/R)*((1/temp[i])-(1/Tref)))
_ktH[site][i] = k0tH[site]*exp((-E[2]/R)*((1/temp[i])-(1/Tref)))
_ktbeta[site][i] = k0tbeta[site]*exp((-E[2]/R)*((1/temp[i])-(1/Tref)))
_kd[site][i] = k0d[site]*exp((-E[2]/R)*((1/temp[i])-(1/Tref)))
data = {
f"site_{k}": {
"ka": list(_ka[k].values()),
"kij1": list(_kij1[k].values()),
"kij2": list(_kij2[k].values()),
"kpij11": list(_kpij11[k].values()),
"kpij12": list(_kpij12[k].values()),
"kpij21": list(_kpij21[k].values()),
"kpij22": list(_kpij22[k].values()),
"ktCoC": list(_ktCoC[k].values()),
"ktij11": list(_ktij11[k].values()),
"ktij12": list(_ktij12[k].values()),
"ktij21": list(_ktij21[k].values()),
"ktij22": list(_ktij22[k].values()),
"ktH": list(_ktH[k].values()),
"ktbeta": list(_ktbeta[k].values()),
"kd": list(_kd[k].values()),
"T": list(temp)
}
for k in range(1, sites+1)
}
site_k = {k: pd.DataFrame(data[k]) for k in data.keys()}
model.vars = pyo.Param([(t,key,site) for t in model.t for key in ["ka", "kij1","kij2", "kpij11","kpij12","kpij21","kpij22","ktij11","ktij12","ktij21","ktij22","ktbeta","ktH","ktCoC","kd"] for site in range(1,sites+1)],mutable=True)
for site in range(1,sites+1,1):
for t in model.t:
model.vars[t,"ka",site] = _ka[site][t]
model.vars[t,"kij1",site] = _kij1[site][t]
model.vars[t,"kij2",site] = _kij2[site][t]
model.vars[t,"kpij11",site] = _kpij11[site][t]
model.vars[t,"kpij12",site] = _kpij12[site][t]
model.vars[t,"kpij21",site] = _kpij21[site][t]
model.vars[t,"kpij22",site] = _kpij22[site][t]
model.vars[t,"ktij11",site] = _ktij11[site][t]
model.vars[t,"ktij12",site] = _ktij12[site][t]
model.vars[t,"ktij21",site] = _ktij21[site][t]
model.vars[t,"ktij22",site] = _ktij22[site][t]
model.vars[t,"ktCoC",site] = _ktCoC[site][t]
model.vars[t,"ktbeta",site] = _ktbeta[site][t]
model.vars[t,"ktH",site] = _ktH[site][t]
model.vars[t,"kd",site] = _kd[site][t]
model.CATdt = pyodae.DerivativeVar(model.CAT, withrespectto=model.t)
model.Cadt = pyodae.DerivativeVar(model.Ca, withrespectto=model.t)
model.Cddt = pyodae.DerivativeVar(model.Cd, withrespectto=model.t)
model.M1dt = pyodae.DerivativeVar(model.M1, withrespectto=model.t)
model.M2dt = pyodae.DerivativeVar(model.M2, withrespectto=model.t)
model.P01dt = pyodae.DerivativeVar(model.P01, withrespectto=model.t)
model.P02dt = pyodae.DerivativeVar(model.P02, withrespectto=model.t)
model.P11dt = pyodae.DerivativeVar(model.P11, withrespectto=model.t)
model.P12dt = pyodae.DerivativeVar(model.P12, withrespectto=model.t)
model.P21dt = pyodae.DerivativeVar(model.P21, withrespectto=model.t)
model.P22dt = pyodae.DerivativeVar(model.P22, withrespectto=model.t)
model.D0dt = pyodae.DerivativeVar(model.D0, withrespectto=model.t)
model.D1dt = pyodae.DerivativeVar(model.D1, withrespectto=model.t)
model.D2dt = pyodae.DerivativeVar(model.D2, withrespectto=model.t)
for i in range(1,sites+1,1):
model.CATdt[t_i,i].fix(-_ka[i][0]*CAT0*CoC)
model.Cadt[t_i,i].fix(_ka[i][0]*CAT0*CoC)
model.Cddt[t_i,i].fix(0)
model.M1dt[t_i,i].fix(0)
model.M2dt[t_i,i].fix(0)
model.P01dt[t_i,i].fix(0)
model.P02dt[t_i,i].fix(0)
model.P11dt[t_i,i].fix(0)
model.P12dt[t_i,i].fix(0)
model.P21dt[t_i,i].fix(0)
model.P22dt[t_i,i].fix(0)
model.D0dt[t_i,i].fix(0)
model.D1dt[t_i,i].fix(0)
model.D2dt[t_i,i].fix(0)
model.rule_cat_site = pyo.ConstraintList()
for site in range(1,site+1,1):
for t in model.t:
model.rule_cat_site.add(model.CATdt[t, site] == model.vars[t,"ka", site]*model.CAT[t, site]*model.CoC[site])
model.rule_ca_site = pyo.ConstraintList()
for site in range(1,site+1,1):
for t in model.t:
model.rule_ca_site.add(model.Cadt[t, site] == model.vars[t,"ka", site]*model.CoC[site]*model.CAT[t,site]+(model.vars[t,"ktCoC",site]*model.CoC[site]+model.vars[t,"ktbeta",site]+model.vars[t,"ktH",site]*model.H2[site])*(model.P01[t,site]+model.P02[t,site]) - model.vars[t,"kij1",site]*model.Ca[t,site]*model.M1[t,site] - model.vars[t,"kij2",site]*model.Ca[t,site]*model.M2[t,site] + model.vars[t,"ktij11",site]*model.M1[t,site]*model.P01[t,site] +model.vars[t,"ktij21",site]*model.M1[t,site]*model.P02[t,site] + model.vars[t,"ktij12",site]*model.M2[t,site]*model.P01[t,site] +model.vars[t,"ktij22",site]*model.M2[t,site]*model.P02[t,site])
model.rule_Cd_site = pyo.ConstraintList()
for site in range(1,site+1,1):
for t in model.t:
model.rule_Cd_site.add(model.Cddt[t, site] == model.vars[t,"kd",site]*(model.P01[t,site]+model.P02[t,site]))
sim = pyodae.Simulator(model, package='casadi')
tsim, profiles = sim.simulate(numpoints=int((t_e-t_i)/dt)+1,integrator="idas")#,varying_inputs=model.var_input
Model_8()
我认为问题出在你使用ConstraintList而不是Constraint。pyomo.dae在将模型转换为casadi期望的格式时需要知道哪些约束是时间索引的，ConstraintList对Pyomo隐藏了时间索引集。相反，请尝试这样做:
@model.Constraint(model.t, range(1, site+1))
def rule_cat_site(m, t, site):
return m.CATdt[t, site] == m.vars[t,"ka", site]*m.CAT[t, site]*m.CoC[site]
@model.Constraint(model.t, range(1, site+1))
def rule_ca_site(m, t, site):
return m.Cadt[t, site] == (
m.vars[t,"ka", site]*m.CoC[site]*m.CAT[t,site] + (
m.vars[t,"ktCoC",site]*m.CoC[site] + m.vars[t,"ktbeta",site]
+ m.vars[t,"ktH",site]*m.H2[site])*(m.P01[t,site] + m.P02[t,site]
) - m.vars[t,"kij1",site]*m.Ca[t,site]*m.M1[t,site]
- m.vars[t,"kij2",site]*m.Ca[t,site]*m.M2[t,site]
+ m.vars[t,"ktij11",site]*m.M1[t,site]*m.P01[t,site]
+ m.vars[t,"ktij21",site]*m.M1[t,site]*m.P02[t,site]
+ m.vars[t,"ktij12",site]*m.M2[t,site]*m.P01[t,site]
+ m.vars[t,"ktij22",site]*m.M2[t,site]*m.P02[t,site]
)
@model.Constraint(model.t, range(1, site+1))
def rule_Cd_site(m, t, site):
return m.Cddt[t, site] == m.vars[t,"kd",site]*(m.P01[t,site] + m.P02[t,site])