1

在看到 Pyomo 存储库中“ampl car example”的良好实现之后，我想继续用新的特性和约束来扩展这个问题，但是我在开发过程中发现了下一个问题。有人能够修复它们吗？

1) 增加了新的约束“电动车”：现在通过坚持限制加速度，直到确定速度，然后使用恒功率模型。我无法像我想的那样实施这个约束。中对此进行了评论，但 Pyomo 抱怨约束与变量有关。（现在 Umax 取决于车速）。

2) 增加了新的舒适加速度和加加速度限制。看起来他们工作正常，但如果 Pyomo 大师监督他们并告诉我他们是否真的以正确的方式实施，那应该很好。

3）关于最后一个，以减少冗长的顺序。有什么方法可以在一个唯一约束中组合加速 L 和加速 U 吗？jerkL 和 jerkU 也一样。

4）最后一个特征是速度限制约束，分为两步。同样，我无法让它工作，所以它在代码中进行了注释。有人敢修吗？

# Ampl Car Example (Extended)
#
# 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
# dx/dt = v
# dv/dt = u - R*v^2
# x(0)=0; x(tf)=L
# v(0)=0; v(tf)=0
# -3 <= u <= 1 (engine constraint)
#
#      {v <= 7m/s ===> u < 1
# u <= {                            (electric car constraint)
#      {v >  7m/s ===> u < 1*7/v
#
# -1.5 <= dv/dt <= 0.8 (comfort constraint -> smooth driving)
# -0.5 <= d2v/dt2 <= 0.5 (comfort constraint -> jerk)
# v <= Vmax (40 kmh[0-500m] + 25 kmh(500-1000m])

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

m = ConcreteModel()

m.R = Param(initialize=0.001)  # Friction factor
m.L = Param(initialize=1000.0) # Final position
m.T = Param(initialize=50.0)   # Estimated time
m.aU = Param(initialize=0.8)   # Acceleration upper bound
m.aL = Param(initialize=-1.5)  # Acceleration lower bound
m.jU = Param(initialize=0.5)   # Jerk upper bound
m.jL = Param(initialize=-0.5)  # Jerk lower bound
m.NFE = Param(initialize=100)   # Number of finite elements

'''
def _initX(m, i):
    return m.x[i] == i*m.L/m.NFE

def _initV(m):
    return m.v[i] == m.L/50
'''

m.tf = Var()
m.tau = ContinuousSet(bounds=(0,1)) # Unscaled time
m.t = Var(m.tau) # Scaled time
m.x = Var(m.tau, bounds=(0,m.L))
m.v = Var(m.tau, bounds=(0,None))
m.u = Var(m.tau, bounds=(-3,1), initialize=0)

m.dt = DerivativeVar(m.t)
m.dx = DerivativeVar(m.x)
m.dv = DerivativeVar(m.v)
m.da = DerivativeVar(m.v, wrt=(m.tau, m.tau))

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

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

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

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

def _accelerationL(m, i):
    if i==0:
        return Constraint.Skip
    return m.dv[i] >= m.aL*m.tf
m.accelerationL = Constraint(m.tau, rule=_accelerationL)

def _accelerationU(m, i):
    if i==0:
        return Constraint.Skip
    return m.dv[i] <= m.aU*m.tf
m.accelerationU = Constraint(m.tau, rule=_accelerationU)

def _jerkL(m, i):
    if i==0:
        return Constraint.Skip
    return m.da[i] >= m.jL*m.tf**2
m.jerkL = Constraint(m.tau, rule=_jerkL)

def _jerkU(m, i):
    if i==0:
        return Constraint.Skip
    return m.da[i] <= m.jU*m.tf**2
m.jerkU = Constraint(m.tau, rule=_jerkU)

'''
def _electric(m, i):
    if i==0:
        return Constraint.Skip
    elif value(m.v[i])<=7:
        return m.a[i] <= 1
    else:
        return m.v[i] <= 1*7/m.v[i]
m.electric = Constraint(m.tau, rule=_electric)
'''

'''
def _speed(m, i):
    if i==0:
        return Constraint.Skip
    elif value(m.x[i])<=500:
        return m.v[i] <= 40/3.6
    else:
        return m.v[i] <= 25/3.6
m.speed = Constraint(m.tau, rule=_speed)
'''

def _initial(m):
    yield m.x[0] == 0
    yield m.x[1] == m.L
    yield m.v[0] == 0
    yield m.v[1] == 0
    yield m.t[0] == 0
m.initial = ConstraintList(rule=_initial)

discretizer = TransformationFactory('dae.finite_difference')
discretizer.apply_to(m, nfe=value(m.NFE), wrt=m.tau, scheme='BACKWARD')
#discretizer = TransformationFactory('dae.collocation')
#discretizer.apply_to(m, nfe=value(m.NFE), ncp=4, wrt=m.tau, scheme='LAGRANGE-RADAU')

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

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

t = []
x = []
v = []
a = []
u = []

for i in m.tau:
    t.append(value(m.t[i]))
    x.append(value(m.x[i]))
    v.append(3.6*value(m.v[i]))
    a.append(10*value(m.u[i] - m.R*m.v[i]**2))
    u.append(10*value(m.u[i]))

import matplotlib.pyplot as plt

plt.plot(x, v, label='v (km/h)')
plt.plot(x, a, label='a (dm/s2)')
plt.plot(x, u, label='u (dm/s2)')
plt.xlabel('distance')
plt.grid('on')
plt.legend()

plt.show()

非常感谢，巴勃罗

4

1 回答 1

0

(1) 您不应将 Pyomo 约束规则视为求解器使用的回调。您应该将它们更多地视为一个函数，用于生成约束对象的容器，在构建模型时为每个索引调用一次。这意味着在语句中使用变量是无效的，if除非你真的只使用它的初始值来定义约束表达式。有一些方法可以表达我认为你正在尝试做的事情，但它们涉及将二进制变量引入问题，在这种情况下你不能再使用 Ipopt。

(2) 不能真正提供任何帮助。语法看起来不错。

(3) Pyomo 允许您从约束规则返回双边不等式表达式（例如，L <= f(x) <= U），但它们不能涉及 L 和 U 位置的变量表达式。看起来您所指的约束不能组合成这种形式。

(4) 见 (1)

于 2016-10-11T18:17:06.157 回答