GEKKO求解最大半径轨道转移OCP的终端速度约束添加方法

核心问题 #

原代码无法正确添加终端圆轨道速度约束v(tf)=√(μ/r(tf))，原因是m.fix()仅支持固定常数值约束，无法传入依赖其他优化变量的表达式；模型构建阶段r[nt-1]无求解值，直接传参会触发报错。
此外原代码缺少推力幅值固定的物理约束，会导致求解结果不符合问题设定。

正确实现方法 #

1. 终端速度约束实现 #

利用已定义的final参数（仅终端时刻值为1，其余时刻为0），通过乘子法实现仅在终端时刻生效的等式约束，采用平方形式避免开方运算，提升数值求解稳定性：

# 替换原注释掉的m.fix(v,...)语句，等价于终端时刻v=sqrt(μ/r)
m.Equation(final * (v**2 - mu/r) == 0)

2. 补全推力物理约束 #

问题中推力幅值固定为thr、仅方向可调，原代码将径向、切向推力分量u1/u2设为独立操纵变量，缺少分量平方和为1的约束，会导致求解器默认两个方向同时满推力输出，结果不符合物理规律，需补充全时域约束：

m.Equation(u1**2 + u2**2 == 1)

3. 初值修正 #

原代码中切向推力分量u2初值设为np.pi，会导致初始推力方向与速度方向反向，不利于收敛，调整为初值1（初始推力沿切向正方向）即可。

完整修正后可运行代码 #

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


if __name__ == '__main__':

    m = GEKKO()

    # 常量定义
    nt = 101
    m.time = np.linspace(0, 193, nt)
    thr = 0.85
    mu = (3.32/193)**2
    m0 = 10000
    mdot = 12.9

    # 操纵变量：推力径向分量、切向分量
    u1 = m.MV(value=0)
    u1.STATUS = 1 
    u1.DCOST = 0

    u2 = m.MV(value=1)
    u2.STATUS = 1 
    u2.DCOST = 0

    # 状态变量
    t = m.Var(value=0, fixed_initial=True)
    r = m.Var(value=1, fixed_initial=True)
    u = m.Var(value=0, fixed_initial=True) # 径向速度
    v = m.Var(value=np.sqrt(mu), fixed_initial=True) # 切向速度

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

    # 轨道动力学方程
    m.Equation(t.dt() == 1)
    m.Equation(r.dt() == u)
    m.Equation(u.dt() == v**2/r - mu/r**2 + thr*u1/(m0-mdot*t))
    m.Equation(v.dt() == -u*v/r + thr*u2/(m0-mdot*t))

    # 约束
    m.fix(u, pos=nt-1, val=0.0) # 终端径向速度为0
    m.Equation(u1**2 + u2**2 == 1) # 推力幅值固定
    m.Equation(final * (v**2 - mu/r) == 0) # 终端圆轨道速度约束

    # 目标：最大化终端轨道半径
    m.Obj(-r*final)

    # 求解配置
    m.options.IMODE = 6
    m.options.NODES = 4
    m.options.MV_TYPE = 1
    m.options.SOLVER = 3
    m.solve(disp=False)

    # 结果可视化（可选）
    plt.figure(figsize=(10,6))
    plt.subplot(221)
    plt.plot(m.time, r.value)
    plt.ylabel('Orbit Radius r')
    plt.subplot(222)
    plt.plot(m.time, u.value, label='Radial velocity')
    plt.plot(m.time, v.value, label='Tangential velocity')
    plt.legend()
    plt.subplot(223)
    plt.plot(m.time, u1.value, label='Radial thrust')
    plt.plot(m.time, u2.value, label='Tangential thrust')
    plt.legend()
    plt.tight_layout()
    plt.show()

求解得到的终端半径约为1.52，与该问题的标准参考结果一致。

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