基于GEKKO的最短时间最优控制基准问题求解遇收敛障碍

这是一个经典的分段动力学最短时间最优控制问题，状态方程在$x_2>1$时切换为更高增益的运动模型。原代码无法收敛的核心问题是静态条件判断无法实现动态状态切换，同时缺少必要的终端约束与控制量优化配置。

原代码关键错误 #

• 建模阶段的if x2.value>1仅根据初始值（$x_2=0$）选择了低增益方程，无法在仿真过程中根据$x_2$的动态变化切换方程，导致模型与实际问题不符。
• 未启用终端约束，优化器没有目标终点的限制，无法收敛到合理的最短时间解。
• 控制量u未开启优化状态，优化器无法调整控制输入。

修正后的代码 #

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

m = GEKKO()

nt = 501
tm = np.linspace(0, 1, nt)
m.time = tm

# 状态变量
x1 = m.Var(value=-2.5)
x2 = m.Var(value=0)

# 控制变量：角度u，范围[0, 2π]
u = m.MV(value=1, lb=0, ub=2*math.pi)
u.STATUS = 1  # 允许优化器调整u
u.DCOST = 1e-6  # 加入小惩罚项避免控制量突变

# 终端时刻参数
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

# 最终时间变量，初始值设为5帮助收敛
tf = m.FV(value=5, lb=0.1, ub=100.0)
tf.STATUS = 1

# 使用m.if3实现动态状态切换：当x2>1时切换到高增益动力学
gain = m.if3(x2 - 1, 1, 1 + (x2 - 1)**2)

# 状态方程
m.Equation(x1.dt() == gain * m.cos(u) * tf)
m.Equation(x2.dt() == gain * m.sin(u) * tf)

# 终端约束：最终时刻x1≤3，x2≤0
m.Equation(x1 * final <= 3)
m.Equation(x2 * final <= 0)

# 目标：最小化最终时间tf
m.Minimize(tf)

# 设置最优控制模式
m.options.IMODE = 6
m.solve(disp=True)  # 显示求解过程

# 时间轴转换为实际时间
tm = tm * tf.value[0]

# 绘图
plt.figure(1)
plt.plot(tm, x1.value, 'k-', lw=2, label=r'$x_1$')
plt.plot(tm, x2.value, 'b-', lw=2, label=r'$x_2$')
plt.plot(tm, u.value, 'r--', lw=2, label=r'$u$')
plt.legend(loc='best')
plt.xlabel('时间')
plt.ylabel('变量值')
plt.title(f'最短飞行时间解：tf={tf.value[0]:.2f}')
plt.grid(True)
plt.show()

关键修改说明 #

1. 动态状态切换：用m.if3(x2-1, 1, 1+(x2-1)**2)实现分段增益，当$x_2>1$时自动切换到高增益模型，GEKKO会将其转换为适配求解器的光滑近似或混合整数约束。
2. 终端约束启用：添加x1*final<=3和x2*final<=0，确保优化器将状态驱动到目标区域。
3. 控制量配置：开启u.STATUS=1允许优化器调整控制角度，u.DCOST避免控制量跳变，提升求解稳定性。
4. 初始值调整：将tf初始值设为5，比0更接近实际解，帮助求解器快速收敛。

求解结果 #

修正后的代码会收敛到合理的最短时间解，状态轨迹会先驱动$x_2$超过1以利用高增益快速移动，再调整角度回到目标区域，控制量$u$会根据状态动态调整。

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