简介一份面向MATLAB学习者的二次规划源码包适合正在研究最优化理论、工程优化、经济建模或信号处理等领域的算法学习者。代码演示了如何用quadprog和自定义算法求解带线性约束的二次规划问题涵盖Hessian矩阵、梯度向量、不等式约束与等式约束的设置方式并提供了增广拉格朗日乘子法、有效集法和路径追踪等多种实现视角方便对照教材理解算法迭代细节。包内共4个文件均为.m脚本压缩包仅2KB轻量易读便于逐行运行和调试。已有293人学习下载。通过研读源码读者既能掌握quadprog的标准调用流程也能学习自定义优化选项、迭代限制和精度控制等高级用法同时可参照脚本修改Q、c、A、b、E、d等矩阵快速套用到投资组合优化、成本最小化等实际场景中是课程设计或论文实验的可靠起点。 最近在帮一个做金融风控的朋友写资产配置模型目标函数带了一堆约束条件本质上就是个典型的二次规划问题。MATLAB里直接把quadprog这个函数调用起来调参、跑通、验证结果整个过程踩了不少坑也积累了一些经验。这篇就围绕“MATLAB 源代码 二次规划”这个主题把从建模到代码落地、再到排查问题的方法完整梳理一遍。适用人群包括正在做毕业设计的学生、刚接触优化算法的工程师以及需要在项目中套用二次规划模型的开发者。1. 先搞明白二次规划到底在优化什么1.1 数学表达式的每个符号都别放过二次规划Quadratic Programming简称QP的标准形式长这样min 0.5 * x * H * x f * x 约束条件 A * x b Aeq * x beq lb x ub这个写法里x是决策变量向量H是二次项的Hessian矩阵f是线性项的系数向量A和b定义不等式约束Aeq和beq定义等式约束lb和ub是决策变量的上下界。为什么要写成0.5 * x * H * x而不是x * H * x因为这样求导之后梯度正好是H * x f配合最优性条件H * x f 0形式上非常干净和KKT条件对接也方便。这是数学推导层面达成的共识几乎所有优化工具箱都遵循这个约定。1.2 为什么工程里到处都是二次规划二次规划在工程中的出场率非常高说几个最常见的场景投资组合优化马克维茨均值-方差模型目标是极小化组合风险方差约束是预期收益达标、权重之和为1、不允许卖空等。目标函数天然是二次的约束全是线性的完美匹配QP框架。模型预测控制MPC控制系统在每个控制周期内求解一个带约束的优化问题代价函数通常是状态和控制量的二次型。汽车自适应巡航、无人机轨迹跟踪背后都在解QP。支持向量机SVMSVM的求解可以转化为一个凸二次规划问题对偶形式非常适合用QP求解器处理。最小二乘拟合加约束普通最小二乘是无约束的但如果你需要拟合出的曲线满足单调性、边界范围等条件就变成了带线性约束的二次规划。生活化类比一下二次规划就像是你手里有一堆可调参数目标是一个“成本函数”——这个函数的形状是一个开口向上的碗凸函数你手里的约束是一些直墙线性约束你要做的就是在墙内找到碗底最低的那个点。2. MATLAB里跑通二次规划的工具选择2.1 quadprog函数默认首选MATLAB求解二次规划的核心函数就是quadprog基本调用格式[x, fval, exitflag, output, lambda] quadprog(H, f, A, b, Aeq, beq, lb, ub, x0, options)其中x是最优解fval是最优目标函数值exitflag是求解结束状态大于0代表收敛到最优解output包含迭代信息lambda是拉格朗日乘子用于分析约束的松紧程度。这里有个容易忽略的点H和f缺一不可。如果你的模型没有线性项f也要写成zeros(n,1)不能省略。x0是初始点对于凸二次规划初始点不影响最终结果但会影响迭代步数对于非凸问题x0的选择就比较关键了。2.2 三种算法怎么选quadprog提供了三种算法不同版本略有差异一般R2020a之后版本用interior-point-convex、trust-region-reflective和active-set。算法适用场景优点注意点interior-point-convex大规模凸问题默认推荐内存占用低适合稀疏大矩阵要求H为半正定trust-region-reflective仅有边界约束的问题收敛快适合中小规模要求H为半正定且不能带线性不等式约束active-set中小规模问题冷启动可以处理非凸目标能输出活跃约束迭代较慢内存占用高实际项目中如果没有特殊原因直接用默认的interior-point-convex最省事。若你的H矩阵不是正定矩阵active-set算法有可能仍然给你一个结果但需要自己验证KKT条件。2.3 什么时候该换成fminconquadprog只能处理线性约束和边界约束。一旦遇到非线性约束比如x(1)^2 x(2)^2 1就只能改用fmincon了。但要注意fmincon用于非线性约束优化时默认算法是interior-point它对非凸问题只保证找到局部最优解。如果你的模型本质上是二次规划只是多了几个非线性约束建议先用quadprog求一个无该约束的近似解再把这个解作为fmincon的初值启动这样收敛速度和稳定性都会好很多。3. 手把手从建模到源代码实现3.1 案例一投资组合优化的完整代码假设你有四类资产预期收益率向量r [0.10, 0.12, 0.15, 0.08]协方差矩阵Q已知要求组合的预期收益不低于12%每类资产权重在0到1之间权重之和为1并且权重允许为0即可以完全不配置某类资产。建模如下决策变量x四维向量各资产权重目标函数最小化组合方差x * Q * x约束sum(x) 1r * x 0.120 x 1转成quadprog标准形式H 2 * Q; f zeros(4,1); A -r; % 因为约定条件是 A * x b需要把 变成 b -0.12; Aeq ones(1,4); beq 1; lb zeros(4,1); ub ones(4,1);注意这里A -r原因是quadprog只接受A * x b形式的约束原约束是r * x 0.12两边同时取负号就变成-r * x -0.12。完整源代码% 投资组合二次规划示例 % 目标: 在预期收益不低于12%的前提下最小化组合方差 r [0.10; 0.12; 0.15; 0.08]; % 预期收益率 Q [0.10 0.02 0.04 0.01; 0.02 0.12 0.03 0.02; 0.04 0.03 0.20 0.05; 0.01 0.02 0.05 0.08]; % 协方差矩阵 H 2 * Q; % quadprog目标函数中二次项系数为 0.5*x*H*x f zeros(4, 1); % 没有线性项 A -r; % r*x 0.12 等价于 -r*x -0.12 b -0.12; Aeq ones(1, 4); % 权重之和为1 beq 1; lb zeros(4, 1); % 权重下界 0 ub ones(4, 1); % 权重上界 1 x0 [0.25; 0.25; 0.25; 0.25]; % 初始点 options optimoptions(quadprog, Display, iter, Algorithm, interior-point-convex); [x, fval, exitflag, output, lambda] quadprog(H, f, A, b, Aeq, beq, lb, ub, x0, options); fprintf(最优权重: \n); disp(x); fprintf(最小组合方差: %.6f\n, fval); fprintf(组合预期收益: %.4f\n, r * x); fprintf(exitflag: %d\n, exitflag);3.2 解读运行结果和lambda的含义跑完代码后x的值就是各资产的最优配置比例fval是组合方差的最小值。lambda字段可以帮助你判断哪些约束是“紧”的lambda.ineqlin对应不等式约束的影子价格如果某个值显著大于0说明该约束对结果影响很大。在本例中如果lambda.ineqlin较大说明“收益不低于12%”这条约束是起作用的此时收益约束通常是紧的——最优解刚好卡在12%收益线上。lambda.lower和lambda.upper对应边界约束某个资产权重为0时对应的lambda.lower会大于0说明“不允许卖空”这个下界约束正在限制你。3.3 案例二带等式约束的最小二乘拟合再举一个带等式约束的例子。工程里经常遇到这样的问题用二次函数y a*x^2 b*x c拟合一组数据点同时要求拟合曲线必须通过某个固定点(x0, y0)。这个问题的常规做法是构造线性方程组然后用polyfit但其实也可以建模成二次规划。设决策变量为[a; b; c]目标函数是最小化残差平方和。对数据点(xi, yi)残差为ri a*xi^2 b*xi c - yi目标函数为sum(ri^2)这是一个关于[a;b;c]的二次函数展开后可以写成QP标准形式。约束就是a*x0^2 b*x0 c y0一条线性等式约束。源代码如下% 带固定点约束的二次拟合转化为二次规划 xdata linspace(0, 3, 20); ydata 1.2 * xdata.^2 - 0.5 * xdata 0.8 0.05 * randn(20, 1); % 固定点要求 xfix 1.5; yfix 1.2 * xfix^2 - 0.5 * xfix 0.8; % 假设该点已知 % 构造QP: 决策变量 d [a; b; c] % 残差向量 r_i a*xi^2 b*xi c - yi % 目标 sum(r_i^2) d * M * d - 2 * q * d const Phi [xdata.^2, xdata, ones(size(xdata))]; H 2 * (Phi * Phi); f -2 * (Phi * ydata); Aeq [xfix^2, xfix, 1]; beq yfix; lb [-inf; -inf; -inf]; ub [inf; inf; inf]; options optimoptions(quadprog, Display, off); d quadprog(H, f, [], [], Aeq, beq, lb, ub, [], options); a d(1); b d(2); c d(3); fprintf(拟合系数: a%.4f, b%.4f, c%.4f\n, a, b, c); % 验证固定点 y_fix_fitted a * xfix^2 b * xfix c; fprintf(固定点拟合值: %.4f (真值: %.4f)\n, y_fix_fitted, yfix);这种代码在实际项目里非常实用比如传感器校准、试验数据处理需要强制曲线穿过某些关键点的时候比手动实现拉格朗日乘子法要省事得多。4. 实战中一定会遇到的坑4.1 H矩阵不是正定矩阵怎么办运行quadprog最常遇到的报警或报错就是Hessian矩阵不正定。注意QP的目标函数要存在全局最优解H必须是对称半正定矩阵如果H存在负特征值问题就是非凸的interior-point-convex算法会直接报错告诉你问题不是凸的。处理方案有这么几个检查你的建模是否有误比如二次项系数正负号搞反了目标函数写成了最大化方差结果算出来的H就是负定。加正则化项在H的主对角线上加一个小量epsilon * eye(n)比如epsilon 1e-6让矩阵变成正定。这在处理数据协方差矩阵近似奇异时很常见。改用active-set算法它能处理非凸QP但只能保证找到局部最优解你需要自己对结果做合理性判断。提醒一点如果H不是对称矩阵quadprog在部分版本会直接忽略上三角或下三角只读取一半导致结果和你预期不一致。建议在建模后显式做一次H (H H) / 2;确保对称。4.2 数据量纲差异太大导致收敛极慢如果决策变量里同时存在量级为1e6和1e-6的参数QP求解器在迭代时会出现明显的数值问题——目标函数梯度各分量差异巨大导致收敛缓慢甚至震荡。解决办法是尽量对变量做归一化处理或者调整单位。比如在投资组合例子里如果收益率的量级是0.1而协方差的量级是0.01这些数值都还算健康但如果某个变量的量级达到1e8最好先缩放到1附近算完再还原单位。4.3 exitflag各个取值的含义exitflag是判断求解过程是否成功的第一手信息建议每次跑完都打印出来检查。exitflag含义1函数收敛到解最优性条件满足0迭代次数或时间达到上限结果可能是次优的-2问题不可行即不存在满足所有约束的解-3问题无界目标函数可以无限减小-6非凸问题无法用当前算法求解-8迭代过程中步长过小无法继续推进遇到exitflag -2时先从约束建模入手排查。最常见的原因是约束自相矛盾比如同时要求x 1和x 0或者Aeq矩阵行之间存在线性相关但beq不一致。4.4 大规模问题一定要用稀疏矩阵当决策变量数量达到几万甚至更高时把H和A存储为稠密矩阵会消耗巨大内存计算速度也会骤降。此时务必用sparse构造矩阵把H、A、Aeq都转成稀疏格式quadprog内部会利用稀疏结构加速求解。H_sparse sparse(H); A_sparse sparse(A);一个小经验是如果H矩阵的稠密度超过50%用稀疏格式反而会因为索引开销变慢但工程中的大规模QP往往来源于有限元或网络优化问题H天然是稀疏的用稀疏格式收益非常明显。5. 提升求解效率的几个实用技巧5.1 善用optimoptions微调求解器optimoptions里几个常用参数Display控制命令行输出off关闭输出iter显示每步迭代信息final只显示最终结果。调试阶段建议用iter看收敛趋势上线跑批任务时改成off。OptimalityTolerance最优性容差默认1e-8。工程问题不需要那么高的精度调到1e-6可以明显减少迭代次数。ConstraintTolerance约束容差默认1e-6。如果你的约束条件存在微小违反是可以接受的可以适当放宽。MaxIterations默认值是200对大规模问题可能不够建议根据问题规模调整到1000以上。options optimoptions(quadprog, ... Algorithm, interior-point-convex, ... Display, final, ... OptimalityTolerance, 1e-6, ... ConstraintTolerance, 1e-6, ... MaxIterations, 1000);5.2 warm start加速重复求解MPC这类场景下每个控制周期要解一个类似的QP只是参数在缓慢变化。此时可以用上一次的最优解作为下一次求解的起始点也就是x0显式传入。虽然interior-point-convex算法并不会严格利用初始点做热启动但配合active-set算法时热启动的加速效果非常明显。x_prev x; % 上一次的最优解 [x, fval] quadprog(H, f, A, b, Aeq, beq, lb, ub, x_prev, options);5.3 用KKT条件验证解的合理性机器给的结果未必是物理上可接受的。拿到最优解后建议手动验证一下是否满足所有约束并检查一阶最优性条件% 验证约束 constraint_violation max([A * x - b; Aeq * x - beq; lb - x; x - ub]); fprintf(最大约束违反量: %.2e\n, constraint_violation); % 检查梯度 grad H * x f; fprintf(梯度范数: %.2e\n, norm(grad));如果约束违反量在1e-6量级以下梯度范数也很小基本可以放心使用这个结果。5.4 和符号计算配合做模型验证对维度不高的小问题可以用MATLAB的Symbolic Math Toolbox做一些对照验证。比如随机生成一个H和f用解析法求无约束最优解x -H\f再和quadprog的带约束结果对比看差距是否符合预期。这样能快速发现建模时正负号、上下界搞反的问题。syms x1 x2 real x [x1; x2]; H [4, 1; 1, 3]; f [1; 2]; obj 0.5 * x * H * x f * x; % 求无约束极值 grad gradient(obj, [x1, x2]); sol solve(grad 0, [x1, x2]); disp(double([sol.x1, sol.x2]));6. 几个能直接复用的经验总结第一写出quadprog之前先花五分钟把数学模型在纸上完整写出来标清每一个矩阵维度。我见过大量报错都是因为矩阵维度不匹配A的行数和b的长度对不上或者H不是方阵。纸面上推导一遍大概率能避开这类低级问题。第二目标函数里的0.5并不是可选项。你如果自己手动推导目标函数得到的是x * Q * x传给quadprog时一定要写成H 2 * Q。反过来如果你读别人代码时看到H后面乘了2也就能理解他为什么这么做。第三调试阶段把Display设为iter观察每个迭代步的目标函数下降趋势。如果目标函数反复震荡大概率是约束冲突或者H矩阵数值病态如果一直缓慢下降但迟迟不收敛考虑放宽OptimalityTolerance。第四求解结束后一定要检查exitflag和约束违反量。尤其在做论文或者交付项目时仅仅输出一个x是不完整的把solver返回的验证信息一并记录下来对复现和审计都非常有帮助。第五遇到大规模稀疏QP优先考虑先把模型写对再优化存储结构。过早用稀疏矩阵优化代码如果模型本身不收敛只会让调试难度加倍。二次规划在MATLAB里的实现实际上是一个“建模-转标准形式-调求解器-验证结果”的标准流程掌握了这个流程很多看起来复杂的带约束优化问题都可以按照同样的套路解决。这份源码和踩坑记录希望能帮你少走一些弯路。本文还有配套的精品资源点击获取