简介这是一份面向最优控制与轨迹优化研究的高斯伪谱法实现——GPOPS II 完整程序包适合航空航天、机械工程等领域需要求解非线性动态系统最优控制问题的工程师与科研人员。压缩包共293个文件、约10.57MB包含204个Matlab源码文件m、41个矢量图文件eps、19个PDF文档及LaTeX源文件、多个平台的mex可执行文件等代码与文档结构清晰便于查阅和二次开发。已有4478人浏览学习是理解高斯伪谱法原理与工程落地的实用参考。资源提供从问题设定、高斯节点插值、状态与控制参数化到非线性规划求解的完整实现并支持多阶段问题与内置导数计算。通过研读源码和示例可快速掌握IPOPT等求解器的调用方式并应用于飞行器轨迹优化、能源系统控制等实际场景显著提升复杂优化问题的建模与求解效率。 我第一次在论文里看到“高斯伪谱法”五个字时其实是不太信的。那段时间我正被一个带路径约束的火箭垂直着陆轨迹优化问题折磨得够呛打靶法初值猜得我怀疑人生非线性规划每一轮迭代都像是在悬崖边走路动不动就发散。后来同一个师门的兄弟甩给我一句话你试试 GPOPS高斯伪谱法那个优化程序解这类问题基本是降维打击。然后我就开始了和 GPOPS II 打交道的几年。这个标题“史上最牛逼的高斯伪普法优化程序 GPOPS II”确实起得浮夸但如果你真在航空航天、机器人运动规划、过程控制这些领域里用最优控制做过工程大概率会认同它的地位——至少在开源/学术可获取的范围内GPOPS II 是把“高斯伪谱法”从论文里抠出来、变成能直接求解工程问题的最成熟工具。它本质上是一个基于 MATLAB 的最优控制求解框架核心思路是把你头痛的连续时间最优控制问题通过高斯伪谱配点离散成一个大规模稀疏非线性规划NLP然后交给 SNOPT 或 IPOPT 这类求解器去收敛。这篇文章不打算讲太抽象的理论重点说清楚它为什么好用、哪些环节最容易翻车、以及我实际跑通一个燃料最优着陆问题的完整过程希望能帮你少走点弯路。1. 伪谱法为什么能在最优控制里站稳脚跟1.1 从打靶法到配点法的思路变迁我刚接触最优控制时教材里讲得最多的还是间接法和打靶法。间接法就是推导哈密顿函数、协态方程、横截条件最后解一个边值问题。这个方法数学上很漂亮但工程上非常脆只要目标函数、动力学或者约束稍微复杂一点解析求导就成了噩梦边界条件稍微给得不合适打靶的初值稍微偏一点整个 shooting 过程就会出现剧烈的数值振荡。后来接触直接法思路完全反过来——不碰协态不推导最优性必要条件直接把状态和控制离散成一段一段的变量让优化器去搜一组能满足动力学和约束、同时让目标函数最小的离散点。这里面的关键问题就变成了怎么离散用欧拉法离散精度太差需要把时间轴切得很碎变量维度一下子涨到几万NLP 求解器也扛不住。用高阶龙格库塔法离散精度上去了但每一步的雅可比矩阵算起来很麻烦且离散点分布并不总能捕捉到轨迹里那些变化剧烈的位置。伪谱法走的是另一条路把状态和控制都用全局插值多项式来逼近在特定的高斯配点上满足动力学约束。因为高斯求积公式的精度极高少量配点就能达到很高的离散精度。换句话说伪谱法用“少而精”的离散点换来整体高精度这是它能压缩问题规模的关键。1.2 高斯伪谱法的核心把连续问题离散成 NLP高斯伪谱法里最典型的做法是 Legendre-Gauss 配点。它把时间区间映射到 [-1, 1]然后在整个区间上用拉格朗日插值多项式来近似状态和控制要求动力学残差在配点处等于零。因为这个过程会把微分方程约束变成一个代数方程组原问题就变成了一个标准的 NLP优化变量配点处的状态值、控制值、初始/终端时间以及可能的静态参数等式约束配点处的动力学残差为零、初始/终端状态与状态变量的连接关系、事件约束不等式约束路径约束如过载、动压限制和控制幅值限制目标函数通常是终端状态函数或积分型指标积分部分用高斯求积近似。GPOPS II 内部做了大量细节处理比如把多个时间段的网格Mesh拼接起来、用高斯求积计算积分约束、自动生成稀疏雅可比矩阵然后输出成 NLP 求解器能吃的格式。你不需要自己去手写配点方程也不需要手动推导导数只需要按照它的接口把动力学函数、端点条件、边界范围写清楚剩下的交给框架去装配。这里我想强调一个容易被忽视的点伪谱法的强项是“近似全局最优解的速度”但它的解是离散配点意义上的最优不是连续解析意义上的严格最优。工程上我们通常会先用 GPOPS II 得到一个精度合格的参考轨迹再把它作为初值喂给更精细的优化器或者做闭环跟踪。理解这个定位你就不会对它的某些数值现象产生误判。2. GPOPS 与 GPOPS II 的进化差异2.1 GPOPS 和 GPOPS II 不是换皮关系很多新手会把 GPOPS 和 GPOPS II 混为一谈觉得后者只是修了几个 Bug 的升级版。实际上 GPOPS II 在算法架构上是重新设计的。第一代 GPOPS 基本是固定配点数量的全局伪谱法一旦整个时间域上的轨迹变化剧烈比如存在 bang-bang 控制或者急转弯固定配点就很容易出现振荡你只能手动加配点数量然后让求解器硬扛。GPOPS II 引入了 hp 自适应网格细化——把时间域分成多个网格段每个段内用伪谱配点优化完一轮后评估误差误差大的网格段会自动加密配点或者细分网格下一轮重新优化。这个机制大大提高了处理非光滑轨迹的能力。2.2 hp 自适应网格细化精度与算力的平衡hp 自适应里的 h 指网格段的尺寸p 指配点阶数两者配合使用达到“把计算量花在刀刃上”的效果。GPOPS II 默认的网格细化方法有 hp-PattersonRao 和 hp-LiuRao区别在于估算误差和决定加密策略的细节。实际使用中我大部分时间就用默认的 hp-PattersonRao它在一阶不连续点附近能自动识别出需要加密网格的位置对 bang-bang 控制轨迹很友好。需要注意的是网格细化不是免费的。每细化一轮就要重新求解一次 NLP如果初始网格给得太差可能要迭代十几轮才能达到设定的容差总时间反而比用更密的初始网格慢。工程上有个土办法先用宽松的setup.mesh.tolerance 1e-3快速跑通一条可行轨迹看看轨迹形状再根据形状手动指定初始网格段和配点数最后收紧容差到 1e-6 做精细求解。一上来就设 1e-8除了让求解器疯狂细化网格、不断报“网格迭代次数超限”之外没有任何好处。2.3 ADiGator 自动微分摆脱手推导数的噩梦GPOPS II 另一个重大升级是集成了 ADiGator 自动微分工具。早期用直接法做轨迹优化最痛苦的环节之一就是给 NLP 求解器提供导数解析推导容易错有限差分算起来慢而且差分步长选不好还会引入数值噪声。ADiGator 会对你的 MATLAB 连续函数和端点函数做源码级别的自动微分生成精确的一阶和二阶导数代码然后直接提供给 NLP 求解器使用。这件事听起来是“自动化”但实际使用有一个关键前提你的动力学函数里不能有非光滑操作。比如abs、floor、ceil、interp1这类函数自动微分处理起来非常麻烦轻则生成错误的导数重则直接报错。做工程问题建模时凡是涉及查表或取整的逻辑尽量改写成光滑近似或者用边界约束去规避否则后患无穷。我第一次连续函数里写了个min(max(...))做推力限幅ADiGator 倒是没报错但求解出来的轨迹一看就是假的后来排查了半天才发现是导数不对。3. GPOPS II 工程架构与文件组织实战3.1 setup 阶段把问题“翻译”成机器能懂的结构GPOPS II 的用法比大多数类似工具要规整一个完整的求解流程通常分三步写 setup 结构体、写连续函数、写端点函数。setup 阶段相当于把问题“翻译”成优化器能读懂的配置表单它告诉 GPOPS II 状态有几个、控制有几个、边界在哪里、用什么求解器、用什么导数模式、网格细化怎么控制。下面是我常用的一套最小配置骨架setup.name VTVL_FuelOptimal; setup.functions.continuous vtvlContinuous; setup.functions.endpoint vtvlEndpoint; setup.auxdata struct(...); % 附加常量参数如 g0, Isp 等 setup.bounds.phase.initialtime.lower 0; setup.bounds.phase.initialtime.upper 0; setup.bounds.phase.finaltime.lower 0; setup.bounds.phase.finaltime.upper 300; setup.bounds.phase.initialstate.lower [1; 0; 1]; setup.bounds.phase.initialstate.upper [1; 0; 1]; setup.bounds.phase.state.lower [0; -1.5; 0.4]; setup.bounds.phase.state.upper [1.2; 1.5; 1]; setup.bounds.phase.finalstate.lower [0; 0; 0.4]; setup.bounds.phase.finalstate.upper [0; 0; 1]; setup.bounds.phase.control.lower 0.1; setup.bounds.phase.control.upper 4; setup.derivatives.supplier adigator; setup.derivatives.derivativelevel second; setup.nlp.solver ipopt; setup.mesh.method hp-PattersonRao; setup.mesh.tolerance 1e-6; setup.mesh.maxiterations 20; setup.mesh.colpointsmin 4; setup.mesh.colpointsmax 15;一个新手很容易忽略的地方是setup.auxdata。它像是一个只读的全局变量容器专门用来传递不参与优化的常量。把它用好有两个好处一是方便批量跑不同参数的算例二是避免在连续函数里写死常量导致后续改工况时到处翻代码。我的习惯是把所有物理常量、模型参数都塞进 auxdata连续函数里需要什么就取什么。3.2 连续函数与端点函数两个回调函数的职责划分连续函数负责定义“微分方程动力学 路径约束 积分被积函数”它在每个配点处被调用。输入是你当前时间网格上的时间input.phase.time、状态input.phase.state、控制input.phase.control输出则是动力学导数output.dynamics、路径约束值output.path和积分被积函数output.integrand。端点函数负责定义“目标函数 事件约束”它只在初始端和终端被调用输入是初始/终端时间和状态以及参数输出是目标值output.objective和事件约束output.eventgroup.group。这里最容易搞混的是初始状态约束和终态状态约束可以直接写在setup.bounds.phase.initialstate和finalstate里不一定非要通过端点函数的事件约束去处理。事件约束适合那些没法直接用状态上下界表达的复杂条件比如“终端质量不能低于初始质量的 40%”“两个阶段之间的某个状态必须平滑连接”这类跨阶段约束。从功能划分上理解连续函数对应系统的“物理规则”端点函数对应“任务目标与最终要求”。物理规则写错了后面再怎么调参数都是白搭任务目标写歪了优化器会非常诚实地给你一个满足约束但完全没用的“最优解”。3.3 从 GPOPS 结构到 NLP 求解器的调用链路GPOPS II 内部并不是直接调用你写的连续函数去求导而是先把你提供的函数和 setup 配置打包成一个内部结构然后由它的核心驱动函数完成三件事将时间域划分成网格段并在每个网格段内生成配点、在配点上计算动力学残差和约束值、通过 ADiGator 生成整个 NLP 问题的雅可比矩阵和海森矩阵如果导数级别设成 second 的话。最后它把稀疏矩阵形式的 NLP 模型输出给 SNOPT 或 IPOPT 求解。这意味着你写的 MATLAB 函数会被自动微分工具“扫描”很多遍所以代码风格会影响求解稳定性。尽量写向量化表达式避免在函数内部搞复杂的循环嵌套文件名和函数名保持一致不要用全局变量传参都用input.auxdata传。这些习惯能让 ADiGator 生成的导数代码更干净也能减少魔改时踩坑的概率。4. 一个可复现的垂直着陆燃料最优问题4.1 问题建模与无量纲化处理这里用一个简化版的垂直起降着陆问题来说明完整流程。假设飞行器在竖直平面内做一维运动状态取无量纲高度 h、无量纲速度 v、无量纲质量 m控制取无量纲推力 T重力加速度设为 1比冲折算成一个无量纲常数。目标是最小化燃料消耗等价于最大化终端质量即目标函数写为负的终端质量。无量纲化这一步看起来多此一举但对求解器来说极其关键。SNOPT 和 IPOPT 这类 NLP 求解器对变量尺度非常敏感。如果高度用米、速度用米每秒、质量用千克量级可能差到 10^5 以上求解器在计算 KKT 条件时容易出现数值病态。无量纲化后所有变量基本都在 0 到几个单位之间收敛速度和稳定性都会有肉眼可见的提升。4.2 setup 和函数的完整代码连续函数定义为function output vtvlContinuous(input) h input.phase.state(:, 1); v input.phase.state(:, 2); m input.phase.state(:, 3); T input.phase.control(:, 1); aux input.auxdata; g aux.g; IspEff aux.IspEff; hdot v; vdot T ./ m - g; mdot -T ./ IspEff; output.dynamics [hdot, vdot, mdot]; output.path T ./ m - aux.amax; % 过载不超过 amax路径约束 output.integrand T; % 记录推力积分可用于附加约束 end端点函数定义为function output vtvlEndpoint(input) mf input.phase.finalstate(3); tf input.phase.finaltime; output.objective -mf; % 最大化终端质量 output.eventgroup.group [mf; tf]; % 事件约束终端质量下限、着陆时间上限 end对应的 setup 里要增加事件约束的边界setup.bounds.eventgroup.group.lower [0.4; 0]; setup.bounds.eventgroup.group.upper [1; 200];初始猜测可以给得很粗糙但不能完全违背物理。比如setup.guess.phase.time [0; 100]; setup.guess.phase.state [1, 0, 1; 0, -0.8, 0.5]; setup.guess.phase.control [1; 1];这里有一个很容易踩的坑初始猜测里的终端速度如果直接给 0且中间没有任何过渡伪谱多项式在边界附近会产生剧烈的龙格现象导致初始迭代崩溃。给一个“大致在减速但还没完全刹停”的猜测反而更容易收敛。我的经验是初值曲线宁可粗糙但整体趋势合理不要人为制造一个与动力学明显冲突的“精确猜测”。4.3 网格细化过程中的误差追踪求解完成后GPOPS II 会把结果放在result.solution和result.solver等字段里。我每次做完优化都会看一下result.meshhistory——它记录了每一轮网格细化时的最大误差估计、网格段数量和配点数量。正常情况下误差应该逐轮下降网格段和配点数在最开始几轮增长后期趋于平稳。如果发现误差在某个网格段反复横跳说明那个区域可能存在不连续可以手动把该区域附近的初始网格段切得更细。比如你在setup.mesh.phase里可以这样初始化网格setup.mesh.phase.fraction [0, 0.3, 0.7, 1]; setup.mesh.phase.colpoints [6, 8, 8];这表示将时间域分成三段每段分别用 6、8、8 个配点。给 GPOPS II 一个分布合理的初始网格会让它在第一轮迭代时就得到一个质量不错的解网格细化迭代次数能少一半以上。5. 边界条件、路径约束与积分量最容易翻车的地方5.1 initialstate/finalstate/eventgroup 的语义区别这三类东西看起来都是“给状态加限制”但语义完全不同。initialstate和finalstate是直接对状态向量的初末值给出上下界适合最简单的情况比如“初始高度为 1、初速为 0、终端速度和高度都必须为 0”。但工程上经常还有更复杂的要求比如“终端质量不小于某个值”或“终端高度和速度之间有线性关系”这时候就没法简单设置上下界了需要把finalstate的范围放开然后在端点函数里通过eventgroup.group写自定义约束。我见过很多人把终端状态上下界写得很死同时又在 endpoint 里重复加同类约束结果导致约束冗余、雅可比矩阵奇异。约束越多越好的想法在这里行不通NLP 求解器需要约束保证线性无关冗余约束轻则降低收敛速度重则直接让求解器退出。保持约束的精简和独立是我在这个项目里学到的很重要一课。5.2 路径约束的缩放与初值可行性路径约束是伪谱法最容易踩坑的环节。以过载约束为例直接写T ./ m - aux.amax 0在数学上没问题但如果 T 的量级很大约束函数值就会非常大求解器在计算约束违反度时可能出现大数吃小数的现象。建议在路径约束函数内部做归一化处理让约束输出和约束边界保持在差不多的量级比如改成T ./ (m * aux.Tmax) - 1。另外路径约束的处理方式实际上是配点采样两个配点之间的轨迹可能轻微违反约束这一点很多新手不知道。GPOPS II 的网格细化会在检测到约束违反时加密网格但如果你初始猜测离可行域太远第一轮优化就可能收敛到一个局部不可行解后面网格细化也救不回来。5.3 积分约束与目标函数的两种写法GPOPS II 里积分量有两种用途放进目标函数或放进事件约束。比如最小化“总冲量”这种指标就可以在连续函数里定义output.integrand T然后在端点函数里取input.phase.integral作为目标函数的一部分。这个积分值是通过高斯求积计算的与配点分布精度一致。我在做任务规划时经常用积分量来限制“累计热流”或者“累计剂量”这类工程指标这时候需要注意如果你在连续函数的输出里定义了一个 integrand但没有在这个阶段的事件约束或目标函数里使用它GPOPS II 依然会去计算它这会白白增加计算量。多个积分量共用很少见但如果你真的需要GPOPS II 的 integrand 可以是多列的每一列对应一个独立的积分量在端点函数里input.phase.integral就是向量。6. 求解器选型、导数模式与常见坑位盘点6.1 SNOPT 与 IPOPT 怎么选GPOPS II 支持的 NLP 求解器有几款最常用的是 SNOPT 和 IPOPT。SNOPT 适合中等规模、约束较多的稀疏 NLP它基于可行点序列方法对不满足约束的初值容忍度不错收敛到局部最优解的成功率很高但它是商业软件需要单独的 license。IPOPT 是开源的内点法实现对大规模问题内存占用更友好缺点是有时候会跑到可行域边界附近再收敛导致解的精度看起来差一点。我的选择标准很简单学术验证、批量跑参、不想折腾 license 就用 IPOPT做工程项目、对稳健性要求更高、有 SNOPT 授权就用 SNOPT。同一个问题这两个求解器给出的最优轨迹基本一致但迭代步数和中间过程差别很大不要在一个求解器上卡住之后死磕同一个配置换另一个求解器尝试经常能直接绕过问题。6.2 NaN、网格不收敛、收敛到错误解的排查链路遇到 NaN 时我有一套固定的排查顺序先看初始猜测是否让动力学函数产生了除零或负数的开方——比如质量 m 接近 0 就会让T./m爆炸再看路径约束是否让状态变量在迭代过程中越界越界之后又反过来让动力学函数出错最后检查导数模式如果用了adigator但函数里有不支持的运算生成的雅可比矩阵可能是错误的表现出来就是优化迭代过程一直走不动或者目标函数乱跳。网格不收敛这件事先确认一下是不是setup.mesh.tolerance设得太小同时maxiterations又不够大其次检查初始网格段的分配是否合理比如最优轨迹前半段很平缓、后半段有剧烈的控制切换你应该一开始就把后半段分得更细而不是让网格细化算法从零开始慢慢试探。收敛到错误解是我认为最隐蔽也最可怕的坑。GPOPS II 是一个局部优化框架它不能保证全局最优。同一个问题换一组初始猜测可能收敛到完全不同的轨迹。我的习惯是对关键算例至少用三组差异较大的初始猜测分别求解如果三组结果的目标函数和轨迹形状差别很大就要警惕问题可能有多局部解或者约束条件设置有问题。6.3 和 DIDO、PROPT 等工具横向对比后的体会市面上同类最优控制工具里DIDO 和 PROPT 也很有名。DIDO 基于 Legendre 伪谱法使用门槛低、盒装体验好但底层实现不透明很多底层的求解细节你无法干预遇到问题就只剩换初始猜测一条路。PROPT 基于 TOMLAB 环境功能丰富对大规模过程优化支持很好但它是商用的license 价格不便宜而且 TOMLAB 的整体风格更偏向传统最优控制研究。GPOPS II 的优势在于透明度和可定制性。它的源码本身是开放的你可以在里面看到网格细化算法、配点生成逻辑、导数装配逻辑这意味着当算法行为不符合预期时你有机会去定位原因而不是只能黑盒重试。对于喜欢把算法用在非标准问题上的工程师来说这是非常宝贵的特性。当然代价是你需要花一些时间理解它的内部结构不可能像用商业软件那样双击图标就能跑通。在我实际的项目里GPOPS II 的定位一直是一个“高精度轨迹快速生成器”。它解决的是“给我一条满足所有约束、且指标足够好的参考轨迹”这个问题而不是“在嵌入式实时环境里在线滚动优化”。把这两件事分清楚之后我对它的期望值就变得很合理用起来也就不容易失望。最后再分享一个小技巧跑新问题之前先在官方自带的月球着陆示例上把你安装的 GPOPS II 环境完整跑一遍确认 ADiGator 和求解器的调用链没问题再动你自己的模型。这一步能省掉大量“明明模型没问题却怎么都收敛不了”的排查时间因为很多时候问题不是出在你的模型上而是出在环境配置和接口调用上。本文还有配套的精品资源点击获取