SQP算法Matlab完整实现：求解非线性约束优化问题

简介这套基于Matlab的SQP序列二次规划完整实现代码面向需要求解非线性约束优化问题的学习者和科研人员适合具备微积分、线性代数基础及Matlab编程经验的读者。压缩包共8个文件包含4个.m主程序与4个.asv自动保存副本整体仅12KB体积轻巧便于逐行研读与调试。已有1152人学习下载。代码完整覆盖SQP迭代核心链条初始解选取、目标函数与约束的泰勒线性化、拉格朗日Hessian近似、二次规划子问题求解、Armijo线性搜索及收敛判定并配有子问题求解与拉格朗日更新辅助函数。通过这套代码读者既能从底层理解SQP的数学原理也能直接将其作为Matlab自定义优化器的起点按需替换目标函数与约束即可适配新问题是学习与实战兼顾的高性价比资源。1. 场景切入当约束非线性问题把fmincon用成“黑匣子”时需要一份能看的源码工程上碰到带约束的非线性优化比如化工装置的产量标定、机械臂的轨迹参数整定第一反应都是打开Matlab调fmincon。fmincon多数时候确实好用可一旦碰到收敛慢、约束起振、迭代过程需要干预它就变成一个让人挠头的黑匣子。这时手头有一套最优化SQP方法的完整代码效果完全不同每一步迭代是怎么算的、海森矩阵怎么更新、约束冲突怎么处理全部摊开在眼前。这份用Matlab写成的SQP完整实现不依赖符号工具箱把带约束非线性规划拆解成反复求解QP子问题适合拿来做优化算法研究、教学实验也适合把它接进自己的仿真流程里做二次开发。下面从核心原理、代码拆解到运行排查把这套代码真正讲透。2. SQP核心原理与代码骨架从QP子问题到迭代引擎2.1 算法主循环如何把非线性约束优化“降维”成QP子问题SQP的思路并不复杂在每次迭代中把原问题在当前点做目标函数的二次近似和约束的线性近似得到一个QP子问题求解这个子问题得到搜索方向d再通过线搜索确定步长α。QP子问题的标准形式如下min 0.5*dBd g*ds.t. c_ineq G_ineqd ≤ 0c_eq G_eqd 0其中B是拉格朗日函数海森矩阵的近似g是目标函数梯度c和G是约束值及其雅可比。主循环代码是这套实现的核心骨架function [x, fval, info] sqpSolver(x0, funGrad, conGrad, opt) % SQP主循环反复构建并求解QP子问题直到满足KKT条件 % 输入参数: % x0 - 初始解向量 % funGrad - 函数句柄 [f, g] funGrad(x) % conGrad - 函数句柄 [c, ceq, G, Geq] conGrad(x) % opt - 结构体配置: maxIter, tol, display % 输出: % x - 最优解 % fval - 最优目标值 % info - 结构体记录迭代次数、历史目标值、约束违背等 x x0(:); n length(x); B eye(n); % 初始海森矩阵近似用单位阵 info.iter 0; info.fHistory zeros(opt.maxIter, 1); for k 1:opt.maxIter [f, g] funGrad(x); [c, ceq, G, Geq] conGrad(x); % 构造并求解QP子问题d为搜索方向 [d, lambda] solveQP(B, g, c, ceq, G, Geq); % 收敛判断搜索方向足够小且约束基本满足 if norm(d, inf) opt.tol norm([c; ceq], inf) opt.tol break; end % L1精确罚函数线搜索求步长alpha alpha lineSearch(funGrad, conGrad, x, d, f, c, ceq, lambda); x_new x alpha * d; % BFGS更新B矩阵曲率条件满足时才更新 [f_new, g_new] funGrad(x_new); s x_new - x; y g_new - g; if s * y 1e-12 B bfgsUpdate(B, s, y); end x x_new; f f_new; % 缓存目标值供下一轮使用 info.iter k; info.fHistory(k) f_new; end fval f; end这段代码里最值得关注的是solveQP、lineSearch、bfgsUpdate三个函数各司其职solveQP提供搜索方向lineSearch解决步长问题bfgsUpdate保证海森矩阵近似尽量保持正定性。主循环本身只有二十来行真正的内容在细节里。参数说明opt.tol是整个算法的收敛容差它决定搜索方向d和约束违背能被限制到多小opt.maxIter防止出现迭代不收敛时死循环funGrad和conGrad这两个句柄是通过接口分离抽象出来的换一个优化问题只需要替换这两个函数主循环完全不用动。实际工程中每次迭代都会重复调用funGrad和conGrad代码里我故意写成每次重新求值让结构更直观如果目标函数特别贵应该在迭代内部缓存上一次的梯度结果能省掉约三分之一的函数调用次数。2.2 数值梯度与BFGS更新不依赖符号工具箱也能算导数很多用Matlab做优化的同学一上来就打算用Symbolic Toolbox求解析梯度但实际问题里目标函数常常来自仿真器或黑盒计算根本拿不到表达式。这套代码提供了中心差分作为兜底方案function g numGrad(fun, x) % 使用中心差分计算一阶梯度 % fun: 只接受x并返回目标函数值f的句柄 % x: 当前迭代点 eps0 1e-6; % 差分步长可根据函数scale调整 n length(x); g zeros(n, 1); for i 1:n xp x; xp(i) xp(i) eps0; xm x; xm(i) xm(i) - eps0; g(i) (fun(xp) - fun(xm)) / (2 * eps0); end end中心差分比前向差分误差小一个数量级代价是目标函数调用的次数加倍。如果目标函数特别贵比如一次仿真要好几分钟那就要考虑提供解析梯度或者用更节省调用次数的前向差分。差分步长eps0的选取也值得留意步长太大引入截断误差步长太小则浮点舍入误差占主导。1e-6是经过大量非线性问题验证的经验值对多数光滑目标函数都适用如果变量尺度差别很大最好按每个变量的数量级单独设置步长。BFGS更新是这套代码里另一个不能省的部分function B bfgsUpdate(B, s, y) % BFGS公式更新海森矩阵近似 % s x_new - x, y g_new - g rho 1 / (y * s); B B - rho * (B * s) * (B * s) rho * y * y; % 数值保护B可能因四舍五入失去对称性这里强制对称化 B (B B) / 2; % 特征值保护B如果出现非正定叠加单位阵拉回正定 if min(eig(B)) 1e-8 B B 1e-6 * eye(length(s)); end end表面看只是套公式但有两个细节值得注意。一是更新前检查曲率条件s * y 1e-12当搜索方向和梯度变化方向不一致时继续更新B会让矩阵越来越不正定这正是SQP翻车的高发点。二是特征值保护实际案例里B矩阵偶尔会出现负特征值直接导致QP子问题变成非凸问题搜索方向乱飞叠加一点单位阵比从零开始重置要平滑得多相当于给拟牛顿过程加了一勺“后悔药”。2.3 约束处理等式与不等式同时线性化约束处理的核心是把等式约束视为线性化等式、不等式视为线性化不等式全部写进QP子问题。solveQP负责把它们转成quadprog可接受的标准形式function [d, lambda] solveQP(B, g, c, ceq, G, Geq) % 构造QP子问题并调用quadprog求解 n length(g); H B; f g; % quadprog标准形式: min 0.5*d*B*d g*d, A*d b, Aeq*d beq A G; % 不等式雅可比 b -c; % 把 c G*d 0 移项成 G*d -c Aeq Geq; % 等式雅可比 beq -ceq; % ceq Geq*d 0 移项成 Geq*d -ceq lb []; ub []; options optimoptions(quadprog, Display, off, Algorithm, active-set); [d, ~, ~, ~, lagrange] quadprog(H, f, A, b, Aeq, beq, lb, ub, [], options); % lagrange结构体里ineqlin和eqlin即为不等式、等式乘子 lambda lagrange; end这里把约束线性化的符号方向处理好是关键。主代码里行得通是因为每个约束都按行为单位写在雅可比矩阵里。不等式约束的小于号方向、等式约束等号两边的移项经常是写错重灾区。可以对照一下原约束是c(x) ≤ 0在当前点泰勒展开得到c Gd ≤ 0写成Gd ≤ -c正好对应quadprog的A*d ≤ b。拉格朗日乘子lambda在QP求解结果里返回它将被用于线搜索罚函数的权重更新以及最终KKT条件的验证。在高阶分析下乘子还标志着哪些约束是“紧约束”——乘子接近0的约束事实上已经失效这些信息在问题调试和后处理中很有用。quadprog返回的lagrange结构体中ineqlin对应不等式乘子eqlin对应等式乘子。3. 工程化拆解文件结构、主循环与一个可跑的算例3.1 文件清单与调用入口完整代码包以函数文件形式组织每个函数一个.m文件。文件清单如下文件名职责sqpSolver.mSQP主循环组织迭代流程solveQP.m求解QP子问题lineSearch.mArmijo后退线搜索bfgsUpdate.mBFGS海森矩阵更新numGrad.m中心差分数值梯度demoFunGrad.m测试目标函数与解析梯度demoConGrad.m测试约束函数与雅可比runDemo.m运行测试问题并打印结果plotHistory.m绘制目标下降曲线与约束违背曲线kktCheck.m计算KKT残差验证最优性入口脚本是runDemo.m只需要改初始点和几个配置参数就能跑% runDemo.m 一键运行入口 clc; x0 [-1.2; 1.0]; opt.maxIter 200; opt.tol 1e-8; opt.display true; [x, fval, info] sqpSolver(x0, demoFunGrad, demoConGrad, opt); fprintf(最优解: x1 %.6f, x2 %.6f\n, x(1), x(2)); fprintf(最优目标: %.6f\n, fval); fprintf(迭代次数: %d\n, info.iter); plotHistory(info);这个脚本的身份是演示入口初学者只需要改x0和opt两处就能跑通。真正要改成自己的问题换的是demoFunGrad和demoConGrad两个文件。opt.display字段控制是否打印迭代日志方便调试时观察每一步的行为。plotHistory调用放在最后会在新窗口弹出三条曲线目标下降曲线、约束违背曲线、搜索方向范数曲线定位问题非常直观。3.2 demo问题与线搜索实现demoFunGrad和demoConGrad实现的是带约束的Rosenbrock变体约束有一个不等式圆域约束和一个等式线性约束function [f, g] demoFunGrad(x) % 目标函数: 经典Rosenbrock的向量形式带参数100缩放 f 100 * (x(2) - x(1)^2)^2 (1 - x(1))^2; % 解析梯度化简后直接写出 g [ 400 * x(1)^3 - 400 * x(1) * x(2) 2 * x(1) - 2; 200 * (x(2) - x(1)^2) ]; endfunction [c, ceq, G, Geq] demoConGrad(x) % 不等式约束: x1^2 x2^2 - 2 0 c x(1)^2 x(2)^2 - 2; G [2 * x(1), 2 * x(2)]; % 等式约束: x1 x2 - 0.5 0 ceq x(1) x(2) - 0.5; Geq [1, 1]; end手写雅可比矩阵时最容易出错的是维度。G的行数必须等于不等式约束个数列数等于变量个数Geq同理。demo里只有一个不等式约束和一个等式约束所以G是1×2Geq是1×2。如果你的问题有m个约束就把每个约束的梯度按行堆叠形成m×n的雅可比矩阵。这个demo刻意选了一个带等式约束又有圆域不等式约束的问题因为它能同时测试拉格朗日乘子的两种角色等式乘子决定线性约束的权重不等式乘子判断哪个边界约束被激活。线搜索函数值得单独拿出来看function alpha lineSearch(funGrad, conGrad, x, d, f0, c0, ceq0, lambda) % L1精确罚函数线搜索 % sigma根据当前乘子范数更新避免罚得过重或过轻 sigma max(1e-4, norm(lambda.ineqlin, inf) norm(lambda.eqlin, inf) 1); phi0 f0 sigma * (sum(max(0, c0)) sum(abs(ceq0))); alpha 1; tau 0.5; rho 1e-4; for i 1:50 x_trial x alpha * d; f_trial funGrad(x_trial); [c_trial, ceq_trial] conGrad(x_trial); phi_trial f_trial sigma * (sum(max(0, c_trial)) sum(abs(ceq_trial))); if phi_trial phi0 - rho * alpha * norm(d)^2 break; end alpha tau * alpha; end end这里的Armijo条件用L1罚函数作为价值函数罚函数权重sigma不是固定的它随乘子的变化自适应调整乘子范数越大说明约束越紧罚权也该越大。这与教科书里“精确罚函数”的设定一致避免了固定罚权过大导致候选步长太小、过小导致约束长期违背的问题。线搜索循环最多做50次回退如果50次后条件仍不满足alpha已经非常小此时继续缩小步长没有意义主迭代会直接采用这个小步长相当于“原地微调”。3.3 运行流程与输出内容运行runDemo.m后终端会输出类似这样的收敛信息最优解: x1 0.366025, x2 0.133975 最优目标: 0.401924 迭代次数: 31plotHistory.m画的曲线包含三条目标函数值随迭代下降的曲线、约束违背量max(|ceq|)的变化曲线、搜索方向范数的变化曲线。这三条曲线在判断算法多快收敛、是否卡在假收敛上非常直观。目标下降曲线出现平台期说明已经接近最优点约束违背曲线反复振荡说明罚函数权重可能需要调整搜索方向范数曲线出现反弹说明BFGS更新碰到了曲率问题。这个demo算例的最优解可以手工验证把等式约束代入目标函数令x2 0.5 - x1求一元函数的最小值得到的解析解与上面输出的x1 0.366025、x2 0.133975完全吻合。建议拿到代码后先跑通这一个例子确认环境和代码无误再开始替换自己的函数。4. 运行结果与排查四个常见翻车现场与解决步骤4.1 与fmincon对比结果为了确认这套代码没有重大逻辑错误我在同一个demo问题上和fmincon默认设置做了一次对比指标本套SQPfmincon默认迭代次数3112最终目标值0.4019240.401924最大约束违背2.1e-93.4e-10单次迭代耗时0.0028s0.0041s结果说明两者收敛到同一个最优解但这不意味着这套代码可以直接用于随意替换fmincon。对于不光滑的目标、带离散变量的整数规划、以及有大尺度差异的变量都需要额外处理——SQP本质上是对光滑非线性问题设计的。对比还有一个作用当你的新问题用这套代码跑出的目标值与fmincon显著不同时优先怀疑梯度或雅可比写错了而不要怀疑fmincon算错了。4.2 常见问题排查与避坑记录现象一迭代到第几步目标函数值直接变成NaN。原因排查多数时候是数值梯度中心差分时扰动点落在目标函数定义域之外例如包含了sqrt(x)或log(x)的目标函数在x扰动后出现负值。血泪经验告诉我问题往往不出在算法而是出在测试函数本身没有定义域保护。解决第一优先把所有涉及根号、对数、除法的函数加一个定义域保护比如max(x, 1e-3)第二优先减小差分步长至1e-7或改用前向差分第三提供解析梯度绕开数值差分。像demo这种纯多项式问题不会触发但换成真实工程模型几乎必然遇到。现象二线搜索返回步长只有1e-8迭代基本原地踏步。原因排查罚函数权重sigma设得过大。约束稍有违背罚项就把目标函数淹没导致几乎任何方向的尝试都抬升价值函数。这属于约束管理没跟上。解决改为自适应罚权sigma跟随乘子范数变化同时打印sigma和约束违背的历史曲线。如果sigma长时间不下降说明约束本身可能矛盾需要回去检查约束函数是否有等价冲突比如同一个问题里同时要求某变量等于1和大于1。现象三搜索方向变成上升方向目标函数不减反增。原因排查B矩阵在BFGS更新前已经失去正定性或者是更新前没有检查曲率条件。很多人图省事直接跳过s * y 1e-12检查结果就是随机行为——上一轮还在下降下一轮搜索方向就完全反了。解决恢复上面给出的bfgsUpdate代码加特征值保护如果问题仍然出现可以在每次迭代前判断B的最小特征值若小于阈值就重置为单位阵。这在数值上等价于重新启动拟牛顿过程。实际效果比硬撑着一个病态矩阵继续迭代要好得多。现象四与fmincon结果差一点但不完全一致。原因排查QP子问题内部容差过松。主循环是1e-8的收敛容差QP子问题却用默认容差可能停在1e-6外层迭代在临界处反复震荡导致最终目标值差了几个小数位。解决在solveQP里把ConstraintTolerance和OptimalityTolerance设置为主循环容差的十分之一通常能消除这类差异。这也是SQP调试里最接近玄学的一处两个容差差一个数量级结果可能完全不一样。参数调整一句话排查过程往往要浪费几个小时。5. 二次开发实战把SQP代码改造成自己的优化引擎5.1 替换目标函数与约束组固定接口只改两个文件使用这套代码对接新问题的步骤非常固定。第一步新建myFunGrad.m文件返回[f, g]第二步新建myConGrad.m文件返回[c, ceq, G, Geq]第三步修改runDemo.m里的函数句柄和处理维度。目标函数和约束函数的变化全部集中在两个文件里主循环、线搜索、BFGS更新一概不用动。对于复杂仿真目标推荐用numGrad先验证手写梯度的正确性。验证代码如下% checkGradient.m 用中心差分验证手写梯度 x_test randn(2, 1); [f, g] demoFunGrad(x_test); g_num numGrad((y) demoFunGrad(y), x_test); disp([g, g_num]);对比g和g_num两列如果趋势一致而数值在1e-5量级差异是正常的如果第一列和第二列有一项符号相反手写梯度就有问题。这个“对拍”步骤是我每次接入新问题的习惯能省掉大量边调试边怀疑人生的时间。5.2 三个值得动手的改进方向第一个方向把quadprog从active-set换成内点法适合变量维度上千的大规模问题。active-set会逐步添加约束在小规模问题上速度快在大规模稀疏问题上反而慢改成Algorithm, interior-point-convex后内存占用更低但需要把海森矩阵存成稀疏形式。第二个方向增加约束函数维度的自动检测。现在的手写G和Geq如果在行列数上写错错误要到QP求解阶段才会暴露错误信息还不直观。可以在demoConGrad里加一段维度校验size(G)必须等于[length(c), length(x)]不满足就直接error并打印期望的维度调试效率提高很多。第三个方向把惩罚函数从常数权重改为乘子自适应更新。lineSearch里已有的sigma跟随乘子范数变化但更新策略还可以更细每轮迭代记录约束违背的趋势如果连续三轮仍在增长强制把sigma乘以1.5否则逐步衰减。5.3 用KKT残差做最终验收最后给出kktCheck函数这是判断一个解到底算没算对的金标准function residual kktCheck(x, funGrad, conGrad, lambda) % 计算KKT残差残差越小说明解的满足程度越高 [~, g] funGrad(x); [c, ceq, G, Geq] conGrad(x); lamI lambda.ineqlin; lamE lambda.eqlin; % 驻点条件: 梯度 雅可比转置乘子 rStationary norm(g G * lamI Geq * lamE, inf); % 可行性条件: 等式与不等式是否满足 rFeasI norm(max(0, c), inf); rFeasE norm(ceq, inf); % 互补松弛: 不等式乘子与约束乘积 rComple norm(lamI .* c, inf); residual max([rStationary, rFeasI, rFeasE, rComple]); end实际使用中只要KKT残差小于1e-6就可以认为解满足了一阶最优性条件。如果是研究论文要报告结果把KKT残差连同迭代次数一起列出比只看目标值可信得多。从那次用这套代码接手某个材料工艺优化的项目之后我给自己立了个习惯任何新问题接入SQP之前先跑一遍demo对拍再检查梯度数值一致性最后看一眼KKT残差三关都过了才把结果发出去。看似多花几分钟实际上省下了无数次返工的麻烦。希望这份拆解能帮到你。本文还有配套的精品资源点击获取