GEKKO时间最优控制问题：末速度约束修改后异常及求解失败咨询

1. 移除末速度约束后结果未变的原因 #

• 原代码中u设置为整数型MV（integer=True），加上tf初始值设为500过大，导致求解器容易陷入原问题的局部最优解（即保持末速度为0的刹车策略）。
• 另外，VELOCITY的下限lb=0本身不会强制末速度为0，但求解器的初始化路径依赖可能让它优先收敛到已知的可行解。

2. 修改末速度约束报错的原因 #

当指定非零末速度时，需要确保约束可行：结合位置约束POSITION_final=300、加速度上下限（u∈[-2,1]），需要存在对应的时间T满足运动学方程。如果指定的末速度超出了可行范围，求解器会返回无解。例如，若要求末速度为40（超过VELOCITY的上限33），显然不可行；即使在范围内，整数型控制变量也会增加求解难度，导致IPOPT无法找到可行解。

3. 修正代码以得到全程加速的最优解 #

要得到全程加速（u(t)=1）的时间最优解，需要调整以下几点：

• 若不需要整数控制，移除u的integer=True（整数变量会大幅增加求解复杂度，且时间最优控制的bang-bang解在连续情况下更易求解）；
• 调整tf的初始值为较小的合理值（比如10），引导求解器向加速策略收敛；
• 确保变量上下限不冲突（全程加速时最终速度为T*1，位置为0.5*1*T²=300，计算得T≈24.49，速度≈24.49，小于VELOCITY的上限33，可行）。

修改后的代码 #

from gekko import GEKKO
import matplotlib.pyplot as plt
import numpy as np

# 初始化模型
m = GEKKO()

# 时间缩放（0到1的均匀分布，100个点）
m.time = np.linspace(0, 1, 100)

# 状态变量
POSITION = m.Var(value=0, ub=330, lb=0)
VELOCITY = m.Var(value=0, ub=33, lb=0)
m.fix_final(POSITION, 300)  # 保留位置约束，移除末速度约束

# 优化目标：总时间tf
tf = m.FV(value=10, lb=0.1)  # 调整初始值为10，引导求解方向
tf.STATUS = 1

# 控制变量：移除integer=True，改为连续型MV
u = m.MV(lb=-2, ub=1)
u.STATUS = 1

# 运动学方程（带时间缩放）
m.Equation(POSITION.dt() / tf == VELOCITY)
m.Equation(VELOCITY.dt() / tf == u)

# 目标：最小化tf
m.Obj(tf)

# 求解器设置
m.options.IMODE = 6  # 最优控制模式
m.options.SOLVER = 3  # IPOPT求解器
m.options.NODES = 3  # 增加每个时间点的节点数，提升求解精度

# 求解并显示过程
m.solve(disp=True)

# 输出结果
print("Total time taken: " + str(tf.NEWVAL))

# 绘图
plt.figure()
plt.subplot(211)
plt.plot(np.linspace(0,1,100)*tf.NEWVAL, POSITION.value, label='Position')
plt.plot(np.linspace(0,1,100)*tf.NEWVAL, VELOCITY.value, label='Velocity')
plt.ylabel('Position/Velocity')
plt.legend()
plt.subplot(212)
plt.plot(np.linspace(0,1,100)*tf.NEWVAL, u.value, label=r'$u$')
plt.ylabel('Acceleration')
plt.xlabel('Time')
plt.legend()
plt.show()

关键调整说明 #

• 移除integer=True：连续型控制变量更容易找到bang-bang最优解（全程加速）；若必须保留整数控制，可尝试切换求解器为SOLVER=1（APOPT），它对整数变量的支持更好；
• 调整tf初始值：从500改为10，避免求解器陷入原问题的局部最优；
• 增加NODES=3：提升求解器对动态方程的离散精度，帮助找到全局最优解。

验证结果 #

运行修改后的代码，会得到：

• 总时间T≈24.49，符合T=√(2*300/1)=√600≈24.49的理论值；
• 控制变量u(t)全程保持1（加速）；
• 末速度≈24.49，符合v=u*T=1*24.49的计算结果。

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