简介GPOPS4.1是一款基于MATLAB的轨迹优化工具箱面向航天航空、机器人导航等领域的科研人员与工程师专门解决含多决策变量和约束条件的高维度非线性优化、弹道设计与路径规划问题。该工具箱以庞特里亚金最大值原理为核心通过构建哈密顿函数将状态方程与优化目标融合即使在复杂动态系统中也能锁定全局最优控制策略并且无需SNOP等外部求解器独立运行即可完成计算任务。压缩包共291个文件以250个m源码文件为主体同时提供pdf文档、tex源代码、readme/install等说明文件整体仅3.22MB文件类型覆盖源码、文档、安装配置与示例程序查找使用都很方便。资源内置了覆盖多种实际场景的丰富示例用户可快速理解核心语法并迁移到自己的研究中也可以自定义成本函数与约束条件灵活适配特定需求。已有459人学习下载适合需要高效处理轨迹优化问题的中高级MATLAB使用者。1. GPOPS 4.1到底在解决什么问题从“手推协态”到“写出状态方程”做火箭上升段、机械臂或者车辆换道这类轨迹优化问题时最折磨人的往往不是建模而是“初始猜测怎么给”。GPOPS 4.1就是直接面对这个问题的轨迹优化软件它基于高斯伪谱法把连续时间最优控制问题离散成大规模非线性规划NLP问题再交给IPOPT等求解器处理。你只需要把状态方程、目标函数和约束写清楚不用手推协态方程也不用逐步打靶猜初值。对做运动规划与制导控制的人来说这意味着能把精力从算法陷进里抽出来放到“问题本身该怎么描述”上。这套思路适合单段轨迹、多段飞行程序以及带路径约束的复杂工况。2. 建模先行把最优控制问题翻译成GPOPS 4.1能读的Setup结构2.1 先把物理过程改写成状态空间形式状态、控制、约束一个都不能少GPOPS 4.1不会理解你的物理图景它只认识“状态变量、控制变量、微分方程、路径约束、端点约束”这几样东西。所以拿到一个轨迹优化任务第一步是把物理过程写成状态空间形式而不是急着打开MATLAB写代码。我一般习惯在纸上先列三张表。第一张是状态变量表比如高度、速度、质量注意它们必须有明确的单位制而且量纲不要混用第二张是控制变量表像推力、舵偏角、加速度第三张是约束表包括状态和控制的上下界、路径约束、起终点状态约束。这三张表列完问题才算是“翻译”了一半。举个例子模拟项目X简化的一维垂直上升问题里状态高度 h速度 v质量 m控制推力 T动力学dh/dt vdv/dt (T - D(h, v)) / m - gdm/dt -T / (Isp · g0)这里路径约束是 0 ≤ T ≤ T_max以及动压或过载限制端点约束是初始高度、速度、质量和终端高度、终端速度。这里的关键是凡是需要在整个时间区间上成立的限制都属于路径约束必须由连续函数承载凡是只作用在起点或终点的限制都是端点约束由端点函数承载。这个区分直接决定了后面的代码结构。状态个数不是越多越好。GPOPS 4.1把问题最终转成NLP状态多了离散变量规模会成倍增加。能合并成参数的不要设成状态能用常数表示的不要额外引入。比如比冲、参考面积这类固定参数直接写进连续函数里不要变成状态。2.2 Setup函数里的五个输入name、funcs、d、bounds、options缺一不可GPOPS 4.1的入口是一个叫setup的函数它返回一个结构体problem然后主脚本调用gpops(problem)开始求解。这个setup结构里主要有五块内容问题名字、函数句柄、离散参数、边界约束和求解选项。模板长这样function problem setup() problem.name demo_vertical_ascent; % 问题名显示和输出文件用 problem.funcs.continuous continuous; % 连续过程函数 problem.funcs.endpoint endpoint; % 端点函数 problem.bounds.phase.initialtime.lower 0; problem.bounds.phase.initialtime.upper 0; problem.bounds.phase.finaltime.lower 10; problem.bounds.phase.finaltime.upper 30; problem.bounds.phase.state.lower [0; 0; m0 - 500]; problem.bounds.phase.state.upper [1e5; 2000; m0]; problem.bounds.phase.control.lower 0; problem.bounds.phase.control.upper Tmax; problem.bounds.phase.initialstate.lower [0; 0; m0]; problem.bounds.phase.initialstate.upper [0; 0; m0]; problem.bounds.phase.finalstate.lower [h_f; 0; 0]; problem.bounds.phase.finalstate.upper [h_f; 0; m0]; problem.guess.phase.time [0; 15]; problem.guess.phase.state [0, 0, m0; h_f, 0, m0 - 200]; problem.guess.phase.control [0.5 * Tmax; 0.5 * Tmax]; problem.options.nlp.solver ipopt; problem.options.mesh.tolerance 1e-4; problem.options.mesh.maxiterations 10; end注意几个容易忽略的地方initialstate和finalstate是严格定义在t0和tf上的边界状态它们和state的下界不冲突但彼此独立control是标量还是向量要与连续函数里用到的一致guess里的行数要和每个维度对齐。guess本质上是给NLP求解器一个起始点不要求满足全部约束但量级不能离谱。problem.funcs里除了continuous和endpoint还可以挂mesh函数用于自定义网格不过大多数情况下用默认的hp自适应网格就够了。options里的mesh.tolerance是判定网格细化和NLP收敛的基准这个值后面会专门说它是新手最容易设置过头的地方。2.3 连续函数和端点函数的分工路径约束、端点约束和目标的写法continuous函数在GPOPS 4.1的每个网格点上被调用它返回的是动力学导数以及需要在整个轨迹上满足的路径约束。endpoint函数只在积分区间两端被调用它负责计算目标函数里与端点相关的部分并给出端点约束的残差。function phaseout continuous(input) t input.phase.time(:, 1); h input.phase.state(:, 1); v input.phase.state(:, 2); m input.phase.state(:, 3); T input.phase.control(:, 1); hdot v; vdot (T - D(h, v)) / m - g; mdot -T / (Isp * g0); phaseout.dynamics [hdot, vdot, mdot]; phaseout.path T ./ Tmax - overload(h, v); % 路径约束示例 end这里的path约束写成 g(t) ≤ 0 的形式GPOPS 4.1自动处理不等式。endpoint函数则长这样function output endpoint(input) tf input.phase(1).finaltime; hf input.phase(1).finalstate(1); vf input.phase(1).finalstate(2); mf input.phase(1).finalstate(3); output.objective -mf; % 比如燃料最省末端质量最大 output.eventgroup(1).event [hf - h_target; vf - v_target]; endobjective里的量可以是Mayer项、Lagrange积分项或者两者之和。如果目标里有积分项常见做法是增加一个状态变量比如让 x4 的导数等于 0.5·u²然后把objective写成 x4(tf)。这样做比在endpoint里手写积分更稳定也方便后续在网格细化中统一缩放。一个小建议目标函数和约束的所有常量在MATLAB里用全局参数或者闭包传参都可以但千万别在continuous函数里直接写魔法数字。同一套物理量在setup、continuous、endpoint里共用了三处一旦单位不一致GPOPS不会给你报错只会给你一条“看起来平滑但物理上离谱”的轨迹。3. 跑通第一个算例双积分器的控制能量最小轨迹3.1 为什么用双积分器验证GPOPS 4.1它包含你以后要用的所有要素第一个算例不要直接上火箭或者机械臂建议先用一个你有解析解的问题把GPOPS 4.1的调用链跑通。我最常用的是双积分器控制能量最小问题系统是 dx1/dt x2dx2/dt u初始状态 x1(0)0x2(0)0终端状态 x1(1)1x2(1)0目标是最小化 0.5 * ∫u² dt控制约束 -1 ≤ u ≤ 1。这个问题的好处是有解析解可以对照GPOPS解出来偏差应该在1e-4以内状态少NLP规模小整个求解在十几秒内结束非常适合调试环境。它同时包含了需要写在continuous里的动力学、需要写在endpoint里的端点约束、需要作为状态引入的目标积分以及控制量边界。可以说以后你写复杂问题会碰到的结构这个例子里全都有。为了把积分目标变成Mayer形式我引入第三个状态 x3使 dx3/dt 0.5 * u²最终目标就是 x3(1) 最小。这样endpoint里不需要手动做数值积分一切都交给配点法和NLP求解器处理。3.2 完整可运行脚本main、setup、continuous、endpoint主脚本部分% main_run_double_integrator.m clc; clear; close all; addpath(genpath(gpops4.1)); % 按你的实际目录修改 problem setup(); solution gpops(problem); % 输出关键结果 t solution.phase(1).time; x1 solution.phase(1).state(:, 1); x2 solution.phase(1).state(:, 2); x3 solution.phase(1).state(:, 3); u solution.phase(1).control(:, 1); fprintf(目标函数 J %.6f\n, x3(end));setup函数function problem setup() problem.name double_integrator_energy; problem.funcs.continuous continuous; problem.funcs.endpoint endpoint; problem.bounds.phase.initialtime.lower 0; problem.bounds.phase.initialtime.upper 0; problem.bounds.phase.finaltime.lower 1; problem.bounds.phase.finaltime.upper 1; problem.bounds.phase.state.lower [-10; -10; 0]; problem.bounds.phase.state.upper [ 10; 10; 10]; problem.bounds.phase.control.lower -1; problem.bounds.phase.control.upper 1; problem.bounds.phase.initialstate.lower [0; 0; 0]; problem.bounds.phase.initialstate.upper [0; 0; 0]; problem.bounds.phase.finalstate.lower [1; 0; 0]; problem.bounds.phase.finalstate.upper [1; 0; 10]; problem.guess.phase.time [0; 0.5; 1]; problem.guess.phase.state [0, 0, 0; 0.5, 1, 0.1; 1, 0, 0.2]; problem.guess.phase.control [0; 0; 0]; problem.options.nlp.solver ipopt; problem.options.nlp.ipopt.tol 1e-7; problem.options.nlp.ipopt.max_iter 500; problem.options.mesh.tolerance 1e-4; problem.options.mesh.maxiterations 10; endcontinuous和endpoint函数function phaseout continuous(input) x1 input.phase.state(:, 1); x2 input.phase.state(:, 2); u input.phase.control(:, 1); dx1 x2; dx2 u; dx3 0.5 * u.^2; phaseout.dynamics [dx1, dx2, dx3]; end function output endpoint(input) output.objective input.phase(1).finalstate(3); output.eventgroup(1).event [ input.phase(1).finalstate(1) - 1; input.phase(1).finalstate(2) - 0 ]; end这里有个值得说的细节控制约束写在bounds里端点约束写在event里而状态的下界在代码里设成 -10 和 10只是给NLP求解器一个安全范围不要让状态乱跑。真实轨迹里x1不会超过0.5x2也不会超过1但你给的范围太紧而初值又不好时求解器会很不舒服。代码里没有显式设置初始网格的配点数GPOPS 4.1默认会用一个较粗的网格开始。第一次跑建议用默认值不要一上来就搞很多配点。3.3 第一次运行的三个注意点求解器选项、初值猜测和缩放第一个注意点是求解器选项。IPOPT是GPOPS 4.1最常用的NLP求解器它稳定、开源而且对大规模稀疏问题支持很好。tol建议先设1e-6到1e-7不需要追求1e-10。GPOPS 4.1的收敛是靠网格细化来保证的NLP求解器本身咬得太死只会增加算力消耗不会明显提高轨迹精度。如果你的机器内存小可以把linear_solver设为mumps然后关掉IPOPT的输出problem.options.nlp.ipopt.linear_solver mumps; problem.options.nlp.ipopt.print_level 5;第二个注意点是初值猜测。GPOPS 4.1对初值比传统打靶法宽容得多但也不是完全免疫。guess里的时间序列、状态序列和控制序列必须行数相同且要覆盖整个时间区间。最省事的办法是用线性插值起终点状态明确就填起终点中间随便插几个点完全没概念的控制量填一个处于约束中位的常数比如0而不是边界值。我在给别人排查初值问题时见过最多的翻车现场是控制猜测填了最大值然后NLP第一步就计算出巨大的斜率直接导致不可行。记住初值的第一使命是让NLP求解器进入可行域不是让你猜准最优解。第三个注意点是缩放。GPOPS 4.1内部会做一部分变量缩放但你的状态量如果跨越五个数量级比如高度到1e5、加速度只有10NLP的雅可比矩阵条件数会很难看。常见做法是在写continuous之前就把物理单位改造成归一化或半归一化量。等模型跑通了再把单位换算放回后处理里。4. 结果解读与参数调优检查可行性、网格细化与可视化4.1 从solution里取出轨迹数据时间和状态、控制怎么对应GPOPS 4.1返回的solution结构体phase(1)里的time、state、control是列向量state的每一行对应time那一行的时刻。注意不是等间隔采样靠近边界和拐点的地方点更密这是伪谱法配点的正常分布。t solution.phase(1).time; x solution.phase(1).state; % 每一列是一个状态变量 u solution.phase(1).control; % 每一列是一个控制变量 figure; plot(t, x(:, 1), LineWidth, 1.5); hold on; plot(t, x(:, 2), LineWidth, 1.5); plot(t, u, k--, LineWidth, 1.2); grid on; legend(x1, x2, u); xlabel(Time (s));GPOPS 4.1也自带绘图函数可以直接把状态和控制画出来。不过我还是习惯自己去plot原因是自带的图一般把多段phase画在一张图上相位之间如果没有明显区分很容易看错段。自己画可以按需调整线型、标注关键点。取数据时有一个高频失误点是控制量取出来之后直接拿去给高保真模型用。GPOPS 4.1的控制是在配点处满足约束的你如果要做一个时间步进的仿真必须对控制做插值或者把轨迹重新细化为等间隔时间序列不能直接把数组当成连续控制律。4.2 判断这次求解可不可信退出标志与约束违反量很多人在跑完GPOPS后只看“Optimization terminated”这一句话就认为万事大吉这是不对的。IPOPT输出的“求解成功”只代表NLP在这组离散点上找到了局部最优不代表原最优控制问题被解出来了。正确做法是看三样东西NLP求解器的exitflag是不是0或者IPOPT的status是不是“Solve_Succeeded”端点约束的event残差比如 finalstate(1) - 1 是不是小于1e-5路径约束的最大违反量在continuous里定义的所有path项是否都满足符号要求。GPOPS 4.1的solution里通常会有约束违反量的记录没有的话就用输出解自己回代验算。路径约束别只看配点上的值因为真正起作用的是配点和它们之间的插值多项式的值。如果配点上满足、点间不满足通常是网格点数不够或者约束本身有突变。表格做一个常见现象的快速定位现象可能原因检查方向exitflag0但轨迹明显不合理初值差落进局部解改guess或加网格点数event残差大endpoint里等式写错检查残差表达式或单位路径约束在配点满足但曲线穿出网格太粗调低mesh.tolerance目标函数比解析解差不少网格细化没跑完看mesh iteration次数4.3 网格细化容差怎么定精度与运行时间的权衡表mesh.tolerance是GPOPS 4.1里最影响求解时间的参数。它决定伪谱法在哪些区间要加密配点容差越小网格越细NLP变量越多单次IPOPT求解越慢。新手最容易犯的错误是一上来设1e-8然后把原本半分钟能解决的问题变成半小时而且经常因为网格过细导致IPOPT在内存或迭代上翻车。我个人的经验是先设1e-3跑通看轨迹形状合理、约束满足、目标函数与解析解的量级一致再改成1e-4或1e-5做最终求解。对于实际工程场景1e-4到1e-3往往已经够用。网格过多的另一个副作用是后处理时你拿到的配点可能多达几百个用这些点做控制插值时反而更容易出现插值振荡。你可以在setup里单独控制网格结构problem.options.mesh.tolerance 1e-4; problem.options.mesh.maxiterations 12; problem.options.mesh.nodesPerInterval.min 4; problem.options.mesh.nodesPerInterval.max 30;网格自适应是会迭代的maxiterations设太小可能导致网格还没细化到位就被强制结束。这个时候看solution.phase(...).mesh里输出的配点数量如果连续两次迭代都还在增加说明精度没到阈值需要给它更多迭代次数。5. GPOPS 4.1常见问题与避坑求解失败、初值敏感与约束错位5.1 现象IPOPT报“Restoration failed”或者迭代直接NaN这个报错是最劝退的。现象是IPOPT在前几次迭代里一步都走不动最后输出Restoration failed或者目标函数变成NaN然后崩掉。原因绝大多数出在continuous函数写得不光滑要么用了abs、max、min、sign这类不可导运算要么在动力学里写了分段的查表逻辑。解决方法是把不可导运算改造成可导形式。比如限幅约束不要写成 T min(10, T_cmd)而是把 T 直接当作控制变量让求解器自己去满足边界阻力系数曲线D(h,v)如果是查表函数先用多项式或样条拟合出光滑表达式再写进continuous里。% 错误写法直接查表 D interp1(table_h, table_D, h, linear, extrap); % 推荐写法事先用多项式拟合 D polyval(p_coeff, h) polyval(q_coeff, v);如果你确认函数本身是光滑的那再检查一下初值和边界范围。状态范围给太窄求解器在边界上反复碰壁也容易出现Restoration failed。5.2 现象初值很差导致结果“看起来对但物理上不对”有时候求解成功了画出来的轨迹也很平滑但你知道从物理上解释不通比如燃料全程不变、推力一直顶在边界上却还满足终端条件。这通常是求解器落进了一个数学上可行的局部解而这个解对应的物理过程是不成立的。这类问题靠调NLP选项解决不了要从初值入手。GPOPS 4.1对全局解没有保证伪谱法把问题转成NLP之后局部性由NLP求解器继承。我的习惯是先做一个“简化版”把要求降低比如去掉过载约束或者放宽终端窗口让目标函数主导一切得到一个粗略解然后用这个粗略解作为完整问题的guess。这样比拍脑袋猜状态序列可靠得多。5.3 现象端点条件不满足而状态曲线穿出边界求解器明明报告成功但你检查output.eventgroup(1).event残差发现并不为零或者状态曲线在起点附近很明显超出了initialstate设定的值。这类问题的根源几乎都是把端点约束写在了错误的位置。有一些写法状态在配点处的边界是满足的但拉格朗日插值多项式在两个配点之间的振荡让曲线穿出上下界这是伪谱法的固有现象。处理办法是加密网格或者在动力学里增加过渡缓冲。但更常见的错误是端点约束写在了continuous的path里让求解器误以为“所有时刻都要满足端点值”结果端点没卡住反而把路径逼弯了。确认一下initialstate和finalstate写在bounds里事件残差写在endpoint里路径约束只保留真正“全程”的限制。5.4 现象控制量剧烈震荡末端出现抖动控制曲线在末端高频小幅震荡是伪谱法里典型的现象尤其在时间自由、终端状态约束多的算例中出现概率很高。原因有两类一是网格在这个区间过粗高频控制分量不被分辨二是终端状态约束迫使控制做剧烈调整数值上容易出现吉布斯式抖动。处理顺序也固定先把mesh.tolerance往下调一档看控件震荡是否减弱如果没减弱检查终端的约束是不是给得过强比如终端速度必须严格为0还同时要求终端加速度为0这会使系统的自由度不足。一个有效手段是给目标函数增加一个控制变化率的惩罚项也就是把u对时间的导数写进目标里让系统自己平衡“省控制”和“平滑控制”。5.5 现象MATLAB报错找不到函数或求解器最常见的是在另一台电脑上打开工程时gpops文件夹没有加入MATLAB路径。这个好解决用addpath(genpath(...))一次性加进来即可。第二个常见问题是IPOPT的mex文件是64位编译的而你的MATLAB是32位或换了架构导致调用求解器时直接报“Invalid MEX-file”。我一般在新环境里做的第一件事是用GPOPS 4.1自带的测试脚本跑通自带算例确认求解器可用再跑自己的模型。还有一个隐蔽坑MATLAB当前工作目录下如果有某个脚本和GPOPS函数重名比如你也写了一个continuous.mMATLAB会优先执行当前目录那个导致行为完全混乱。建议把自定义函数名都改成带前缀的名字比如demo_continuous而不是裸的continuous。6. 进阶用法多phase拼接、解导出与GPOPS 4.1的使用习惯6.1 多段任务的phase连接连接点状态连续的条件实际工程任务很少是一个连续段走到底的比如火箭上升段结束后要滑行机械臂抓取前要切换运动模式。GPOPS 4.1支持把整个任务拆成多个phase每个phase有自己的时间、状态和控制phase之间用连接条件互相约束。problem.bounds.phase(2).initialtime.lower problem.bounds.phase(1).finaltime.lower; problem.bounds.phase(2).initialtime.upper problem.bounds.phase(1).finaltime.upper; problem.bounds.phase(2).initialstate.lower problem.bounds.phase(1).finalstate.lower; problem.bounds.phase(2).initialstate.upper problem.bounds.phase(1).finalstate.upper;这样第二个phase的起点被绑定到第一个phase的终点。连接条件还可以只约束部分状态连续比如位置连续但速度允许跳变用于撞击或抛投类任务。多phase的好处是每个phase可以独立设网格精度和计算量的分配更灵活。6.2 把GPOPS 4.1的结果变成其它工程工具的初始轨迹把配点解导出成等间隔时间表是GPOPS结果最常被用到的落地方式。我习惯在后处理里做两步第一步把时间和控制量用样条插值第二步按固定步长重新采样保存成结构体或表格。t_fine linspace(t(1), t(end), 2000); u_fine interp1(t, u, t_fine, spline); state_fine interp1(t, x, t_fine, spline);插值后一定要重新检查控制约束边界。样条插值很容易在穿越上下界的地方冲出一点点如果你之后要拿这个表去驱动仿真模型别忘了加饱和限幅或者改用pchip插值避免振荡。有一点我想特别强调GPOPS 4.1的价值是帮你快速得到一个可行的、接近最优的基准轨迹它不代替高保真仿真。我的收尾习惯是每跑完一组算例就在工程笔记里记下三件事问题规模、网格容差、最终约束违反量。下次调参数时直接翻这个表省掉大量重复试错。这个习惯让我少踩了很多坑也希望帮到你。本文还有配套的精品资源点击获取