Scarlett Orangin
Scarlett Orangin

Reputation: 1

Pyomo: Getting keyerror: 4903576400

I'm trying to solve an optimal control problem using pyomo, but I always get this keyerror which seems like a series of random numbers.

KeyError                                  Traceback (most recent call last)
<ipython-input-28-926e5210dc1f> in <module>
     69 
     70 solver = SolverFactory('ipopt')
---> 71 res = solver.solve(m)

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/pyomo/opt/base/solvers.py in solve(self, *args, **kwds)
    568             initial_time = time.time()
    569 
--> 570             self._presolve(*args, **kwds)
    571 
    572             presolve_completion_time = time.time()

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/pyomo/opt/solver/shellcmd.py in _presolve(self, *args, **kwds)
    207         self._define_signal_handlers = kwds.pop('use_signal_handling',None)
    208 
--> 209         OptSolver._presolve(self, *args, **kwds)
    210 
    211         #

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/pyomo/opt/base/solvers.py in _presolve(self, *args, **kwds)
    665             write_start_time = time.time()
    666             (self._problem_files, self._problem_format, self._smap_id) = \
--> 667                 self._convert_problem(args,
...
pyomo/repn/plugins/ampl/ampl_.pyx in genexpr()

pyomo/repn/plugins/ampl/ampl_.pyx in genexpr()

KeyError: 4903576400

Here is my code, my objective function is an integral, and the constraints are 4 ode, I am trying to use 'ipopt' to solve it. I don't really know what's wrong with it...

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

# parameter
r1 = 1.3195 * 10**(-6)
r2 = 1.885 * 10**(-5)
r3 = 0.19
r4 = 0.2

k = 0.00054

# initial value
S0 = 102028
E0 = 8251
I0 = 1843
R0 = 40435

m = ConcreteModel()

m.t = ContinuousSet(bounds=(0,1000))
m.S = Var(m.t,initialize = S0)
m.E = Var(m.t,initialize = E0)
m.I = Var(m.t,initialize = I0)
m.R = Var(m.t,initialize = R0)
m.u = Var(m.t)

m.dSdt = DerivativeVar(m.S,wrt=m.t)
m.dEdt = DerivativeVar(m.E,wrt=m.t)
m.dIdt = DerivativeVar(m.I,wrt=m.t)
m.dRdt = DerivativeVar(m.R,wrt=m.t)

def ode1(m,t):
    if t == 0:
        return Constraint.Skip
    return m.dSdt[t] == -r1*m.S[t]*(1-m.u[t])*m.E[t] - r2*m.S[t]*(1-m.u[t])*m.I[t]
m.odeS = Constraint(m.t,rule=ode1)

def ode2(m,t):
    if t == 0:
        return Constraint.Skip
    return m.dEdt[t] == r1*m.S[t]*(1-m.u[t])*m.E[t] + r2*m.S[t]*(1-m.u[t])*m.I[t] - r3*m.E[t]
m.odeE = Constraint(m.t,rule=ode2)

def ode3(m,t):
    if t == 0:
        return Constraint.Skip
    return m.dIdt[t] == r3*m.E[t] - r4*m.I[t]
m.odeI = Constraint(m.t,rule=ode3)

def ode4(m,t):
    if t == 0:
        return Constraint.Skip
    return m.dRdt[t] == r4*m.I[t]
m.odeR = Constraint(m.t,rule=ode4)


def myintegral(model,t):
    return k * model.I[t] + model.u[t]**2

m.o = Integral(m.t, wrt=m.t, rule=myintegral) 

def myintegral(m,t):
    return m.o
 
m.obj = Objective(rule=myintegral,sense=minimize)
m.pprint()

solver = SolverFactory('ipopt')
res = solver.solve(m)

I learned all the things by myself but was actually not that thorough, and I really need your help...

Upvotes: 0

Views: 80

Answers (1)

jsiirola
jsiirola

Reputation: 2599

The KeyError is coming from the Pyomo "NL writer" (it translates your model into an AMPL NL file that the ipopt solver can read and optimize. The specific cause is that your model has constraints that reference objects other than Vars, Params and native constants (in particular, you have DerivativeVar and Integral objects). This is admittedly a poor error message, and has been resolved by a rewrite of the NL writer (available in Pyomo >= 6.4.2, but will not become the default until 6.5).

You can resolve the error by transforming your model from a DAE into an NLP before sending it to a solver. There are several transformations available; see the DAE extension documentation for more information and examples.

Upvotes: 1

Related Questions