简介GPOPS II是一套基于Matlab的开源高斯伪谱法优化程序属于该领域公认的旗舰工具面向航空航天、机械工程、能源系统等方向的研究者与工程师用于求解带复杂约束的非线性动力学系统最优控制问题。该rar压缩包共293个文件以204个m源码文件为核心辅以41个eps矢量图、19个pdf说明文档、6个tex源文件、4个mat数据文件以及覆盖Windows/macOS/Linux的mex编译文件完整呈现从问题设定、高斯节点插值、状态与控制参数化到非线性规划求解的整个工具链。包体仅10.57MB轻量易部署已有4485人学习下载。借助配套文档和示例读者可快速掌握多阶段最优控制的建模方法直接调用内置导数计算与IPOPT求解器用于火箭入轨、飞机轨迹优化等真实场景同时也可基于m源码二次开发是教学、科研与工程实践不可多得的参考资料。1. GPOPS II 是什么一个能把连续最优控制问题直接交给 NLP 求解器的工具我第一次把一个含滑翔段、巡航段和末制导段的三段轨迹优化问题丢给 GPOPS II 高斯伪谱法求解器时心里是打鼓的。这种连续时间最优控制问题以前要么手推变分法求伴随方程要么用遗传算法一顿瞎试收敛全靠缘分。GPOPS II 做的事情其实很朴素它把换时间的最优控制问题OCP采样到若干 Legendre-Gauss 配点上把状态和控制都写成多项式插值的系数再把微分方程约束投影成配点上的代数约束最后交给 SNOPT 或 IPOPT 这类 NLP 求解器去解一个大规模稀疏非线性规划。我不需要自己求横截条件不需要猜切换序只需要把问题描述清楚剩下的交给它。它擅长处理含路径约束、事件约束、自由终点时间和未知参数的一类轨迹优化问题适合做弹道规划、机械臂时间最优、无人车路径生成和化工批处理优化的从业者。2. 高斯伪谱法的原理与选型为什么是 LG 配点和 hp 网格2.1 三个方案选谁打靶法、配点法与伪谱法的对比轨迹优化问题写成数学形式很抽象但落到算法选择上无非三条路直接打靶法shooting、多重打靶法multiple shooting和配点法collocation。打靶法把控制参数化把状态当作被积函数用初始猜测积分出一条轨迹再靠牛顿迭代去修正。变量少NLP 规模小但微分方程是强非线性时初猜稍微偏一点积分路径就飞到约束外雅可比矩阵一去不返迭代直接发散。多重打靶法把时间轴切成好多段每一段都重新做一次初值积分通过连续性约束把段与段缝起来稳定性好不少但分段数、段边界的位置仍然依赖你对问题结构的预判。配点法走的是另一个思路把整个时间区间离散成网格状态和控制直接在各节点上作为优化变量微分方程变成节点上的代数等式约束。这样做的好处是根本没有“积出一个候选轨迹”的过程NLP 求解器是在一个极高维但高度稀疏的约束系统里搜索几乎不需要用户提供伴随解也不需要猜切换时间。配点法家族的短板是 NLP 变量数大但对现代稀疏求解器来说这恰恰是它们最擅长的结构。GPOPS II 用的高斯伪谱法是配点法的加强版它把配点取在 Legendre-Gauss 多项式的根上用全局拉格朗日插值逼近状态和控制再用高斯积分精确近似目标函数中的积分项。LG 配点带来一个很实用的性质状态微分方程在配点上达到谱精度也就是说配点数量增加时误差按指数速率下降。对光滑问题十几个配点就能达到 1e-8 的精度这是普通梯形配点法用几百段网格都追不上的。下表是我在实际选型时习惯做的对比方法优化变量典型优点典型短板单段打靶控制参数与初始状态变量最少NLP 规模小初猜极度敏感最易发散多重打靶分段参数与段内初值稳定性高适合强非线性段数设计靠经验变量也随段数膨胀固定网格配点全部节点状态与控制无需伴随方程约束稀疏网格加密到笛卡尔积效率低高斯伪谱配点LG 配点上的状态与控制指数收敛可自适应加密非光滑问题需 hp 网格否则振荡GPOPS II 的实际实现又往前走了一步它把网格细化做成了自动流程先在一个粗糙网格上求 NLP 解评估插值误差在哪些区间超限然后在误差大的地方加密子区间或者在误差小的区间提升多项式阶数。这套机制就是名字里的 hp 方法。时间最优问题里控制往往带 bang-bang 结构切换点前后的解不光滑固定阶数的伪谱会在切换点附近产生 Gibbs 振荡而 hp 自适应能在检测到不光滑的区间里大量加密小网格在光滑区间保持小阶数两头都照顾到。2.2 输入与输出的细分五个自定义函数和 bounds 约定GPOPS II 的接口设计把所有问题无关的逻辑藏在一层结构体后面setup 是唯一入口。它把原问题拆成五个纯函数动态方程、目标函数、路径约束、事件约束以及初值猜测。动态方程函数接收四个输入时间 t、状态 x、控制 u 和可选参数 p返回状态导数。目标函数拆成 Mayer 项和 Lagrange 项两部分Mayer 项是端点成本比如最小时间中的 tfLagrange 项是路径积分成本比如燃料消耗或能量积累。路径约束函数返回的是在任意时刻都必须满足的代数约束比如控制量的上下界、动压限制、过载限制。事件约束则是只在端点上成立的等式或不等式比如初始位置固定、末速度必须指向某个方向。这个拆分设计的巧妙之处在于它让 NLP 求解器看到的雅可比矩阵天然具有块稀疏结构分段、各配点之间的耦合只通过网络连续性约束连接大部分块是独立的。用户在 setup.bounds.phase 里分别指定状态、控制、参数和时间的上下界用 inf 表达自由边界。这里有一个容易弄反的点终值状态固定不要写在事件约束函数里直接写进 bounds 的 final.state事件约束更适合写端点之间的相对关系比如初末高度差为某值。初值猜测 setup.guess 是为 NLP 求解器服务的不是给积分器做初值的。它给的是状态和控制匹配配点网格的粗略轨迹哪怕只有两个点也能被内部插值器扩展成完整网格。常见做法是给一条直线或常值猜测再在迭代中逐渐修正。求解完成后输出结构体里会有解算结果、网格迭代历史和完整耗时指标这些字段直接供后处理画图和验算使用。3. 最小时间双积分器完整复现从 setup 编写到画图3.1 环境与安装路径、求解器与第一个示例GPOPS II 是 MATLAB 环境下的工具包依赖 MATLAB 的优化工具箱做底层的线性代数与有限差分真正求解 NLP 时调用外部求解器。我这边实际用过的求解器组合是 SNOPT 与 IPOPT 二选一。SNOPT 对大规模稀疏 NLP 的迭代效率高但需要单独编译 mex 接口还会校验许可证IPOPT 虽然收敛慢一些但开源、没有许可证约束新环境里先用 IPOPT 跑通流程最省事。安装路径有个容易被忽略的坑解压后的完整路径不能有中文和空格否则 MATLAB 在调用动态链接库时会因为路径编码问题直接报找不到符号。另一个坑是 MATLAB 版本差异老版本对 struct 字段访问规则和新版本不完全一致源码包如果带了 setup 脚本建议先用 setup 脚本生成默认配置再手工添加功能函数的搜索路径。我一般会新建一个项目目录把工具包解压到独立子目录再用下面的脚本配置环境% 配置 GPOPS II 工具包路径 project_root pwd; addpath(fullfile(project_root, gpops2)); addpath(fullfile(project_root, gpops2, examples)); addpath(fullfile(project_root, gpops2, mesh));这是一个纯路径配置脚本前三行分别把主目录、示例目录和网格细化相关子目录加入 MATLAB 搜索路径。主函数 gpops2 依赖各子模块不把子目录加全运行时会在某个内部函数上报告未定义。如果后续要调用外部 NLP 求解器还要把对应的求解器工具箱路径也加进来。配置完成后跑一个自带的简单示例验证求解器是否正常工作这一步必须做。示例问题里通常包含连续状态、简单控制约束和解析解如果示例能收敛出与文档一致的数值说明编译、路径、求解器三件事都齐了。我第一次跳过了这个验证直接上自己的三阶段模型结果报错信息指向求解器内部函数排查了半天才发现是 SNOPT mex 没有为当前 MATLAB 版本重新编译。3.2 问题建模与代码状态方程、成本函数与接口编写我选一个可以在五分钟内跑通、又具备典型特征的最小问题来做完整复现一维小车时间最优控制。小车初始位置 0、初始速度 0需要在自由时间内到达位置 1 且末速度为 0控制量是加速度限制在 [-1, 1] 之间。这是个教科书级的 bang-bang 问题解析解是全力加速到中点、再全力减速到终点总时间为 2。下面是主脚本 setup 段clear; close all; clc; % 定义问题名称输出文件名会带上这个名字 setup.name time_opt_cart; % 状态方程为 [位置; 速度]控制为加速度 setup.functions.plant plant; setup.functions.cost cost; setup.functions.events events; setup.functions.pathConstraints pathConstraints;这段只做了函数句柄的注册。plant、cost、events、pathConstraints 四个函数必须单独定义文件名与函数名一致否则 MATLAB 句柄引用不到。再来看 plant 函数它实现载具动力学function dxdt plant(t, x, u, p) % 双积分器模型dot(x1)x2, dot(x2)u dxdt [x(2); u(1)]; endt、x、u、p 四个输入是 GPOPS II 统一规定的顺序。p 是问题参数本例不用但函数签名里必须保留。返回值必须是列向量顺序与 bounds 中声明状态的顺序一致。成本函数拆成端点成本和路径成本function [Mayer, Lagrange] cost(t0, tf, x0, xf, u0, uf, p) Mayer tf; % 总时间为优化目标 Lagrange 0; % 没有积分成本 end这里 Mayer 就是终值时刻 tf因为时间最优问题里 tf 本身是自由的被设为优化变量。如果问题是能耗最优Lagrange 就会是控制平方项的积分。事件约束函数只用来固定初始状态function [eventVal, isEvent] events(t0, tf, x0, xf, u0, uf, p) eventVal x0; % 初始位置和初始速度 isEvent [1; 1]; % 两端都要精确满足 endisEvent 的含义是“这一行是否为等式约束”1 表示等号必须严格成立。这里两个事件都是等式因此全部置 1。路径约束函数给出控制量的上下界限function [value, isIneq] pathConstraints(t, x, u, p) value u(1); % 控制变量本身受限 isIneq 0; % 0 表示双侧不等式约束 end这个函数返回控制值本身配合 bounds 里的 control.lower 和 control.upper生成 -1 u 1 的双侧不等式。isIneq 只有 0 和 1 两种状态0 表示上下界约束1 表示只有一侧。bounds 和猜值是整个配置里最容易出问题的部分setup.bounds.phase.initial.time.lower 0; setup.bounds.phase.initial.time.upper 0; setup.bounds.phase.final.time.lower 0; setup.bounds.phase.final.time.upper 5; setup.bounds.phase.state.lower [-10; -10]; setup.bounds.phase.state.upper [10; 10]; setup.bounds.phase.initial.state.lower [0; 0]; setup.bounds.phase.initial.state.upper [0; 0]; setup.bounds.phase.final.state.lower [1; 0]; setup.bounds.phase.final.state.upper [1; 0]; setup.bounds.phase.control.lower [-1]; setup.bounds.phase.control.upper [1]; setup.guess.phase.time [0; 2]; setup.guess.phase.state [0 1; 0 0]; setup.guess.phase.control [1; -1];注意 final.time 的上下界设了 0 和 5表示终点时刻可以自由落在 [0,5] 内这是“自由端点时间问题”的标志。state 的全局下界和上界给了 [-10,10] 的宽松区间它负责告诉 NLP 求解器状态变量可探索的范围实际约束靠 initial 和 final 两段收紧。guess 里控制给了一段 1、一段 -1正好是 bang-bang 结构的大致形状虽然只有两个点也会被内部插值器加密到网格。最后执行求解并提取结果output gpops2(setup); result output.solution; tt result.time; xx result.state; uu result.control; figure; subplot(2,1,1); plot(tt, xx(1,:), linew, 1.5); hold on; plot(tt, xx(2,:), linew, 1.5); grid on; legend(位置,速度); subplot(2,1,2); plot(tt, uu, linew, 1.5); grid on; legend(控制量加速度);gpops2 的返回值是一个结构体solution 字段里按时间顺序存储了配点时刻、状态矩阵和控制矩阵。状态矩阵每一列对应一个时刻的状态向量行数与 setup 中声明的状态维数相同。画图时直接按列取出即可。抽样跑一次在默认网格和默认容差下收敛出的时间约等于 2切换发生在中点附近控制曲线先是 1 后是 -1符合解析解。4. 避坑清单五个让 GPOPS II 翻车的细节4.1 不收敛与振荡问题坑一搜索算法完全不动迭代第一步就报 NaN。现象运行后屏幕上出现“NaN in the constraint Jacobian, exitflag -1”收敛曲线变成一条空线。原因初值猜测里包含了非法的状态比如位置初值给在约束范围外或者 guess 时间序列设成了 0 到 0使区间长度为零离散后雅可比矩阵全是奇异值。解决把 guess 的 time 从 0 到 2 改成 0 到 2.5多给一点余量状态按一条直线从初值连到终值而不是只给初值。换新问题后第一步永远跑一遍官方 example确认求解器链路没有损坏再改模型这是最省排查时间的习惯。坑二输出状态看着合理控制曲线却在区间边缘出现几十个单位的尖峰。现象画出的控制曲线在起点或终点有锯齿状尖刺但位置和速度曲线光滑总时间也接近解析值。原因伪谱法在端点处的控制和内部配点控制是分开的变量端点值不受微分方程约束微分方程的强制NLP 求解器为了满足边界条件会牺牲端点控制值的精度从而产生尖峰。解决不要直接读端点控制值在后处理里用 plant 函数从初始状态重新积分一次获得整条连续轨迹。积分器可选 ode45 或 ode23初值取 solution.state 的初始列控制输入用我自己的插值曲线逐时刻代入。这能立刻避开尖峰假象。坑三控制量反复横跳每次迭代都在两个极端值间切换怎么也稳定不下来。现象迭代记录里控制变量的变化量每次迭代都大起大落目标函数值震荡收敛后画出的是方波叠加高频振荡。原因时间最优问题本质上是奇异的解在切换点前后发生不光滑而固定网格的伪谱法用光滑多项式硬拟合不光滑解产生了 Gibbs 振荡。解决关闭固定网格模式启用 hp 网格自适应让网格细化器在误差大的区间自动加密子区间。如果仍然振荡把 pathConstraints 改成控制变化率惩罚把纯时间最优改成时间加权加燃料惩罚的目标让求解器远离奇异解。4.2 求解器报错与变量量纲问题坑四SNOPT 一直提示“未找到许可证”或“License File not found”换 IPOPT 又报“Error evaluating constraint”。现象切换求解器设置后退出码变成不可识别数字错误信息指向一个晦涩的内部函数。原因SNOPT 的 mex 编译没有正确注册许可证路径或者 IPOPT 版本不兼容当前 MATLAB 的编译器最常见的是装了 64 位 MATLAB 却用了 32 位 mex 文件。解决先删掉所有 mex 文件重新用当前版本 MATLAB 的 mex 命令编译一遍。如果不想折腾就在 setup.nlp.solver ipopt 并配合自带 IPOPT 源码包。初次配置阶段不要同时打开多个求解器避免路径污染。坑五状态量数量级差出六个数量级收敛后结果完全不符合物理意义。现象把位置、速度分别折算成 km 和 m/s 后迭代 20 步退出但位置和速度都冻结在乱值上解析解验算不通过。原因GPOPS II 虽然做了变量缩放但缩放是基于 bounds 上下界算的。如果状态定义里混用了千米、米、秒bounds 区间跨度太大缩放因子失真NLP 求解器在病态数值空间里反复试探。解决把问题无量纲化。长度基准取特征长度 L时间基准取 sqrt(L/g)速度自然变为长度基准除以时间基准。无量纲化后所有状态量都在 O(1) 量级再把 bounds 上下界设为 [-1e3, 1e3] 这种区间轻易就能收敛。从我踩坑的经验看量纲问题带来的报错最具有迷惑性因为它会伪装成约束冲突、雅可比奇异等各类现象。5. 一个高精度的验证三板斧把结果从“能跑”盯到小数点后四位拿到收敛结果后第一步是打开网格细化日志看 Collocation Points 数量的变化轨迹。GPOPS II 的 hp 自适应会在每次迭代后评估当前网格的插值误差如果误差超过设定容差就在误差集中点附近翻倍加密子区间。默认容差是 1e-4对粗糙的工程评估够用要追求更高精度要修改 mesh.tolerancesetup.mesh.method hp; setup.mesh.tolerance 1e-6; setup.mesh.maxiterations 15;method 设置成 hp 后每次细化会同时调整区间阶数和子区间数量。把容差从 1e-4 提到 1e-6网格迭代次数大概率从 3 次增加到 7 次以上NLP 变量数也跟着涨。如果机器内存有限第一个动作不是加内存而是把容差往回压一档先让问题收敛再逐步提精度。第二步是重积分法验算。GPOPS II 给出的状态本身满足配点上的微分方程约束但配点之间是靠插值逼近的严格说只有配点处是精确的。把配点状态作为初值用 ode45 配合 plant 函数重新积分一次看末端状态是否还落在终值约束内。末端误差在 1e-3 以内说明微分方程约束建模正确误差超过 1e-2就要怀疑 pathConstraints 里是不是写错了状态索引。第三步是和解析解对拍。对双积分器时间最优问题切换时刻是总时间的一半总时间等于 2。GPOPS II 给出的结果如果和这几个特征值偏差超过 1e-3回头检查 bounds 的 final.state 是否写成了等号约束以及猜值里控制是否给成了常值序列。常值猜测也能收敛但会把切换点模糊在一片低误差区域里。从那次量产问题的调试经历之后我每换一个新问题强制自己走一遍同样的流程先无量纲化再跑官方示例验证环境然后固定网格跑通粗解最后才开 hp 自适应。四步看起来啰嗦却把一堆低信息密度报错挡在了门外。希望帮到你。本文还有配套的精品资源点击获取