简介：面向航空航天、机械工程等领域的科研与工程人员，开源最优控制求解器 GPOPS II 的完整程序包以高斯伪谱法为核心，将连续时间最优控制问题离散化，适合处理多阶段约束、多状态变量的轨迹优化与控制策略求解。包内共 293 个文件，压缩包为 rar 格式，约 10.57MB；主体为 204 个 Matlab 源程序，另有 pdf/tex 文档、eps 图、mat 数据以及多种平台 mex 编译文件，便于跨环境运行和对照原理阅读。已有 4485 人学习下载。使用者可从中获得高斯节点插值、状态与控制参数化、Collocation 优化及后处理等完整代码，以及多阶段问题设定、状态轨迹绘制等辅助脚本和内置导数计算模块；配合文献图表理解 IPOPT 等非线性规划求解的调用方式，便于迁移到火箭、飞行器轨迹优化等实际项目中。
1. GPOPS II 到底解决什么问题：最优控制里的“把微分方程直接丢给优化器”
“史上最牛逼”这种话听听就好，但高斯伪普法优化程序 GPOPS II 这套工具本身，确实把一类让人头疼的问题变成了可以反复跑的任务。做机器人时间最优轨迹、飞行器入轨、车辆能耗最优、机械臂动作规划时，传统做法是手写间接法去猜协态初值、猜切换时间，猜错一次就要重推一轮；GPOPS II 的做法是直接把你写好的状态方程、边界条件和路径约束打包成一个大规模稀疏非线性规划，配点、插值、网格细化全自动完成，你只需要关心物理模型建得对不对、初值给得合不合理。适合的读者很明确：在 MATLAB 里做最优控制或动态优化的人，以及被终端边界和路径约束反复折磨的菜鸟与老手。
2. 高斯伪谱法为什么是主流方案：从打靶法到配点法的取舍
2.1 你写不出协态初值时，直接法是唯一活路
最优控制问题教科书里的经典解法是间接法：先写哈密顿函数，对状态求偏导得到协态微分方程，再结合横截条件把问题变成两点边值问题。逻辑很漂亮，工程上却很现实——协态初值没有物理意义，你根本不知道它大概在什么量级，只能靠猜。遇到 bang-bang 控制还要先猜切换结构，遇到奇异弧更是雪上加霜。我见过有开发者在一个带积分约束的问题上卡了两周，最后发现是协态初值给错了方向。
直接法绕开协态猜测的思路是：把连续时间最优控制问题在时间轴上取一批配点，用插值多项式近似状态和控制，把微分方程改成配点处的代数方程，目标函数和约束也全部离散成代数形式。这样求出来的就是数值解，虽然没有闭式解析表达式，但对工程足够。高斯伪谱法就是直接法的一种具体实现，GPOPS II 是这类实现里最成体系的一套。
2.2 LG 和 LGR 配点：伪谱法把配点放在哪、凭什么准
普通直接配点法喜欢用均匀网格，优点是代码简单，缺点是精度低。要让结果准，就得不断加密网格，问题规模随之膨胀。伪谱法的思路是换一种配点分布：用正交多项式的根作为配点。GPOPS II 内部用的是 Legendre-Gauss-Radau 配点，简称 LGR，配点分布在区间内部以及终止端。
这套设计有两个直接好处。第一，终端状态就是配点上的值，不需要像很多方法那样用外推公式去估计终点，最终的边界约束可以直接写进离散方程组；第二，积分运算能用高斯求积公式的权重精确处理，目标函数里的积分约束和代价项不会白白引入额外误差。很多刚接触的读者会把这个当成黑匣子，其实可以在你自己的问题里做对比实验：同一个问题分别用均匀配点和 LGR 配点跑，LGR 通常用更少的点就能达到同样的精度。
2.3 h 与 p 的博弈：hp 自适应网格细化的基本盘
GPOPS II 的网格细化策略叫 hp 自适应。每个阶段求解完毕，它会估计当前网格下的误差，再决定下一步怎么改。h 是指把时间区间切成更多段，每段长度更短；p 是指增加某一段的配点数量，提高该段插值多项式的阶数。两者都能提高精度，但代价不同，实际效果也随问题性质变化。
GPOPS II 会根据误差分布决定哪些段该 h、哪些段该 p。你不需要手动干预，但一定要会读它打出来的迭代日志。日志里一般会显示每一段的最大误差估计值，如果某一段反复出现相同量级的误差却迟迟不收敛，多半是动力学在那一段有非光滑变化，或者是函数里有不可导操作，后面避坑章我再细说。网格配置里有两个参数直接决定这一步的走向：网格容差和最大迭代次数，一动一静之间就是收敛与翻车的距离。
2.4 GPOPS II 与手写打靶相比，多付出的成本在哪
直接法并非没有代价。状态、控制、路径约束全部离散成配点上的变量后，NLP 规模比打靶法大得多。比如一个 10 段网格、每段 8 个配点的两状态问题，光状态变量就有 160 个，加上控制、时间、参数，量级会涨到几百甚至上千。好在这个规模的稀疏 NLP 已经有成熟的 IPOPT、SNOPT 这类求解器兜底，GPOPS II 帮你生成雅可比矩阵、约束稀疏结构，省掉的推导时间远大于求解器多花的那几秒。
从我的使用习惯看，手写打靶法的上限是两到三个状态、无路径约束的简单问题；一旦加入路径约束、多段、积分约束，直接法就是唯一能在一晚上给出结果的选择。这也是为什么很多做轨迹规划的团队最终都会转向这类配点法工具——不是因为它数学上更华丽，而是因为它让你把精力留在建模上。
3. 用 GPOPS II 跑通最短时间问题：建模、离散到结果提取的完整流程
3.1 问题定式：一维小车最短时间到达
用一个经典例子把流程串起来：一维小车，状态是位置 x 和速度 v，控制是推力 u，满足 x 对时间的导数等于 v，v 的导数等于 u。小车从静止出发，要到达位置 x=1，终点速度必须是 0，目标是最短时间，控制约束在 [-1, 1]。这个问题有解析解，最优控制是 bang-bang：先油门全开加速，再油门全开减速，非常适合用来验证 GPOPS II 的建模和求解是否正确。
写代码之前先说一个常见误区：不要把连续时间问题空格离散写进去。GPOPS II 的处理方式是先给一个初始网格，再用自适应细化逼近连续解。你需要做的是把问题说清楚，也就是写三个东西：状态方程和路径约束、目标函数、边界条件。这里的路径约束就是控制 u 的限幅，目标函数是最终时间，边界条件是首末状态固定。
3.2 写 continuous 和 endpoint 函数：GPOPS II 的主战场
GPOPS II 的问题定义集中在两个回调函数上。第一个是 continuous 函数，在配点处被调用，返回状态导数和路径约束值。第二个是 endpoint 函数，在整段积分结束后被调用，返回目标函数和事件约束。以下是最小时间小车的 continuous 函数：
function [statesdot, path] = cart_continuous(input) state = input.phase(1).state; % N x 2 矩阵，第一列是 x，第二列是 v control = input.phase(1).control; % N x 1 矩阵，推力 u x = state(:, 1); v = state(:, 2); u = control(:, 1); % 状态方程：dx/dt = v, dv/dt = u statesdot = [v, u]; % 路径约束：u 的上下界，这里直接返回 u 本身 path = u; end常规逻辑是：statesdot 的列数必须和状态数一致，列顺序要和 state 的列顺序一致。path 的列数对应路径约束的数量，约束的实际上下界不在函数里写，而是在 setup 的 bounds 里声明。这里把 u 返回给 path，然后在 bounds 里约束 u 在 [-1,1]，比在函数里手动截断 u 正确得多——后者会引入不可导操作，解出来的不是最优控制。
endpoint 函数负责目标函数和事件约束。终点位置和终点速度这类边界条件，可以直接写进 bounds 的 finalstate，也可以放进 endpoint 里当事件约束。我一般优先用 bounds，因为更直观，也少写一组 eventgroup 的上下界。下面的 endpoint 只返回目标：
function [objective, eventg] = cart_endpoint(input) objective = input.phase(1).finaltime; % 最小化最终时间 eventg = []; % 事件约束为空，终态交给 bounds end如果你的 GPOPS II 版本不接受空的 eventg，就把终态放进事件约束组，用 eventgroup 的 lower 和 upper 限定，两种方式结果一致。区别是 bounds.finalstate 的实现更轻量，事件约束则适合多段衔接的场景。新手建议先用前者，少一个出错点。
3.3 定义 bounds、初始猜测并调用求解器：一个完整的 setup 脚本
下面这一段是主脚本，把函数句柄、边界、初始猜测、网格和 NLP 配置串起来：
setup = struct; setup.name = 'cart_mintime'; setup.functions.continuous = @cart_continuous; setup.functions.endpoint = @cart_endpoint; % 时间边界：起点固定 0，终点要优化，给一个合理范围 setup.bounds.phase(1).initialtime.lower = 0; setup.bounds.phase(1).initialtime.upper = 0; setup.bounds.phase(1).finaltime.lower = 0.01; setup.bounds.phase(1).finaltime.upper = 10; % 状态边界：位置和速度的物理范围 setup.bounds.phase(1).initialstate.lower = [0; 0]; setup.bounds.phase(1).initialstate.upper = [0; 0]; setup.bounds.phase(1).finalstate.lower = [1; 0]; setup.bounds.phase(1).finalstate.upper = [1; 0]; setup.bounds.phase(1).state.lower = [-0.5; -2]; setup.bounds.phase(1).state.upper = [1.5; 2]; % 控制边界 setup.bounds.phase(1).control.lower = -1; setup.bounds.phase(1).control.upper = 1; % 初始猜测：NLP 需要一个非零的出发点 nGuess = 20; setup.guess.phase(1).time = linspace(0, 1, nGuess)'; setup.guess.phase(1).state = [linspace(0, 1, nGuess)', zeros(nGuess, 1)]; setup.guess.phase(1).control = 0.2 * ones(nGuess, 1); % 网格设置 setup.mesh.method = 'hp'; setup.mesh.tolerance = 1e-6; setup.mesh.maxiteration = 15; setup.mesh.colpoints = 4 * ones(1, 10); % NLP 设置 setup.nlp.solver = 'ipopt'; setup.nlp.tolerance = 1e-7; % 求解 sol = gpops2(setup);这里有两处容易踩坑。第一，initialtime 的上下界都写 0，表示时间起点固定，自由的是最终时间；如果只想优化时间区间，上下界可以给不同值。第二，初始猜测的时间轴只有 0 到 1，而 finaltime 上界给的是 10，两者不一致没关系，NLP 求解器会自己调整网格时间分布，但初始猜测不能给全零数组，否则雅可比矩阵容易奇异。
求解完成后，结果都存在 sol.solution 里。取数据用下面的方式：
phase = sol.solution.phase(1); t = phase.time; % 配点时间 x = phase.state(:, 1); % 位置 v = phase.state(:, 2); % 速度 u = phase.control(:, 1); % 最优控制 lambda = phase.costate; % 协态，用于校验，见第 6 章 figure; plot(t, x, 'b-', t, v, 'r--', t, u, 'k:'); legend('x(t)', 'v(t)', 'u(t)'); xlabel('time');这个例子的解析解是 t_final=2，u 先 +1 再 -1，位置单调上升到 1。如果你跑出来的 max 位置超过 1，不用慌，只要终态落在 x=1、v=0 就是可行解，过冲程度和网格细化精度有关。把网格 tolerance 从 1e-6 放到 1e-8 再跑一轮，轨迹会更贴解析解，这也是后文验证网格设置是否合理的最直观方法。
4. 收敛全靠细节：GPOPS II 必调的 5 个参数与初始猜测策略
4.1 参数总表：从 NLP 求解器到网格容差
GPOPS II 的默认配置能解决简单问题，但真正棘手的问题，参数调不好就是死活不收敛。我总结了五个优先级最高的参数，按推荐调整顺序列成下表：
| 参数 | 位置 | 常见坑 | 推荐起点 |
|---|---|---|---|
| NLP 求解器 | setup.nlp.solver | 默认解不出就换 | ipopt，稀疏问题可试 snopt |
| NLP 容差 | setup.nlp.tolerance | 设到 1e-9 容易卡 | 1e-6 起步，验证后再收紧 |
| 导数模式 | setup.nlp.derivatives | 复杂函数自动微分易崩 | 先 numeric，再 automatic |
| 网格容差 | setup.mesh.tolerance | 过严导致网格迭代爆炸 | 1e-6 起步 |
| 网格最大迭代 | setup.mesh.maxiteration | 过小会提前终止 | 10 到 15 |
IPOPT 和 SNOPT 都是成熟求解器，没有绝对优劣。IPOPT 对内存占用更温和，适合普通电脑；SNOPT 在大型稀疏问题上通常更稳，但对可行性容差的处理更严格。如果你的问题在 IPOPT 下卡在约束违反度下不来，换 SNOPT 重新跑同一份 setup 是常见做法。
4.2 变量缩放：状态、控制、时间三个维度的“一个数量级”
经验法则是让参与优化的所有变量名义值落在 0.1 到 10 之间。GPOPS II 对时间会做归一化处理，但状态和控制的缩放完全靠你自己。以飞行器轨迹优化为例，高度可能是几万米，速度可能是每秒几千米，质量可能是几百公斤，三个变量差三四个数量级，如果不缩放，NLP 的雅可比矩阵条件数会非常难看。
我一般不在 bounds 里手动换算，而是在建模函数里除以标称值，例如把速度写成 v = state(:,2) * vScale，然后把 v 的 bounds 设置成无标称值的倍率。这样力学方程里的系数也跟着缩放，方程形式变了但物理不变。另一个量纲是控制力，如果控制上限本身有物理含义，先除上限再约束到 [-1,1]，会让梯度数值更平衡。记住一个原则：宁可花十分钟做缩放，不要省这几分钟然后对着不收敛日志发呆一下午。
4.3 初始猜测策略：先用低配点加宽松容差跑热启动
很多新手一上来就设 maxiteration=50、网格 tolerance=1e-9，结果卡在网格细化阶段反复迭代，看起来像死循环。正确做法是分两步走。第一步，把网格 tolerance 放到 1e-3，colpoints 用小规模 3 或 4，先跑出一个粗解。第二步，把粗解作为初始猜测，带着更紧的 tolerance 再跑一轮。
GPOPS II 支持把上一轮的解直接塞进 guess 字段，操作是读取 sol.solution.phase(1) 的 time、state、control，赋给 setup.guess.phase(1)。这个热启动路径能省掉大量冷启动的随机搜索。对 bang-bang 问题，粗解往往已经能看到控制曲线的开关结构，你再决定要不要手动把猜测改成两段常值来加速收敛。初始猜测别给零，但也别给精确解——给精确解会让求解器觉得问题太简单，反而掩盖了建模里的错误。
4.4 网格最大迭代数与网格细化日志怎么读
网格迭代日志是 GPOPS II 最值钱的输出。每轮它都会打印当前网格各段的误差估计值，以及本轮把哪些段做了 h 加密、哪些段做了 p 升阶。读日志的关键是看误差是不是在逐轮下降，而不是看它迭代了多少轮。如果两轮之间误差完全没有变化，说明问题出在函数本身，而不是网格不够密。
maxiteration 卡住时的另一个常见原因是路径约束在配点之间悄悄超限。路径约束在连续函数里只检查配点位置，配点之间的包络线可能穿出约束边界，网格细化算法会根据插值残差捕捉到这一点并不断加密。这时你应该先检查是不是状态或控制本身存在剧烈振荡，而不是一味调大迭代上限。把 mesh.colpoints 调大试试，有时误差不降是因为某一段需要更高阶插值，而不是需要更多段。
5. GPOPS II 常见问题与避坑记录：不收敛、NaN、网格爆炸从哪查
5.1 现象：网格细化卡在同一个阶段，日志误差半天不动
原因大多是动力学里存在非光滑项。比如摩擦模型用了 sign 函数、控制器模型里夹了 if-else、或者路径约束用了不可导的取最大值函数。伪谱法的插值多项式是全局光滑的，遇到非光滑点会在该段产生吉布斯式振荡，误差始终压不下去。
解决方法是把非光滑项改成光滑近似。摩擦力用 tanh 代替 sign，速度死区用平滑过渡函数，切换逻辑用连续可导的分段函数拼接。改完之后通常一到两轮网格迭代误差就会大幅下降。这是最值得记住的一条：GPOPS II 对函数光滑性的要求，比普通打靶法苛刻得多。
5.2 现象：continuous 函数偶发 NaN，但单独测试一切正常
原因通常发生在配点插值之后。GPOPS II 在网格细化阶段会尝试新的配点分布，插值出的状态可能落在 bounds 之外，比如速度变负、质量变小、对数函数括号里出现负数。单独测试 continuous 函数时只测了你给的初值，不会覆盖所有插值路径。
解决方法是先看 NaN 发生时的配点状态，打印 input.phase(1).state 的值，找到是哪个变量越界。然后在建模函数里对越界做防御，比如对数输入加一个 floor 限制，再配合 bounds 把状态约束收窄。我不建议用 try-catch 把 NaN 吞掉，那样求解器会失去方向的指引，收敛更加随机。
5.3 现象：加入积分约束后目标值异常，最优解明显不对劲
积分约束的正确写法是把积分量通过输入提供给模型，而不是自己在状态向量里再造一个状态。GPOPS II 的 continuous 函数可以返回第三个输出 integrals，对应积分约束的导数，然后需要在 bounds 里定义 integral 的下界与上界。如果你把积分量写成了普通状态，相当于额外增加了一条微分方程，还要自己处理它的初值和末值，很容易漏掉一致性条件。
解决方式是按工具包自带多段积分示例的写法，把积分量放到 integral 维度，不要让状态维度膨胀。同时检查 bounds.phase(1).integral 是否定义全，漏掉任何一维都会导致约束数量不匹配，求解器要么报错要么给出离谱结果。
5.4 现象：自由最终时间求出个离谱的大数
自由终端时间问题里，最终时间的上界给太宽，会导致时间变量和状态变量的尺度严重失衡。另外，如果目标函数对最终时间的变化不敏感，例如最终时间出现在动力学方程里但目标函数里的权重很小，求解器可能就直接接受一个大时间。
解决方案是给最终时间一个物理上能接受的窄范围。以小车最短时间为例，解析解是 2，我给上界 10 已经够宽松；如果某个问题物理上只能跑 1 秒左右，上界给 1000 就是给自己找麻烦。先固定最终时间求解一次，再用得到的状态轨迹作为自由时间问题的初始猜测，也是一种有效收敛手段。
5.5 现象：换了 NLP 求解器，最优解变了
同一份 setup 用 IPOPT 和 SNOPT 跑出不同结果，不一定是模型错了，可能只是求解器的可行性容差和迭代终止条件不同。IPOPT 倾向于在目标函数方向上走得更深，SNOPT 更重视可行性，两者对约束违反度的判断标准也有细微差异。
解决方法是先固定同一容差，再对比两者结果的约束违反度和目标函数值。如果目标函数值差很多但约束都满足，大概率是多解或奇异性问题，不是求解器 bug。每次只改一个变量，是调参的基本纪律，同时把两个求解器的日志都存下来，方便对比终止原因。
6. 用伴随变量验证最优性：GPOPS II 结果的最后一道校验
6.1 协状态轨迹怎么从结果里取
GPOPS II 在求解过程中会生成完整的协态变量估计，保存在结果结构的 costate 字段中。取法和状态一样：
phase = sol.solution.phase(1); t = phase.time; u = phase.control(:, 1); costate = phase.costate; lambda_v = costate(:, 2); % 第二列对应速度分量直接画 u(t) 和 lambda_v(t) 两条曲线，是判断解是否满足最优性条件最快捷的方式。对小车最短时间这类 bang-bang 问题，最优条件是控制方向由协态符号决定：u=+1 时对应速度协态为负，u=-1 时对应速度协态为正。看到符号匹配，再结合终点状态约束，心里就踏实了。
6.2 对平滑问题数值校验 dH/du 是否接近零
对没有控制饱和或饱和区间较短的平滑问题，最优控制应该让哈密顿函数对控制的梯度接近零。你不需要手推协态微分方程，直接用数值梯度算 H 对 u 的值即可：
% 重构哈密顿函数 H = L + lambda' * f，这里 L=0，f=[v; u] H = lambda_v .* u; % 控制相关项，其余项为常数 dH_du = gradient(H, t); % 数值近似 dH/du 的趋势由于 GPOPS II 的协态符号约定可能与教材差一个负号，第一次使用前先拿自带的 bang-bang 示例对比一次符号方向，确认约定后再应用到自己的问题。校验偏差较大时，优先怀疑网格不够细，其次是协态插值在边界配点处的精度损失。三分钟内完成这一步，能过滤掉大半粗心错误。
6.3 我的使用习惯：校验脚本固化下来
不管问题是 bang-bang 还是平滑最优，我现在跑完 GPOPS II 都会把协态和控制画在同一张图上。数值上我不会追求严格为零，而是看趋势和符号是否合理。协态是分段的就让它在段内连续，段间有跳跃是正常的，反而说明多段设置是对的。
这套方法的可取之处在于它不依赖解析解，也不要求你手工重推 KKT 条件，只是把求解器已经算好的信息用起来。我常在怀疑自己建模出错时先看这张图——如果协态曲线出现无法解释的毛刺，多数情况下是建模函数在某个切换点附近不光滑，回到第 5 章的平滑策略处理一遍，比重新猜初值有效得多。希望帮到你。
本文还有配套的精品资源，点击获取
