简介GPOPS II 是求解非线性动态系统最优控制问题的高斯伪谱法优化程序可用于飞行器轨迹优化、机械控制和能源系统等工程场景适合具备一定数值优化基础的研究人员与工程师。资源包共 293 个文件压缩后约 10.57MB其中以 204 个 MATLAB 源文件为主体包含实例脚本和辅助函数另配有 19 个 PDF 文档、LaTeX/BibTeX 排版源码、41 个 EPS 矢量图以及 MAT 数据文件还提供 Windows/Linux/Mac 下的 MEX 编译版本便于跨平台直接调用目录结构较清晰。内容涉及火箭级间分离、飞行攻角与航迹角等典型示例完整展现了从问题建模、高斯节点插值、状态与控制参数化到非线性规划求解的流程。已有 4478 人学习使用其源码、文档与图表齐全适合作为 GPOPS II 的入门学习参考、论文复现素材和工程优化分析的实用工具基础。 做最优控制这么多年求解器换了一茬又一茬但每次有人问我“现在用什么工具解轨迹优化”我大概率会先反问一句你接触过GPOPS II吗倒不是说我非要推荐它而是这个工具在“高斯伪谱法”这个技术路线上确实做到了某种程度的极致——你给它一个非线性最优控制问题它能在配点策略、网格细化、求解器衔接这些环节替你省掉一大半的功夫。网上有人叫它“史上最牛的高斯伪普法优化程序”这个称号带点夸张但就实际体验来说GPOPS II 在同级别开源/半开源工具里确实有很高的完成度。这篇文章我不想写那种“安装点下一步、运行看结果”的说明书式教程而是想从原理和实操两个维度聊聊高斯伪谱法到底解决了什么问题GPOPS II 的关键设计是什么以及你在用它建模、调试、踩坑时最该注意哪些细节。全文内容偏实战适合正在做飞行器轨迹优化、机器人运动规划、过程控制优化以及任何跟非线性最优控制打交道的人参考。1. 高斯伪谱法到底解决了什么难题要理解 GPOPS II 的厉害之处得先回到它背后的数学方法。最优控制问题本质上是一个无穷维优化问题你要求解的不是一组静态参数而是一条随时间变化的状态轨迹和控制轨迹同时要满足动力学方程、边界条件和路径约束。1.1 从“无限维”到“有限维”的转化经典变分法需要手推哈密顿函数、伴随方程和横截条件问题稍微复杂一点比如路径约束带不等式、状态变量有上下界推导量就会爆炸。数值方法阵营里直接打靶法和多重打靶法把状态和控制按时间网格离散然后用常微分方程积分器向前传播思路很直观但对初值极其敏感而且遇到刚性问题时收敛性很差。伪谱法走了另一条路把连续时间的状态和控制函数用全局插值多项式逼近然后在离散点上强制满足动力学方程。这样就把微分方程约束变成了代数约束最优控制问题变成了非线性规划问题NLP。核心区别在于打靶法是“猜初始值→积分→修正”伪谱法是“把整个时间域离散化后一次性优化”前者是迭代推进逻辑后者是全局配平逻辑所以伪谱法在处理复杂约束时天然更稳定。1.2 为什么偏偏是“高斯”配点伪谱法里有几个分支切比雪夫伪谱法用 Chebyshev-Gauss-Lobatto 点Legendre 伪谱法用 Legendre-Gauss-Lobatto 点而 GPOPS II 用的是 Legendre-GaussLG配点。这里有个看似反直觉的点LG 点不包含时间域的两个端点 -1 和 1而 LGL 点包含端点。为什么不直接用包含端点的配点原因藏在协态映射定理costate mapping theorem里使用 LG 配点时离散化后的一阶最优性条件KKT 条件与连续时间的变分最优性条件之间有一种精确的映射关系协态变量在配点处的值可以直接从 NLP 的拉格朗日乘子还原出来不需要额外的后处理步骤。而 LGL 和 CGL 配点由于端点被强制满足约束这种映射关系会出现退化导致协态估计精度下降。换句话说GPOPS II 选用 LG 配点不只是一种数值偏好而是为了保证“离散解的最优性条件”和“连续问题的最优性条件”在数学上严格一致。这一点在需要事后验证解的最优性、或者需要提取协态变量做进一步分析时非常关键。1.3 GPOPS II 与 GPOPS 的区别GPOPS II 不是简单地在第一代基础上修修补补它在算法层面做了三个重要升级第一采用 hp 自适应网格细化策略。这一点稍后会展开讲是它拉开与一代差距的核心能力。第二对稀疏有限差分雅可比矩阵做了大幅优化配合 SNOPT 和 IPOPT 这类稀疏 NLP 求解器时效率提升明显。第三内置了更干净的接口层用户只需要提供连续时间的最优控制问题描述不需要手动处理配点和离散化。所以回到“史上最牛”这个评价它牛的不是某一项指标而是在“求解能力、易用性、数学严谨性”这三者之间取得了很好的平衡。2. GPOPS II 的整体设计与方案选型逻辑我在实际使用 GPOPS II 时最直观的感受是它把用户的思考重心从“怎么离散化、怎么求梯度”拉回到了“怎么描述问题本身”。这种设计哲学值得拎出来说说。2.1 四个核心函数定义问题的方式决定了上手指数GPOPS II 的顶层接口设计得很干净。你定义一个最优控制问题需要完成四件事描述连续时间动力学、声明边界条件与路径约束、给出性能指标目标函数、提供初始猜测。对应到代码上就是四个函数或数据结构continuous函数定义微分方程 dx/dt f(x, u, t) 和路径约束 g(x, u, t)endpoint函数定义初始状态、末端状态、事件约束以及目标函数中的 Mayer 项和 Lagrange 项setup结构体整合上述函数并指定时间区间、状态/控制变量的边界、初始猜测求解器调用GPOPS II 内置了与 IPOPT 和 SNOPT 的接口你只需要把 NLP 的 option 传进去这里有个很关键的设计理念GPOPS II 把“连续问题”和“离散化求解”彻底解耦了。你在continuous函数里写的方程就是连续时间域的原始方程不需要关心 LG 配点上有多少个点、离散化后的 NLP 变量如何排序。所有网格相关的信息会作为附加参数传给 continuous 函数你只要按格式读取 t、x、u 就能算。2.2 为什么用 hp 自适应网格而非固定网格这是 GPOPS II 设计中最精髓的部分值得多说几句。固定网格伪谱法的思路是预先选定配点数量一次求解完事。这个方法在光滑问题上表现不错但遇到状态轨迹突变比如滑行段和机动段的切换、bang-bang 控制、路径约束激活/失活时单个固定网格很难兼顾“分辨率”和“计算量”。hp 自适应方法借鉴了有限元里 h-refinement 和 p-refinement 的思路在每个网格区间内检查残差如果某个区间的插值误差超过阈值就对这个区间重新细分网格h 细化或者增加该区间内的配点数目p 细化。GPOPS II 的底层有一套误差估计机制能够自动判定“这里需要更多点数”而不是“那里需要更密的网格”。这套机制带来的实际好处是什么我举个例子你就明白了。解一个带有路径约束的最优控制问题约束只在某个时间窗口内激活如果用手动固定网格你得提前猜出激活区间位置然后手动在那个区间加密网格猜错了就得反复试。GPOPS II 会自动发现残差集中区域并完成加密这个能力在复杂工程问题里几乎不可替代。2.3 底层 NLP 求解器的衔接策略GPOPS II 本身不做非线性规划求解它是把离散化后的 NLP 问题交给外部求解器处理。我常用的是 SNOPT 和 IPOPT两者各有侧重求解器优势劣势适用场景SNOPT稀疏 SQP 方法内存占用小处理大规模稀疏问题效率高需要商业许可证对非光滑问题的收敛性依赖初值大规模轨迹优化、传统 aerospace 领域IPOPT开源免费内点法框架对初值敏感度较低鲁棒性好大规模问题时内存占用较高选项参数较多需要调中小规模问题、免费方案的优先选择GPOPS II 在衔接层面做了不少功夫它会把 NLP 问题的雅可比矩阵结构信息以稀疏模式传给求解器并允许用户通过 setup 结构体指定求解器选项。我的建议是能搞到 SNOPT 许可证就优先用 SNOPT搞不到就先用 IPOPT 把方案流程跑通再到换 SNOPT 做精度验证。3. 实操剖析从建模到求解的完整链路理论聊完了下面进入真正的“动笔环节”。我用一个非常经典的问题做载体——月球软着陆燃料最优问题。这不是我随便选的它是一个典型的状态和控制都有约束、目标函数非平凡、而且能明显看出“网格自适应”发挥作用的测试案例。3.1 月球软着陆问题的连续时间数学描述设飞船质量 m高度 h速度 v推力 T 作为控制输入。动力学方程为dh/dt v dv/dt -g T / m dm/dt -T / (Isp * g0)其中 g 是月球重力加速度Isp 是比冲g0 是地球标准重力加速度。路径约束是推力幅值限制 T_min ≤ T ≤ T_max。目标是最小化燃料消耗等价于最大化终端质量 m(tf)或者最小化 ∫T dt。边界条件初始时刻 hH0, v0, mm0终端时刻 h0, v0终端质量自由。这组方程看起来简单但它有非常鲜明的控制切换特性——最优解通常包含最大推力段和零推力段控制剖面有间断点这正好可以用来考验 hp 自适应网格的识别能力。3.2 GPOPS II 的代码实现与关键配置第一步是写动力学方程。GPOPS II 的 continuous 函数签名是dx continuous(x, u, t, auxdata)注意 t 在这个函数里不是向量而是一个标量因为 GPOPS II 将所有配点“摊平”后逐个调用也可以通过 vectorized 模式批量计算。核心代码如下function dx lunar_continuous(x, u, t, auxdata) % 状态向量 x [h; v; m] h x(1,:); v x(2,:); m x(3,:); T u(1,:); g auxdata.g; Isp auxdata.Isp; g0 auxdata.g0; dx zeros(size(x)); dx(1,:) v; dx(2,:) -g T ./ m; dx(3,:) -T ./ (Isp * g0); end这里提前把物理常量放在 auxdata 里而不使用全局变量是为了避免 MATLAB 全局变量在 GPOPS II 自动求雅可比时引发性能问题。第二步是写 endpoint 函数。这个函数负责处理只在初始/终端时刻出现的约束和目标函数。目标函数里有一个 0 时刻到 tf 时刻的时间积分项GPOPS II 中称之为 Lagrange 项需要额外定义积分被积函数function [objective, con] lunar_endpoint(x0, xf, t0, tf, u, t, auxdata) % 目标最小化燃料消耗即最大化终端质量 objective -xf(3); % 因为求解器默认最小化 % 边界条件 con []; con [con; x0(1) - auxdata.h0]; % h0 con [con; x0(2) - 0]; % v0 con [con; x0(3) - auxdata.m0]; % m0 con [con; xf(1) - 0]; % hf 0 con [con; xf(2) - 0]; % vf 0 end function lag lunar_lagrange(x, u, t, auxdata) T u(1,:); lag T; % 积分 ∫T dt end有个细节需要注意目标函数写成-xf(3)是因为 GPOPS II 默认是最小化问题最大化终端质量等价于最小化负终端质量。这个符号问题在最优控制建模里极其容易搞错我见过不少人把目标函数符号写反结果求解器输出一个乱七八糟的解还以为是算法不收敛。第三步是配置 setup 结构体。这里要指定状态和控制变量的边界、初始猜测、网格设置和求解器选项setup.name lunar_landing; setup.functions.continuous lunar_continuous; setup.functions.endpoint lunar_endpoint; setup.functions.lagrange lunar_lagrange; % 如果 objective 有积分项 setup.auxdata.g 1.622; setup.auxdata.Isp 311; setup.auxdata.g0 9.81; setup.auxdata.h0 1000; setup.auxdata.m0 2000; setup.bounds.phase.initialtime.lower 0; setup.bounds.phase.initialtime.upper 0; setup.bounds.phase.finaltime.lower 0; setup.bounds.phase.finaltime.upper 200; setup.bounds.phase.state.lower [0; -10; 100]; setup.bounds.phase.state.upper [1000; 100; 2000]; setup.bounds.phase.control.lower [0]; setup.bounds.phase.control.upper [5000]; setup.bounds.phase.integral.lower [0]; % Lagrange 项的积分范围 setup.bounds.phase.integral.upper [1e7]; % 初始猜测给一条从初始状态到终端状态的线性插值 setup.guess.phase.time [0; 100]; setup.guess.phase.state [[1000; 0; 2000], [0; 0; 1500]]; setup.guess.phase.control [[1000], [1000]]; setup.guess.phase.integral [50000]; setup.mesh.method hp; setup.mesh.tolerance 1e-6; setup.mesh.maxiterations 30; setup.nlp.solver ipopt; setup.nlp.options.tol 1e-8; setup.nlp.options.print_level 3; % 调用求解 output gpops2(setup);3.3 初始猜测的艺术有件事我必须单独拿出来说因为 GPOPS II 的最终求解成败很大程度上不取决于后端 NLP 求解器而取决于你给的初始猜测靠不靠谱。最优控制问题的 NLP 是非凸的伪谱法并不会解决非凸性问题它只负责把问题转成一个“好的”有穷维逼近。如果你的初始猜测离真实解太远SNOPT 的 SQP 迭代很可能在某个局部极小点就停住了。我常用的三种初始猜测构造思路第一物理直觉猜测法。像月球软着陆这种问题至少可以猜一条“先悬停、再下降”的轨迹哪怕是粗糙的线性插值也行。第二先解简化问题。把路径约束去掉、或者把动力学简化成线性模型先用简单模型求出近似解再用这个解作为 GPOPS II 的初始猜测。第三步逐步递增网格。先在粗网格上求一个解再把解作为细网格的初值。GPOPS II 内部 hp 自适应也类似但如果你手动控制网格演化对理解问题性质很有帮助。我在做工程问题时最常用的组合是“物理直觉 简化问题先求一版”效率最高。纯靠 GPOPS II 自己从零开始探索的非线性问题收敛情况很难保证。4. 常见问题与排查技巧实录这部分是从我自己反复折腾 GPOPS II 的过程中总结出来的每条都是踩过坑换来的经验。4.1 状态轨迹出现“锯齿状”振荡但目标函数值已经收敛这是我新手期遇到最多的现象。看起来解出来了画出状态曲线却发现高频抖动。第一次遇到这种情况千万不要急着怀疑求解器先检查一下你给的状态或控制边界是否设得过于宽松。最优控制问题中如果控制变量的上下界范围远大于实际需要的范围NLP 求解器有可能在可行域内部找到一个“目标函数值几乎不变但控制振荡”的解——这在数学上是近似最优解但工程上不可用。解决办法是在保证可行性的前提下尽量收紧边界或者检查 mesh 容差是否需要提高。另外检查你是不是忘了在 continuous 函数里对状态变量做归一化。伪谱法的配点定义在 [-1, 1] 区间GPOPS II 内部会自动做时间缩放但状态变量的量级它不会自动帮你处理。如果 h 是 1e3 量级、v 是 1e2 量级、m 是 2e3 量级而控制 T 是 5e3 量级NLP 的雅可比矩阵数值范围会很差。我的做法是先用无量纲化把状态缩放到 O(1)这个坎过了后面的求解难度直接降一个档。4.2 “Too many function evaluations”或 IPOPT 报 restoration failed这个错误本质上说明 NLP 求解器在可行域外面找不到回退方向了。优先排查三件事第一你的 continuous 函数里是否有不可微的分支逻辑比如用了 if-else 来切换不同的动力学模型GPOPS II 的有限差分雅可比对间断点极其敏感你可能需要把动力学改写成平滑的近似形式。第二路径约束是否在最优点处不可满足比如我要求状态在终端时刻严格等于某个值但控制边界不够大导致动力学上根本无法实现——这种问题无解就是无解求解器卡死是正常的。第三初始猜测是否“物理一致”比如速度边界设了 [-10, 10]初始猜测却给了 50虽然猜测理论上只需要在边界内就行但很多 NLP 求解器对无效初值的鲁棒性并没有你想象的那么好。4.3 协态变量无法还原或还原后与解析解对不上虽然我说 LG 配点有严格的协态映射定理但实际使用中很多人还是翻车。最典型的原因是目标函数里的积分项Lagrange 项没有正确地通过setup.functions.lagrange传给求解器或者在 endpoint 函数里把积分项和 Mayer 项组装错了。GPOPS II 对目标函数的处理有明确分工endpoint 函数返回的objective只是 Mayer 项端点的函数如果你还有积分型的 Lagrange 项必须额外提供 lagrange 函数并且 setup 结构体里要声明setup.bounds.phase.integral的范围。如果你把积分项手动写成端点函数里的一个积分近似值传进去协态映射关系就完全破坏了还原出来的协态必然对不上。4.4 网格自适应一直不收敛残差卡在某个区间hp 自适应方法也有自己的盲区如果真实最优控制在一个极小的时间窗口内剧烈跳变比如 bang-bang 控制在切换点附近的抖振网格细化策略可能一直递归细分也拿不到高精度。这时候不要死磕自动网格手动干预它先跑一遍画出解的残差分布找到残差集中区域然后手动在 setup.mesh.phase.colpoints 里设置较密的初始网格分布让 GPOPS II 从一个更合理的起点开始。另外提醒一句mesh 容差别设得太激进。setup.mesh.tolerance 1e-10听起来很爽但可能让你的求解器陷入无限网格细化的死循环我一般从 1e-4 试起逐步收紧到 1e-6够用了再往高调。5. 常用工具与辅助手段GPOPS II 的生态相对封闭不是那种社区文档满天飞的开源项目所以掌握一些周边辅助工具能让你少走弯路。5.1 学会用 SNOPT 的 options 调收敛性SNOPT 有相当多的选项可以调最常用的是Major feasibility tolerance、Major optimality tolerance和Minor feasibility tolerance。我的一般做法是先把三个容差都放宽到 1e-4 级别保证能够快速找到一个粗解确认问题建模无误后再把容差逐步收紧到 1e-6并配合 GPOPS II 的 mesh 细化一起提高精度。如果你一开始就上严格容差很多问题会因为初始迭代就卡住而让你误判为“模型不收敛”。5.2 用辅助数值验证解的正确性任何最优控制求解器的结果都不应该盲目信。拿到 GPOPS II 输出的解后我的习惯是把它“降回”连续时间域做一次正向积分验证把解出来的控制轨迹存成时间序列然后在高精度 ODE 求解器比如 MATLAB 的 ode45 或 ode15s里重新积分动力学方程看终端状态是否与 GPOPS II 声称的一致。这一步能抓出很多数值解表面上收敛、实际不满足动力学的问题。还可以做一次必要的“一阶最优性验证”如果问题足够简单手推或者用符号计算工具算出最优控制的结构比如 bang-bang 结构然后和 GPOPS II 的输出做对比。对于复杂问题可以检查控制轨迹是否满足庞特里亚金极小值原理的必要条件虽然这一步比较费劲但在发论文或做工程决策前很值得做。5.3 环境兼容性的一些现实提醒GPOPS II 是 MATLAB 平台的工具包不同 MATLAB 版本下的兼容性有细微差别。我踩过的坑是在 R2021a 之后的版本里某些内置稀疏矩阵求解函数的行为发生了变化导致 GPOPS II 在求稀疏雅可比时报错。这里有个偏方更新到较新版本的 GPOPS II官方 GitHub 上维护的版本如果不行尝试把 MATLAB 的optimoptions配置从默认改回 legacy 模式。另外GPOPS II 的 license 模式和文件路径管理也比较特殊建议用addpath显式添加而非依赖startup.m自动加载方便排查路径冲突。6. 一条绕不开的建议最后还是想说点实际的。GPOPS II 有门槛但它的门槛不在操作而在理解最优控制问题本身。你能不能用它解出漂亮的结果取决的不是你会不会点按钮而是你能不能把工程问题抽象成一个规范的最优控制数学模型能不能鉴别出那些看起来收敛但实际不合理的解。我的建议是如果你刚接触这个工具不要直接拿自己的复杂工程问题开刀先用一个解析解已知的经典问题比如月球软着陆、车辆换道、双积分系统时间最优控制走通全流程再把工具替换成你自己的问题模型。这个“带参考答案练习”的过程能帮你省下至少一个月的试错时间。等你在两三个经典问题上把 GPOPS II 的脾性摸透了再回到自己的项目里你会明显感觉到这工具为什么能被称为“史上最牛”之一——它把最难的离散化和网格适配问题藏在了底层而把最需要物理直觉的部分留给了你。这种取舍恰恰是它最成熟的地方。本文还有配套的精品资源点击获取