简介：面向需要求解最优控制问题的MATLAB用户，这份资源提供完整的GPOPS工具箱及配套示例库，覆盖最小爬升、运载火箭上升、高灵敏边界值等多类经典问题，每个案例均包含带详细中文注释的脚本、可运行的Main入口、问题描述txt及求解输出out，并配有一份pdf官方手册和docx安装方法，帮助用户跨越工具配置门槛。压缩包共55个文件，以m源码为核心，辅以txt、out、eps、docx、pdf、mat等格式，整体仅1.55MB，轻量便携。目前已有1180人学习下载，适合研究生、科研工程师及竞赛选手在轨迹优化、最优控制仿真中直接参考。通过研读brysonDenham、hyperSensitive等示例，读者可以理解直接配点法原理、SNOPT求解器调用链路、事件约束与Cost函数写法，并能将示例代码改造迁移至自己的控制问题中。目录按问题场景分模块组织，便于按需检索与对照学习。
1. GPOPS 解决的是 MATLAB 里最容易被写坏的一类优化问题
在 MATLAB 里做轨迹优化，多数人第一反应是打开 matlab优化工具箱 里的 fmincon，或者用 ode45 套打靶。小规模问题都能跑，一旦出现路径约束、状态约束、终端时间自由、控制量带硬限幅，就会陷入网格取点、初值发散、约束只在采样点满足的连环坑。GPOPS 工具箱（General Purpose Optimal Control Software）用 Radau 伪谱法把最优控制问题整体转录成非线性规划，交给 IPOPT 或 SNOPT 求解，再按误差自动做 hp 网格加密。使用者只写微分方程、目标函数和约束，离散化交给工具箱。适合飞行器制导、机器人轨迹规划、车辆与过程控制方向的工程师。压缩包里那本官方手册对字段说明最完整，版本差异造成的坑，最后靠它兜底。
2. GPOPS 建模范式：Bolza 标准形与 setup 结构体
2.1 最优控制问题先写成 Bolza 标准形
GPOPS 直接处理的是 Bolza 形式的最优控制问题：目标函数同时含终端项（Mayer）和积分项（Lagrange），动力学与约束全部写进一组标准形式的表达式里，数学上记为：
min J = E(x(t0), t0, x(tf), tf) + ∫ L(x, u, t) dt
约束包括状态方程 ẋ = f(x, u, t)、路径约束 C(x, u, t) ≤ 0、事件约束 b(x(t0), t0, x(tf), tf) = 0，再加上状态、控制、初末时间和初末状态的上下界。这套标准形的意义在于，它把"飞行器怎么飞、机器人怎么动"这类物理问题压缩成几个固定签名的函数问题。伪谱法在此基础上把时间区间归一化到 [-1, 1]，在 Legendre-Gauss-Radau（LGR）配点上用多项式逼近状态和控制，再用差分矩阵把微分方程变成配点上的代数约束。相比直接打靶，伪谱转录的精度对配点数不敏感，少量配点就能拿到较高的逼近阶；hp 自适应再根据误差把区间切小或升阶，这是它处理 bang-bang 控制和状态约束弧线的底气所在。
2.2 setup 结构体：五个字段把问题描述完整
GPOPS-II 的入口只有一个调用：output = gpops2(setup)。setup 至少需要五个部分：name、functions、bounds、guess、mesh，字段功能先看这张对照表：
| setup 字段 | 内容 | 注意点 |
|---|---|---|
| setup.functions.continuous | 导数、被积函数、路径约束 | 必须是独立 m 文件或脚本局部函数 |
| setup.functions.endpoint | 目标终端项、事件约束 | 每阶段只在首尾调用一次 |
| setup.bounds.phase | 时间、状态、控制、积分上下界 | 多段问题按 phase 下标展开 |
| setup.guess.phase | 初始猜测的时间、状态、控制 | 行数与网格无关，会自动插值 |
| setup.mesh | 方法、容差、迭代上限、配点数 | 收敛性主要靠它控制 |
一个通用骨架如下，字段注释对照压缩包里的例子读更清楚：
setup.name = 'my_problem'; % 问题名，决定日志前缀 setup.functions.continuous = @probContinuous; % 微分方程与积分被积函数 setup.functions.endpoint = @probEndpoint; % 目标与事件约束 setup.bounds.phase.initialtime.lower = 0; % 初始时间固定 setup.bounds.phase.initialtime.upper = 0; setup.bounds.phase.finaltime.lower = 0; % 终端时间留自由区间 setup.bounds.phase.finaltime.upper = 5; setup.bounds.phase.state.lower = [-10, -10]; % 状态上下界，每个状态一列 setup.bounds.phase.state.upper = [ 10, 10]; setup.bounds.phase.control.lower = -1; % 控制量限幅 setup.bounds.phase.control.upper = 1; setup.bounds.phase.integral.lower = 0; % Lagrange 代价总和的范围 setup.bounds.phase.integral.upper = 1e4; setup.guess.phase.time = [0; 2.5; 5]; % 三个时间点就够 setup.guess.phase.state = [0 0; 1 0; 2 0]; % 每行对应一个时刻的状态 setup.guess.phase.control = [0; 0; 0]; setup.mesh.method = 'hp'; % hp 自适应网格 setup.mesh.tolerance = 1e-6; setup.mesh.maxiterations = 15; setup.mesh.colpointsmin = 4; setup.mesh.colpointsmax = 10; output = gpops2(setup);代码里有几个容易被忽略的点：functions 绑定的是函数句柄，但 GPOPS 在内部要对函数做有限差分求导，所以不能写匿名函数或带参数的闭包，必须放独立文件或脚本末尾的局部函数；guess 只需要给出粗略状态轨迹，工具箱会自己映射到网格上，但 guess 不能和硬边界矛盾；integral 上界别设得太小，否则积分代价被边界截断，目标函数值会失真。
2.3 最小用例：先证明你的工具箱真的能跑
拿到解压目录后，第一件事是跑一个自己写、有解析解的最小问题：min ∫₀¹ u² dt，ẋ = u，x(0)=0，x(1)=1。理论解是 u(t) ≡ 1，目标值 J = 1。这个用例专门验证安装、路径配置和函数签名三件事，任何一步写错，结果都不会是干净的 J=1。
clear; clc; setup.name = 'smoke_test'; setup.functions.continuous = @smokeContinuous; setup.functions.endpoint = @smokeEndpoint; setup.bounds.phase.initialtime.lower = 0; setup.bounds.phase.initialtime.upper = 0; setup.bounds.phase.finaltime.lower = 1; setup.bounds.phase.finaltime.upper = 1; setup.bounds.phase.state.lower = -5; setup.bounds.phase.state.upper = 5; setup.bounds.phase.control.lower = -5; setup.bounds.phase.control.upper = 5; setup.bounds.phase.initialstate.lower = 0; setup.bounds.phase.initialstate.upper = 0; setup.bounds.phase.finalstate.lower = 1; setup.bounds.phase.finalstate.upper = 1; setup.bounds.phase.integral.lower = 0; setup.bounds.phase.integral.upper = 100; setup.guess.phase.time = [0; 0.5; 1]; setup.guess.phase.state = [0; 0.5; 1]; setup.guess.phase.control = [0; 0; 0]; setup.mesh.method = 'hp'; setup.mesh.tolerance = 1e-8; setup.mesh.maxiterations = 10; setup.mesh.colpointsmin = 4; setup.mesh.colpointsmax = 10; output = gpops2(setup); fprintf('J = %.8f\n', output.result.objective); function phaseout = smokeContinuous(input) u = input.phase.control(:,1); phaseout.dynamics = u; % ẋ = u phaseout.integrand = u.^2; % 被积函数 L = u² end function output = smokeEndpoint(input) output.objective = input.phase.integral; % 目标值由积分项汇总 end如果输出 J 偏离 1 超过 1e-6，优先怀疑三个位置：integrand 写成了标量而不是列向量；dynamics 的列数和状态维数不匹配；bounds 上下限里有矛盾。脚本内局部函数要求 MATLAB R2016b 以上版本，并且局部函数要写在脚本末尾。这个用例跑通之后再去翻 examples 目录里的现成例子，才能分清是例子的问题还是环境的问题。
2.4 GPOPS 与 GPOPS-II：同名不同命的两个时代
网上流传的"gpops工具箱"既有第一代 GPOPS，也有彻底重写的 GPOPS-II，两者 API 完全不同。第一代主函数叫 gpops.m，setup 里挂的是 dae、path、objective、events 四个函数字段；GPOPS-II 主函数叫 gpops2.m，函数体系合并成 continuous 和 endpoint 两个。判断解压包里是哪一代，直接看有没有 gpops2.m 即可。第一代没有 hp 网格自适应，配点数要手动试，遇到控制切换容易漏拐点，所以拿到旧版资料时按 GPOPS-II 的注释去改字段必然报错。官方手册会写明对应版本号，不确定就先看手册版本说明页。
3. Bryson-Denham 问题：状态约束下的三段解（带详细注释的完整例子）
3.1 为什么选 Bryson-Denham 做真实用例
Bryson-Denham 是最经典的带状态不等式约束的最优控制测试问题：
min J = 0.5 ∫₀¹ u² dt，ẋ₁ = x₂，ẋ₂ = u，x(0) = (0, 1)，x(1) = (0, -1)，约束 x₁(t) ≤ 1/9。
它的价值在于解析解已知。不加约束时，最优解是 u ≡ -2，x₁ = t - t²，峰值 1/4，超过 1/9，所以状态约束必然激活。加上约束后的最优解分成三段：0 到 2/9 秒 u = -9/2，2/9 到 7/9 秒贴着约束边界走且 u = 0，7/9 到 1 秒 u = -9/2，最优代价 J = 4.5。这个解有三个清晰的控制角点，正好用来观察 hp 网格在哪里加密，也用来验证状态约束的数值处理是否正确，是一个比跑通安装用例高一个档次的 benchmark。
3.2 continuous 函数：dynamics、integrand、path 三个输出
continuous 函数负责给出每个配点上的状态导数、被积函数和路径约束。输入结构里 phase.time、phase.state、phase.control 的行数都是配点数，列数分别对应时间、状态和控制维数。下面这份代码的注释密度就是压缩包里"详细注释"那种写法：
function phaseout = bdContinuous(input) % input.phase.time : N x 1 列向量，配点时间 % input.phase.state : N x 2 矩阵，状态按 bounds 定义的列序排列 % input.phase.control: N x 1 列向量 x1 = input.phase.state(:,1); x2 = input.phase.state(:,2); u = input.phase.control(:,1); phaseout.dynamics = [x2, u]; % 第一列是 dx1/dt = x2 % 第二列是 dx2/dt = u phaseout.integrand = 0.5 * u.^2; % Lagrange 代价的被积函数 % 路径约束的另一种写法： % phaseout.path = x1 - 1/9; % 再在 bounds.phase.path.upper 里设 0 enddynamics 的列数必须与状态维数一致，行数必须与输入行数一致，这是 GPOPS-II 最常见的报错来源。integrand 会沿时间积分，在 endpoint 里变成 input.phase.integral。路径约束如果用 path 输出，就把约束写成 C(x,u,t) 的形式，边界统一放 bounds.phase.path，这样可以混合状态和控制；本例是纯状态约束，直接设状态上界 1/9 更简洁，网格加密后约束在配点上满足得也更干净。
3.3 endpoint 函数：目标值与事件约束的归口
endpoint 函数在每个阶段的首尾各调用一次，用来组装目标函数的 Mayer 项、汇总 Lagrange 项，以及定义事件约束：
function output = bdEndpoint(input) x0 = input.phase.initialstate; % 1 x 2，初始状态行向量 xf = input.phase.finalstate; % 1 x 2，终端状态行向量 output.objective = input.phase.integral; % 目标值 = 积分项 % 若要 Mayer 项，直接叠加，例如： % output.objective = xf(1) + input.phase.integral; % 事件约束放 eventgroup，例如要求终端速度为初始位置的两倍： % output.eventgroup(1).event = xf(2) - 2*x0(1); end本例目标函数只有 Lagrange 项，所以 objective 直接取 integral。注意 initialstate 和 finalstate 是 1 x 2 的行向量，而 continuous 里的 state 是 N x 2 矩阵，两个函数对状态的操作维度习惯不一样，初学最容易在这里把维度写拧。
3.4 完整 setup 与结果核对
把 3.2 和 3.3 的两个函数配上下面的 setup，就是一套完整可复现的例程：
setup = struct(); setup.name = 'bryson_denham'; setup.functions.continuous = @bdContinuous; setup.functions.endpoint = @bdEndpoint; setup.bounds.phase.initialtime.lower = 0; setup.bounds.phase.initialtime.upper = 0; setup.bounds.phase.finaltime.lower = 1; setup.bounds.phase.finaltime.upper = 1; setup.bounds.phase.initialstate.lower = [0, 1]; setup.bounds.phase.initialstate.upper = [0, 1]; setup.bounds.phase.finalstate.lower = [0, -1]; setup.bounds.phase.finalstate.upper = [0, -1]; setup.bounds.phase.state.lower = [-inf, -inf]; setup.bounds.phase.state.upper = [1/9, inf]; % 状态约束挂在状态上界 setup.bounds.phase.control.lower = -20; setup.bounds.phase.control.upper = 20; setup.bounds.phase.integral.lower = 0; setup.bounds.phase.integral.upper = 100; setup.guess.phase.time = [0; 0.5; 1]; setup.guess.phase.state = [0 1; 0.12 0; 0 -1]; % 猜测贴着约束边界 setup.guess.phase.control = [-5; 0; -5]; setup.mesh.tolerance = 1e-8; % 有角点，误差压紧一点 setup.mesh.maxiterations = 20; output = gpops2(setup); sol = output.result.solution.phase; fprintf('J = %.6f\n', output.result.objective); subplot(3,1,1); plot(sol.time, sol.state(:,1), '.-'); hold on; plot([0 1], [1/9 1/9], 'r--'); ylabel('x_1'); grid on; subplot(3,1,2); plot(sol.time, sol.state(:,2), '.-'); ylabel('x_2'); grid on; subplot(3,1,3); plot(sol.time, sol.control, '.-'); ylabel('u'); xlabel('t'); grid on;跑完后用下面的核对表逐项验证：
| 核对项 | 期望值 | 失败时的排查方向 |
|---|---|---|
| 目标函数 J | 4.500000，1e-6 量级偏差 | integrand 的 0.5 系数是否漏掉 |
| 控制角点 | t ≈ 2/9 与 7/9 两处 | 网格是否在角点附近自动加密 |
| x₁ 峰值 | 不超过 1/9 | 状态上界是否写进 bounds |
| 网格迭代次数 | 3~6 次收敛 | maxiterations 是否给够 |
提示：output.result.solution.phase 里的 time、state、control 行数等于最后一次网格迭代的配点数，不要用这个行数反推 guess 的维度。个别早期版本输出字段叫 output.solution 而不是 output.result.solution，以官方手册的输出附录为准。
4. 收敛相关的 6 个参数：网格、导数与 NLP 求解器怎么配
4.1 hp 自适应网格迭代在做什么
每次网格迭代分两步：先在当前网格上解 NLP，再估计解的误差。GPOPS-II 把当前解插值到更密的点上评估动力学残差，残差超标的区间要么增加配点数提高多项式阶数（p 方向），要么把区间一分为二（h 方向）。Bryson-Denham 的两个控制角点处控制不可导，p 方向提阶的收益递减，最终靠 h 方向把网格压到角点附近，这就是状态约束问题网格看起来"一头密一头疏"的原因。理解这个迭代机制，就能解释为什么 mesh.tolerance 和 mesh.maxiterations 比 NLP 内部选项更值得先调：tolerance 决定最终解的误差水平，maxiterations 决定给网格多少次机会去达到这个误差。
4.2 六个核心参数与推荐起步值
下面的表格是多个最优控制问题上沉淀下来的起步值。工程上先求收敛，再谈精度：
| 参数 | 起始值 | 何时调整 |
|---|---|---|
| setup.mesh.tolerance | 1e-6 | 做 benchmark 验证解析解时压到 1e-8~1e-10 |
| setup.mesh.maxiterations | 10 | 状态约束弧线和 bang-bang 控制加到 20~30 |
| setup.mesh.colpointsmin | 4 | 启动就 infeasible 时提到 6~8 |
| setup.mesh.colpointsmax | 10 | 超过 15 收益小且容易振荡 |
| setup.derivatives | 'finite-difference' | 状态维数几十以上时换 'ad' 或稀疏模式 |
| setup.solver | 'ipopt' | 有 SNOPT 授权且问题稀疏性好时切换 |
设完这些参数，判断收敛看两个信号：output.meshhistory 显示最后一次网格迭代的误差估计低于 tolerance；NLP 求解器返回的 exitflag 表示最优。两者缺一不可，网格误差达标但 NLP 没收敛，结果不可信；反过来 NLP 收敛但网格误差没达标，说明解在网格间不一致，需要继续加密。meshhistory 是结构体数组，每一行对应一次网格迭代，直接 disp(output.meshhistory) 看一遍比猜快。
4.3 NLP 求解器与导数的选择
GPOPS 本身不附带 NLP 求解器二进制，IPOPT 是缺省推荐，开源、免费、对非线性约束处理成熟；SNOPT 是商业软件，需要独立授权，优势在冷启动稳定性和大规模稀疏问题的收敛性。IPOPT 的行为通过 setup.nlp.ipoptoptions 下面的字段传入，常见需要动的是输出级别、最大迭代次数和线性求解器：
setup.solver = 'ipopt'; setup.nlp.ipoptoptions.print_level = 5; % 0 静默，12 全量 setup.nlp.ipoptoptions.max_iter = 2000; % 防复杂问题迭代耗尽 setup.nlp.ipoptoptions.linear_solver = 'ma57'; % 有 HSL 库时比 mumps 稳提示：IPOPT 选项的封装字段名在不同版本略有出入，动手前先在官方手册的 NLP 设置一节确认拼写，拼错的字段通常会被静默忽略。
导数方式默认 'finite-difference' 对中小规模问题足够；状态维数几十以上时，有限差分会变成瓶颈，可以尝试设置 setup.derivatives 为 'ad'，但需要对应工具支持，配置复杂度随之上升。对大多数轨迹优化应用，先把网格参数调好，导数保持默认。另一个实用习惯是用 setup.auxdata 传常数（重力、质量、气动系数），在函数里通过 input.auxdata.xxx 读取，不要写全局变量，参数扫描时才能避免改一处全局导致缓存失效。
4.4 三个高频失败模式与现场诊断
第一个是网格第 0 次迭代就 infeasible，九成原因是 guess 与 bounds 矛盾：初始状态上下界写错、guess 越过控制限幅，或者动力学符号给反。诊断办法是把 bounds 全部放松，先解一个无约束版本，再把解作为下一轮 guess。第二个是目标值停在某个精度不再下降，常见于 integrand 写错，Bryson-Denham 里 J 收敛到 9 而不是 4.5，就是 0.5 系数漏了。第三个是路径约束处状态超出允许值，多见于 tolerance 太松，先把 mesh.tolerance 降一个量级，再看是不是 path 约束和 finalstate 边界互相矛盾。这三个问题的共同点是看 meshhistory 里误差曲线的形态就能定位：误差平稳下降是网格问题，误差震荡不降是 NLP 或建模问题。
5. 多段问题的 eventgroup 写法与解的可信度验证
5.1 什么情况下必须拆多段
当动力学、控制界或约束在某个时刻发生本质变化，比如飞行器级间分离、机器人抓取前后负载变化、过程中控制权限改变，单段问题就表达不了。GPOPS-II 的多段机制是把整个时间轴切成 phase，每个 phase 有独立的 bounds、guess 和网格，phase 之间用事件约束（eventgroup）连接。拆段的核心收益是每段可以单独控制网格，段的边界正好落在动力学切换点上，避免用一个高次多项式去逼近本身不光滑的解。
5.2 两段问题的最小写法
一个最小但完整的例子：两段都要满足 ẋ₁ = x₂、ẋ₂ = u，阶段 1 控制限幅 ±10，阶段 2 限幅 ±2，模拟分离后控制能力下降；连接点处位置和速度连续，终端状态固定，总时间固定为 1。continuous 里用循环处理两个 phase：
function phaseout = mpContinuous(input) for p = 1:2 x1 = input.phase(p).state(:,1); x2 = input.phase(p).state(:,2); u = input.phase(p).control(:,1); phaseout(p).dynamics = [x2, u]; phaseout(p).integrand = 0.5 * u.^2; end endendpoint 里把连接条件写进第一个 eventgroup，把终端速度条件写进第二个：
function output = mpEndpoint(input) x1f = input.phase(1).finalstate; % 阶段 1 末端状态 x2i = input.phase(2).initialstate; % 阶段 2 初始状态 output.eventgroup(1).event = x1f - x2i; % 连接点连续 output.eventgroup(2).event = input.phase(2).finalstate(2); % v(T)=0 output.objective = input.phase(1).integral + input.phase(2).integral; end对应的 setup 里 bounds 和 guess 都要带 phase 下标，例如 setup.bounds.phase(1).finaltime.lower、setup.guess.phase(2).state。这里容易踩的坑是：eventgroup(1) 只约束了状态连续，没有约束时间连续，phase(1) 的 finaltime 和 phase(2) 的 initialtime 要分别约束为同一个值，或者在 bounds 里把两个时间固定为相等数值。
5.3 验证技巧：用 Hamiltonian 连续性检验解的质量
多段问题算完别急着用。最优控制理论里，自由连接时间条件下 Hamiltonian 在 phase 边界两侧必须连续。如果工具箱输出了 costate（output.result.solution.phase.costate），就可以做数值验证：在每个 phase 边界内侧取一个配点，按 H = L + λᵀf 计算两侧 H，差值在网格容差量级，说明连接点处给出的是真正的最优切换时间。若版本不返回 costate，退而求其次的做法是把 mesh.tolerance 依次降到 1e-7 和 1e-9，观察目标函数变化小于 1e-8 且网格迭代次数稳定，基本可以认定解已收敛。这个 Hamiltonian 连续性检查对任何多段问题都成立，我每次拆完 phase 都会先跑这一步，比肉眼看轨迹平滑可靠得多。
本文还有配套的精品资源，点击获取
