1

Pyomo リポジトリで「ampl car example」の素晴らしい実装を見た後、新しい機能と制約で問題を拡張し続けたいと思っていますが、開発中に次の問題を発見しました。誰かがそれらを修正できますか？

1) 新しい制約「電気自動車」が追加されました: 現在、加速度は一定の速度になるまで順守することによって制限され、一定の電力モデルが使用されます。私が思うように、この制約を実装することはできません。でコメントされていますが、ピョーモは制約が変数に関連していると文句を言っています。(現在、Umax は車速に依存します)。

2) 新しいコンフォート アクセラレーションとジャークの制約が追加されました。それらは正しく機能しているように見えますが、ピョーモの第一人者がそれらを監督し、それらが本当に正しい方法で実装されているかどうかを教えてくれればいいと思います.

3) 冗長性を減らすために、最後のものについて。一意の制約でaccelerationLとaccelerationUを組み合わせる方法はありますか? jerkL と jerkU も同様です。

4) 最後の特徴は、2 段階に分かれた速度制限制約です。繰り返しますが、私はそれを機能させることができないので、コードでコメントされています。誰かがそれを修正することを敢えてしますか？

# 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 制約ルールを、ソルバーによって使用されるコールバックと考えるべきではありません。モデルの構築時にインデックスごとに 1 回呼び出される、制約オブジェクトのコンテナーを生成する関数と考える必要があります。ifつまり、変数の初期値のみを使用して制約式を定義していない限り、ステートメントで変数を使用することは無効です。あなたがやろうとしていることを表現する方法はいくつかありますが、バイナリ変数を問題に導入する必要があるため、Ipopt を使用できなくなります。

(2) 本当に何の助けにもならない. 構文は問題ないようです。

(3) Pyomo では、制約ルールから両側の不等式 (例: L <= f(x) <= U) を返すことができますが、L と U の場所に変数式を含めることはできません。あなたが参照している制約をこの形式に組み合わせることができるようには見えません。

(4) (1)参照

于 2016-10-11T18:17:06.157 に答える