Solve DAE with Pyomo and class

Question

I'm trying to solve a car problem. first, I have an original code of car problem:

# Ampl Car Example
#
# Shows how to convert a minimize final time optimal control problem
# to a format pyomo.dae can handle by removing the time scaling from
# the ContinuousSet.
#
# min tf
# dxdt = 0
# dvdt = a-R*v^2
# x(0)=0; x(tf)=L
# v(0)=0; v(tf)=0
# -3<=a<=1

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

m = ConcreteModel()

m.R = Param(initialize=0.001) #  Friction factor
m.L = Param(initialize=100.0) #  Final position

m.tau = ContinuousSet(initialize=[0.0, 0.80, 1.0])  # Unscaled time
m.time = Var(m.tau)  # Scaled time
m.tf = Var()
m.x = Var(m.tau,bounds=(0,None))
m.v = Var(m.tau,bounds=(0,None))
m.a = Var(m.tau, bounds=(-3.0,1.0),initialize=0)

m.dtime = DerivativeVar(m.time)
m.dx = DerivativeVar(m.x)
m.dv = DerivativeVar(m.v)

m.obj = Objective(expr=m.tf)

def _ode1(m,i):
    if i == 0 :
        return Constraint.Skip
    return m.dx[i] == m.tf * m.v[i]
m.ode1 = Constraint(m.tau, rule=_ode1)

def _ode2(m,i):
    if i == 0 :
        return Constraint.Skip
    return m.dv[i] == m.tf*(m.a[i] - m.R*m.v[i]**2)
m.ode2 = Constraint(m.tau, rule=_ode2)

def _ode3(m,i):
    if i == 0:
        return Constraint.Skip
    return m.dtime[i] == m.tf
m.ode3 = Constraint(m.tau, rule=_ode3)

def _init(m):
    yield m.x[0] == 0
    yield m.x[1] == m.L
    yield m.v[0] == 0
    yield m.v[1] == 0
    yield m.time[0] == 0
m.initcon = ConstraintList(rule=_init)

discretizer = TransformationFactory('dae.collocation')
discretizer.apply_to(m,ncp=1, scheme='LAGRANGE-RADAU')

solver = SolverFactory('ipopt')
solver.solve(m, tee=True)

print("final time = %6.2f" %(value(m.tf)))

Now, I want to use class to express a car,then I could instantiate two cars. So I write like this:

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

m = ConcreteModel()

class Car():
    def __init__(self,friction):

        self.friction = friction

        self.R = Param(initialize = self.friction)  # Friction factor

        self.tau = ContinuousSet(bounds=(0, 1))  # Unscaled time
        self.time = Var(self.tau)  # Scaled time
        self.tf = Var()
        self.x = Var(self.tau, bounds=(0, None), initialize=0)
        self.v = Var(self.tau, bounds=(0, None))
        self.a = Var(self.tau, bounds=(-3.0, 1.0), initialize=0)

        self.dtime = DerivativeVar(self.time)
        self.dx = DerivativeVar(self.x)
        self.dv = DerivativeVar(self.v)

        def _ode1(m, i):
            if i == 0:
                return Constraint.Skip
            return self.dx[i] == m.tf * self.v[i]
        self.ode1 = Constraint(self.tau, rule=_ode1)

        def _ode2(m, i):
            if i == 0:
                return Constraint.Skip
           return self.dv[i] == m.tf * (self.a[i] - self.R * self.v[i] ** 2)

        self.ode2 = Constraint(self.tau, rule=_ode2)

        def _ode3(m, i):
            if i == 0:
                return Constraint.Skip
            return self.dtime[i] == m.tf

        self.ode3 = Constraint(self.tau, rule=_ode3)

m.car1 = Car(0.001)

m.obj = Objective(expr=m.car1.tf)

def _init(m):
    yield m.car1.x[0] == 0
    yield m.car1.x[1] == 100
    yield m.car1.v[0] == 0
    yield m.car1.v[1] == 0
    yield m.car1.time[0] == 0
m.car1.initcon = ConstraintList(rule=_init)

discretizer = TransformationFactory('dae.finite_difference')
discretizer.apply_to(m, nfe=10, scheme='BACKWARD')

solver = SolverFactory('ipopt')
solver.solve(m, tee=True)

print("final time = %6.2f" % (value(m.car1.tf)))

However, I get this:

Traceback (most recent call last):
File "D:/pyo/pyomoceshi/ceshi3/car/classcar3.py", line 79, in <module>
solver.solve(m, tee=True)
File "D:\python\m\lib\site-packages\pyomo\opt\base\solvers.py", line 582, in solve
self._presolve(*args, **kwds)
File "D:\python\m\lib\site-packages\pyomo\opt\solver\shellcmd.py", line 196, in _presolve
OptSolver._presolve(self, *args, **kwds)
File "D:\python\m\lib\site-packages\pyomo\opt\base\solvers.py", line 661, in _presolve
**kwds)
File "D:\python\m\lib\site-packages\pyomo\opt\base\solvers.py", line 729, in _convert_problem
**kwds)
File "D:\python\m\lib\site-packages\pyomo\opt\base\convert.py", line 110, in convert_problem
problem_files, symbol_map = converter.apply(*tmp, **tmpkw)
File "D:\python\m\lib\site-packages\pyomo\solvers\plugins\converter\model.py", line 164, in apply
io_options=io_options)
File "D:\python\m\lib\site-packages\pyomo\core\base\block.py", line 1646, in write
io_options)
File "D:\python\m\lib\site-packages\pyomo\repn\plugins\ampl\ampl_.py", line 357, in __call__
include_all_variable_bounds=include_all_variable_bounds)
File "D:\python\m\lib\site-packages\pyomo\repn\plugins\ampl\ampl_.py", line 783, in _print_model_NL
list(self_varID_map[id(var)] for var in ampl_repn._linear_vars),
File "D:\python\m\lib\site-packages\pyomo\repn\plugins\ampl\ampl_.py", line 783, in <genexpr>
list(self_varID_map[id(var)] for var in ampl_repn._linear_vars),
KeyError: 68767416L

I want to know how to solve it or use other ways.

Answer 1

Below is a working version of your script. I changed things so that instead of a Car class there is a Car function that returns a Pyomo Block representing the car. By having a Car class you were essentially trying to create a subclass of Block and running into several subtle challenges that go along with that. You can see the blog post here for more information. The second change I made was in your declaration of the initial conditions, I changed the name of the ConstraintList from m.car1.initcon to m.car1_initcon. The difference is whether you want the ConstraintList to live on the car1 Block or the model. In your code, the 'dot' in the name meant you were trying to put it on the car1 Block but the constraints yielded in the rule were relative to the model. I changed the name to resolve this inconsistency.

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

m = ConcreteModel()

def Car(model, friction):

    def construct_car_block(b):

        b.R = Param(initialize = friction)  # Friction factor
        b.tau = ContinuousSet(bounds=(0, 1))  # Unscaled time
        b.time = Var(b.tau)  # Scaled time
        b.tf = Var()
        b.x = Var(b.tau, bounds=(0, None), initialize=0)
        b.v = Var(b.tau, bounds=(0, None))
        b.a = Var(b.tau, bounds=(-3.0, 1.0), initialize=0)

        b.dtime = DerivativeVar(b.time)
        b.dx = DerivativeVar(b.x)
        b.dv = DerivativeVar(b.v)

        def _ode1(b, i):
            if i == 0:
                return Constraint.Skip
            return b.dx[i] == b.tf * b.v[i]
        b.ode1 = Constraint(b.tau, rule=_ode1)

        def _ode2(b, i):
            if i == 0:
                return Constraint.Skip
            return b.dv[i] == b.tf * (b.a[i] - b.R * b.v[i] ** 2)
        b.ode2 = Constraint(b.tau, rule=_ode2)

        def _ode3(m, i):
            if i == 0:
                return Constraint.Skip
            return b.dtime[i] == b.tf
        b.ode3 = Constraint(b.tau, rule=_ode3)

    return Block(rule=construct_car_block)

m.car1 = Car(m, friction=0.001)

m.obj = Objective(expr=m.car1.tf)

def _init(m):
    yield m.car1.x[0] == 0
    yield m.car1.x[1] == 100
    yield m.car1.v[0] == 0
    yield m.car1.v[1] == 0
    yield m.car1.time[0] == 0
m.car1_initcon = ConstraintList(rule=_init)

discretizer = TransformationFactory('dae.finite_difference')
discretizer.apply_to(m, nfe=15, scheme='BACKWARD')

solver = SolverFactory('ipopt')
solver.solve(m, tee=True)

print("final time = %6.2f" % (value(m.car1.tf)))