如何用Gekko求解NLP形式轨迹优化问题？遇方程长度超限错误

1. 批量优化初始与目标条件约束 #

无需逐个元素循环设置约束，利用Gekko的数组操作和m.Equations()可批量添加约束，简化代码：

# 批量设置初始状态约束
x0_vars = Z[self.x_index[:,0]]
m.Equations(x0_vars == xic)

# 批量设置目标状态约束
xf_vars = Z[self.x_index[:,N-1]]
m.Equations(xf_vars == xgoal)

Gekko支持数组级别的等式约束，会自动拆解为单个元素的约束，避免冗余循环。

2. 解决「Max Equation Length」错误 #

该错误源于将所有时间步的目标项合并为一个超大表达式，超出Gekko单方程长度限制。解决核心是逐个时间步拆分目标项，而非合并成单一表达式：

# 逐个时间步添加目标项，避免超长方程
for i in range(N-1):
    x = Z[self.x_index[:,i]]
    u = Z[self.u_index[:,i]]
    # 计算当前时间步的状态代价（用m.sum替代numpy矩阵运算）
    cost_x = 0.5 * m.sum([Q[j,j]*(x[j]-xgoal[j])**2 for j in range(nx)])
    # 计算当前时间步的控制代价
    cost_u = 0.5 * m.sum([R[j,j]*u[j]**2 for j in range(nu)])
    # 逐个累加目标函数
    m.Minimize(cost_x + cost_u)

# 添加终端代价（如果需要）
xf = Z[self.x_index[:,N-1]]
cost_xf = 0.5 * m.sum([Qf[j,j]*(xf[j]-xgoal[j])**2 for j in range(nx)])
m.Minimize(cost_xf)

关键说明： #

• 用m.sum()替代numpy的矩阵乘法（@），Gekko变量不支持直接通过numpy运算生成超长表达式，m.sum()会将求和拆解为多个小项，规避长度限制。
• 多次调用m.Minimize()时，Gekko会自动累加所有代价项到目标函数，无需手动合并。

完整修改后代码示例 #

import numpy as np
from gekko import GEKKO

nx = 13
nu = 4
tf = 10.0  # 示例时间总长
h = 0.1    # 示例步长
N = int(tf/h)+1
nz = nx*N + nu*(N-1)  # 状态+控制总变量数

Q = np.diag(np.ones(nx))
R = 0.1*np.diag(np.ones(nu))
Qf = 10*np.diag(np.ones(nx))    

# 示例初始状态与目标状态
xic = np.zeros(nx)
xgoal = np.ones(nx)

# 示例索引数组生成（替换为你的实际索引）
x_index = np.arange(nx*N).reshape(nx, N)
u_index = np.arange(nx*N, nx*N+nu*(N-1)).reshape(nu, N-1)

# 初始化模型
m = GEKKO(remote=False)
m.options.IMODE = 3  # 稳态优化模式
m.options.SOLVER = 3              
Z = m.Array(m.Var, nz, value=0, lb=-1e5, ub=1e5)

# 设置控制变量下界
for i in range(N-1):
    for j in u_index[:,i]:
        Z[j].lower = 0

# 批量约束
x0_vars = Z[x_index[:,0]]
m.Equations(x0_vars == xic)
xf_vars = Z[x_index[:,N-1]]
m.Equations(xf_vars == xgoal)

# 目标函数构建
for i in range(N-1):
    x = Z[x_index[:,i]]
    u = Z[u_index[:,i]]
    cost_x = 0.5 * m.sum([Q[j,j]*(x[j]-xgoal[j])**2 for j in range(nx)])
    cost_u = 0.5 * m.sum([R[j,j]*u[j]**2 for j in range(nu)])
    m.Minimize(cost_x + cost_u)

# 终端代价
cost_xf = 0.5 * m.sum([Qf[j,j]*(xf_vars[j]-xgoal[j])**2 for j in range(nx)])
m.Minimize(cost_xf)

# 求解
m.solve(disp=True)

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