Gekko求解航天器交会轨迹优化问题遇无解错误求助

安全阿岳

2026-8-2

航天器交会轨迹优化问题求解失败分析与改造方案

我是Arthur，刚接触Gekko，正在求解航天器交会（追踪器→目标器）的轨迹优化问题：

需满足最终相对距离、速度约束及航行过程中的距离约束
尝试对最终路径采用软约束，以推力的L1范数作为目标函数
计算持续至最大迭代次数（2000次）后终止，未找到解。想知道是否可通过离散化搜索空间在精度与计算时间间权衡，使问题可解？同时未找到infeasibilities.txt文件，恳请协助分析无解原因并提供问题可行化改造方案。

问题代码

import matplotlib.pyplot as plt
import numpy as np

from gekko import GEKKO

# create GEKKO model
m = GEKKO()
print(m.path)
nt = 501
m.time = np.linspace(0,500,nt)

# Variables

# Initial Position and Velocity - Chaser
r1_c = -1182.339411   # [km] 
r2_c = 6816.939420    # [km] 
r3_c = 904.891745     # [km] 
v1_c = 0.1175776      # [km/s] 
v2_c = -0.963776      # [km/s] 
v3_c = 7.494102       # [km/s] 

# Initial Position and Velocity - Target
r1_t = -1182.959348    # [km] 
r2_t = 6817.396210     # [km] 
r3_t = 904.495486      # [km] 
v1_t = 0.1175776       # [km/s] 
v2_t = -0.963776       # [km/s] 
v3_t = 7.494102        # [km/s] 

# initial values
x =  m.Var(value = r1_c)
y =  m.Var(value = r2_c)
z =  m.Var(value = r3_c)
vx =  m.Var(value = v1_c)
vy =  m.Var(value = v2_c)
vz =  m.Var(value = v3_c)


# initial values
x2 =  m.Var(value = r1_t)
y2 =  m.Var(value = r2_t)
z2 =  m.Var(value = r3_t)
vx2 =  m.Var(value = v1_t)
vy2 =  m.Var(value = v2_t)
vz2 =  m.Var(value = v3_t)

# Cost Initial - Integrated thrust (fuel usage)
J =  m.Var(value = 0)

# Manipulated Variable
U = m.MV(value=0, lb=-0.001, ub=0.001)

U.STATUS = 1

# Mask time
p = np.zeros(nt) # mask final time point
p[-1] = 500
final = m.Param(value=p)

# Parameters
RT =  m.Const(6378.139);            # [km] 
mu_t = m.Const(3.9860064*10**(5))   # [km^3/s^2] 


# Equations

# Define intermediate quantities - Chaser
Rc = m.Intermediate((x**2 + y**2 + z**2)**0.5)
v = m.Intermediate((vx**2 + vy**2 + vz**2)**0.5)

ax = m.Intermediate(x * -mu_t / Rc**3 + U *vx / v )
ay = m.Intermediate(y * -mu_t / Rc**3 + U *vy / v )
az = m.Intermediate(z * -mu_t / Rc**3 + U *vz / v )


# Define intermediate quantities - Target
Rt = m.Intermediate((x2**2 + y2**2 + z2**2)**0.5)
v2 = m.Intermediate((vx2**2 + vy2**2 + vz2**2)**0.5)

ax2 = m.Intermediate(x2 * -mu_t / Rt**3 )
ay2 = m.Intermediate(y2 * -mu_t / Rt**3 )
az2 = m.Intermediate(z2 * -mu_t / Rt**3 )


# Governing equations
m.Equations((vx.dt() == ax, vy.dt() == ay, vz.dt() == az))
m.Equations((x.dt() == vx, y.dt() == vy, z.dt() == vz))

m.Equations((vx2.dt() == ax2, vy2.dt() == ay2, vz2.dt() == az2))
m.Equations((x2.dt() ==  vx2, y2.dt() ==  vy2, z2.dt() == vz2))

# Equation relating thrust to fuel usage
m.Equation(J.dt() == m.abs2(U))


# Path Constraints

# specify endpoint conditions

m.fix(x, pos=len(m.time)-1, val=x2)
m.fix(y, pos=len(m.time)-1, val=y2)
m.fix(z, pos=len(m.time)-1, val=z2)
m.fix(vx, pos=len(m.time)-1, val=vx2)
m.fix(vy, pos=len(m.time)-1, val=vy2)
m.fix(vz, pos=len(m.time)-1, val=vz2)

# Constraints 

m.Equations((m.abs2(Rc - RT) > 200, m.abs2(Rc - RT) < 1000, m.abs2(Rc*final - Rt*final) == 0, m.abs2(v*final - v2*final) == 0))


# Objective, minimizing fuel usage
m.Obj(J * final)

# Set solver mode - Optimal Control
m.options.IMODE = 6

# Increase maximum number of allowed iterations
m.options.MAX_ITER = 2000

# Set number of nodes per time segment
m.options.NODES = 3

# Run solver and display progress
m.solve(disp=True)

求解器输出

---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :    3663.92820000000      sec
 Objective      :    185346862893.104     
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found
---------------------------------------------------------------------------

无解原因分析

约束冲突与硬约束过强：用m.fix()强制追踪器最终状态与目标器完全一致，同时叠加重复的最终状态约束，易导致约束不可行；m.abs2(Rc - RT) > 200这类非光滑约束会大幅增加求解难度，且需先验证初始状态是否满足路径约束。
控制量与目标函数设计错误：仅设置单一标量推力U，限制了控制自由度，无法灵活调整轨迹；用m.abs2(U)计算的是L2范数，并非需求的L1范数。
数值求解负担过重：nt=501的时间节点加上NODES=3，导致变量规模过大，延长求解时间且易陷入局部无解；infeasibilities.txt生成在m.path输出的临时路径中，本地求解时需到该路径查看。
软约束未正确实现：代码声称使用软约束，但实际用了m.fix()硬约束，未引入松弛变量或惩罚项处理约束违反情况。

可行化改造方案

1. 修正约束设置

替换硬约束为软约束：移除m.fix()，引入松弛变量和惩罚项，允许最终状态存在小偏差：

# 定义松弛变量
slack_pos = m.Var(value=0, lb=0)
slack_vel = m.Var(value=0, lb=0)
# 最终状态约束
m.Equation(m.abs2(x[-1]-x2[-1]) <= slack_pos)
m.Equation(m.abs2(y[-1]-y2[-1]) <= slack_pos)
m.Equation(m.abs2(z[-1]-z2[-1]) <= slack_pos)
m.Equation(m.abs2(vx[-1]-vx2[-1]) <= slack_vel)
m.Equation(m.abs2(vy[-1]-vy2[-1]) <= slack_vel)
m.Equation(m.abs2(vz[-1]-vz2[-1]) <= slack_vel)
# 惩罚松弛变量，权重可根据需求调整
m.Obj(1e6*slack_pos + 1e6*slack_vel)

修正路径约束为光滑形式：将m.abs2(Rc - RT) > 200替换为Rc - RT >= 200或RT - Rc >= 200（根据安全距离方向），避免非光滑约束；若初始状态不满足路径约束，需调整初始条件或放松约束范围。

2. 优化控制量与目标函数

改为三维推力控制：设置三个独立的推力MV，增加控制自由度：

Ux = m.MV(value=0, lb=-0.001, ub=0.001); Ux.STATUS=1
Uy = m.MV(value=0, lb=-0.001, ub=0.001); Uy.STATUS=1
Uz = m.MV(value=0, lb=-0.001, ub=0.001); Uz.STATUS=1
# 修正加速度方程
ax = m.Intermediate(x * -mu_t / Rc**3 + Ux )
ay = m.Intermediate(y * -mu_t / Rc**3 + Uy )
az = m.Intermediate(z * -mu_t / Rc**3 + Uz )

实现真正的L1范数目标：使用m.abs()计算推力的L1范数积分：
```
m.Equation(J.dt() == m.abs(Ux) + m.abs(Uy) + m.abs(Uz))
```

3. 降低求解难度

减少时间节点数量：先从nt=101、NODES=2开始测试，验证可行性后再逐步增加节点数量。
调整求解器参数：切换到APOPT求解器（m.options.SOLVER=3），它对非线性和非光滑问题的处理更友好；适当降低精度要求（m.options.RTOL=1e-4、m.options.ATOL=1e-4），帮助求解器找到可行解。
使用相对坐标系：将追踪器状态改为相对目标器的位置和速度，简化动力学方程，减少变量数量。例如采用Clohessy-Wiltshire方程（适用于近圆轨道相对运动），目标器的轨道可预先计算为参数，无需作为变量求解。

4. 分步验证模型

先求解无约束的交会问题，验证动力学模型正确性；
逐步加入路径约束，观察求解器是否能找到可行解；
调整初始猜测值，例如让追踪器初始状态完全跟随目标器，再逐步调整到实际初始值。

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

Gekko求解航天器交会轨迹优化问题遇无解错误求助

问题代码 #

求解器输出 #

无解原因分析 #

可行化改造方案 #

1. 修正约束设置 #

2. 优化控制量与目标函数 #

3. 降低求解难度 #

4. 分步验证模型 #