0

I've the following code which uses python pyomo. I think my code is right, but it can't run. I don't understand what's wrong. False statement casadi/core/integrator.cpp:1417: Assertion "!de_in[DE_X].is_empty()" failed:Ill-posed ODE - no state。

this is my code,and for simple, I just put part of code:

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()

i hope deal with the problem. thanks very much

jsiirola
  • 2,259
  • 2
  • 9
  • 12
shuyuan fan
  • 1

1 Answers1

0

I believe the problem is in your use of ConstraintList instead of Constraint. pyomo.dae needs to know which constraints are time-indexed when it is converting the model to the format expected by casadi, and the ConstraintList hides the time indexing set from Pyomo. Instead, try something like:

@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])
jsiirola
  • 2,259
  • 2
  • 9
  • 12