Pyomo+Ipopt扩展DAE优化问题时无法找到最优解的咨询

我使用Pyomo搭配Ipopt求解器时，发现原本可收敛至最优解的DAE优化问题，在提升问题复杂度（如汽车示例中增大行驶距离）并为保证精度增加有限元数量后，求解器输出：

EXIT: Solved To Acceptable Level.
而非此前的“Optimal solution found”。

修改自Pyomo官方汽车示例的代码如下：

# 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 = v
# 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=1000000.0) #  Final position

m.tau = ContinuousSet(bounds=(0,1)) # Unscaled time
m.time = Var(m.tau) # Scaled time
m.tf = Var()
m.x = Var(m.tau,bounds=(0,m.L+50))
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.finite_difference')
discretizer.apply_to(m,nfe=5000,scheme='BACKWARD')

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

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

x = []
v = []
a = []
time=[]

for i in m.tau:
    time.append(value(m.time[i]))
    x.append(value(m.x[i]))
    v.append(value(m.v[i]))
    a.append(value(m.a[i]))
  
import matplotlib.pyplot as plt

plt.subplot(131)
plt.plot(time,x,label='x')
plt.title('location')
plt.xlabel('time')

plt.subplot(132)
plt.plot(time,v,label='v')
plt.xlabel('time')
plt.title('velocity')

plt.subplot(133)
plt.plot(time,a,label='a')
plt.xlabel('time')
plt.title('acceleration')

plt.show()

可行解决措施 #

1. 调整Ipopt收敛容差（可直接实现） #

是的，你可以通过降低Ipopt的收敛容差，迫使求解器继续寻找更优解。具体做法是在调用solve()时传入选项参数，修改核心容差阈值：

solver = SolverFactory('ipopt')
solver.solve(m, tee=True, options={
    'tol': 1e-10,                # 主最优性容差，默认1e-8
    'constr_viol_tol': 1e-10,    # 约束违反容差，默认1e-6
    'acceptable_tol': 1e-12,     # 可接受解的最优性容差，默认1e-6
    'acceptable_constr_viol_tol': 1e-12  # 可接受解的约束违反容差，默认1e-4
})

注意：容差设置过小会显著增加求解时间，甚至可能导致求解器无法收敛，需要根据问题需求逐步调整。

2. 优化离散化策略 #

• 替换离散化方法：放弃固定节点数的有限差分法，改用正交配置法（collocation）。配置法在相同精度下所需节点数更少，能有效降低问题规模，提升数值稳定性：

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

• 自适应离散化：使用Pyomo的自适应细化功能，只在误差较大的区域增加节点，避免盲目增大NFE导致的问题爆炸：

from pyomo.dae.plugins.finitedifference import FiniteDifference
discretizer = FiniteDifference()
discretizer.apply_to(m, refinement_method='adaptive', max_refinement=5)

3. 变量与约束缩放优化 #

问题中行驶距离L=1e6，变量量级差异过大可能导致数值条件数恶化，影响Ipopt的收敛性。可以对变量进行无量纲化处理：

• 例如将位置变量x除以L，转化为x_norm = x/L，对应的约束和目标函数也做相应调整，缩小变量间的量级差距。

4. 调整Ipopt其他求解参数 #

• 增大最大迭代次数：避免因迭代次数不足提前终止，添加max_iter=10000选项。
• 更换线性求解器：默认的mumps在大规模问题下效率较低，尝试使用linear_solver='ma57'或ma97（需提前安装对应线性代数库），提升线性系统求解的稳定性。
• 启用热启动：如果有相近问题的最优解，将其作为初始值传入模型，帮助Ipopt更快定位到最优解区域。

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