用 Pyomo 和类解决 DAE

Question

我正在尝试解决汽车问题。 首先，我有一个汽车问题的原始代码：

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

现在，我想用类来表示一辆汽车，然后我可以实例化两辆汽车。 所以我是这样写的：

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

但是，我明白了：

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

我想知道如何解决或使用其他方法。

Answer 1

以下是您的脚本的工作版本。我改变了一些东西，而不是 Car 类，而是一个 Car 函数，它返回一个代表汽车的 Pyomo 块。通过拥有 Car 类，您实际上是在尝试创建 Block 的子类，并遇到随之而来的几个微妙挑战。您可以查看博文here 了解更多信息。我所做的第二个更改是在您声明初始条件时，我将 ConstraintList 的名称从 m.car1.initcon 更改为 m.car1_initcon。区别在于您是否希望 ConstraintList 存在于 car1 块或模型上。在您的代码中，名称中的“点”表示您试图将其放在 car1 块上，但规则中产生的约束与模型相关。我更改了名称以解决此不一致问题。

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