使用Pyomo求解最优控制问题时出现KeyError:4903576400求助

问题描述 #

使用Pyomo求解含积分目标函数和4个常微分方程约束的最优控制问题，采用IPOPT求解器时出现KeyError:4903576400报错，无法定位问题。

报错堆栈 #

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

原代码 #

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)

错误分析与修正方案 #

核心错误点 #

1. 重复定义函数覆盖积分规则：后续定义的myintegral函数覆盖了之前用于生成积分项的函数，导致Objective的规则逻辑混乱。
2. Objective规则错误：Objective的rule返回的是整个Integral对象m.o，但Pyomo要求Objective直接接收标量表达式，无需额外规则包装。
3. 缺失初始条件约束：虽然状态变量设置了初始值，但t=0时的ODE被Constraint.Skip跳过，需要明确添加初始值约束确保问题边界条件完整。
4. 未离散化DAE系统：Pyomo的DAE问题必须先离散化才能使用IPOPT等NLP求解器，原代码缺少离散化步骤。
5. 控制变量无界：控制变量u未设置合理上下界，可能导致求解器数值不稳定。

修正后的代码 #

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

# 参数定义
r1 = 1.3195 * 10**(-6)
r2 = 1.885 * 10**(-5)
r3 = 0.19
r4 = 0.2

k = 0.00054

# 初始值
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, bounds=(0, 1))  # 添加控制变量合理边界

# 导数变量
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)

# ODE约束（无需跳过t=0，Pyomo会自动处理初始点的导数约束）
def ode1(m,t):
    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):
    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):
    return m.dIdt[t] == r3*m.E[t] - r4*m.I[t]
m.odeI = Constraint(m.t, rule=ode3)

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

# 初始条件约束
def init_S(m):
    return m.S[0] == S0
m.init_S = Constraint(rule=init_S)

def init_E(m):
    return m.E[0] == E0
m.init_E = Constraint(rule=init_E)

def init_I(m):
    return m.I[0] == I0
m.init_I = Constraint(rule=init_I)

def init_R(m):
    return m.R[0] == R0
m.init_R = Constraint(rule=init_R)

# 积分目标函数
def integral_rule(model,t):
    return k * model.I[t] + model.u[t]**2

m.obj_integral = Integral(m.t, wrt=m.t, rule=integral_rule) 

# 目标函数直接使用积分项
m.obj = Objective(expr=m.obj_integral, sense=minimize)

# 离散化DAE系统（必须步骤）
discretizer = TransformationFactory('dae.collocation')
discretizer.apply_to(m, nfe=50, ncp=3)  # 可调整nfe和ncp平衡精度与速度

# 求解
solver = SolverFactory('ipopt')
res = solver.solve(m, tee=True)

# 打印求解结果
print("求解状态:", res.solver.status)
print("求解终止条件:", res.solver.termination_condition)

关键修正说明 #

• 移除重复函数定义：重新命名积分规则函数避免覆盖，Objective直接通过expr参数使用积分项，无需额外规则包装。
• 添加初始条件约束：明确设置t=0时各状态变量的初始值，确保问题边界条件完整。
• 离散化DAE系统：采用配置法离散化连续时间系统，这是Pyomo求解DAE问题的必要步骤，可通过调整nfe（有限元数量）和ncp（每个单元的配置点数）优化求解精度和速度。
• 控制变量加边界：设置u ∈ [0,1]符合问题中(1-u)的物理意义，避免求解器出现数值发散。

内容的提问来源于stack exchange，提问作者Scarlett Orangin