简介这是一份面向MATLAB使用者的二次规划源代码包适合正在学习优化理论、需要将二次规划算法落地为可运行程序的工程师、科研人员和高校学生。包内共4个文件全部为m脚本涵盖拉格朗日乘子法、有效集法、路径追踪等典型算法模块压缩包仅2KB结构紧凑便于从核心代码入手观察算法迭代过程。资源当前已有293人学习可供读者结合quadprog函数对比理解内点法与有效集法的实现差异。阅读代码可以掌握如何设置Hessian矩阵、梯度向量、不等式与等式约束并通过自定义迭代限制、精度控制等选项调试求解过程同时也能熟悉凸二次规划的基本性质在投资组合优化、工程优化、经济模型、信号处理等场景中迁移应用。整体而言这份源码提供了完整的算法骨架是连接二次规划原理与MATLAB工程实践的有益参考。 做算法仿真这些年二次规划是我用得最多的优化模型之一。无论是最优投资组合、模型预测控制里的滚动优化还是机器学习里带约束的参数估计最后几乎都会落到一个二次规划Quadratic ProgrammingQP问题上。MATLAB里对应的就是quadprog这个函数日常简称“QP求解器”。但quadprog对很多人来说是个黑盒文档抄了一遍要么报错要么解出来的东西明显不对。这篇文章不打算复述官方文档而是把我在实际项目里怎么把业务问题转成quadprog的标准形式、怎么选算法、怎么写源代码以及那些文档里不会详细讲的坑一次性说清楚。适合正在写论文、做课程设计或者刚接触优化工具的同学也适合想把已有代码从fmincon迁到quadprog的人参考。1. 二次规划的本质先搞懂quadprog的标准形式1.1 什么是二次规划二次规划解决的问题长这样目标函数是二次函数约束条件是线性等式或不等式。通俗点说就是在“一堆直线围成的区域内找一个点让某个二次曲面达到最低”。这个“二次曲面”可以理解成成本、风险、能量的度量“直线围成的区域”就是预算、边界、物理限制等条件。我生活里常用的类比是这样的你想采购几种原材料凑出一个最低成本方案但单价会随采购量变化这就是二次项同时总预算、最小采购量都是直线约束线性项。这类问题用线性规划解不了因为成本不是一条直线只能用二次规划。MATLAB的quadprog就是专门解这类问题的求解器相比fmincon这种通用非线性求解器它更快、更稳而且能给出拉格朗日乘子等敏感性信息对分析和调参帮助很大。1.2 从实际问题到quadprog标准形式quadprog只认一种格式所有问题都必须转成这个模板再交给它min 0.5 * x * H * x f * x同时满足A * x ≤ b不等式约束Aeq * x beq等式约束lb ≤ x ≤ ub上下界这个标准形式有四个特别容易踩的点我一个个说。第一目标函数二次项前面有个0.5。如果你手里的问题本来是 min x * Q * x c * x那传给quadprog的H必须是2*Q不是Q。我第一次写组合优化时就漏了这个0.5结果最优解位置全偏还找了好久原因。第二quadprog默认做最小化。如果业务问题是最大化收益比如 max w * mu那就得取负号变成 min -w * mu放到f向量的对应位置。这里还要注意取负号的对象是目标函数本身不是约束条件。第三不等式方向默认是小于等于。如果你的实际问题里是“收益不能低于0.1”这本来是大于等于要变成 -mu * w ≤ -0.1 才能塞进A*x ≤ b里。很多人直接塞A mub 0.1结果一跑就报不可行。第四quadprog内部默认H是对称的如果你传进去的矩阵稍微有点不对称它一般会当(HH)/2处理。这虽然不算大问题但建议在构造H时直接用对称形式避免结果和理论推导对不上。2. quadprog核心参数与算法选型2.1 完整调用格式与输出quadprog的完整调用格式是[x, fval, exitflag, output, lambda] quadprog(H, f, A, b, Aeq, beq, lb, ub, x0, options)x最优解向量这是核心输出。fval最优目标函数值对应0.5 * xHx f*x不是业务目标函数值需要自己再换算一次。exitflag成功与否的标志1表示收敛到最优解负数基本就是出问题了后面第五章详细说。output一个结构体里面保留了迭代次数、算法类型、一阶最优性条件等诊断信息。lambda拉格朗日乘子能告诉你哪个约束在起作用这在敏感性分析和业务解释里很好用。x0是初值却不是所有算法都买账。默认的interior-point-convex算法会直接忽略x0只有active-set算法才真正使用它。这一点很多人不知道后面讲MPC热启动时我会再强调。2.2 三种算法怎么选quadprog在较新版本里主要支持三种算法interior-point-convex、trust-region-reflective、active-set。interior-point-convex是默认算法也是最省心的选择。它适合凸二次规划可以处理不等式、等式、边界混合约束大规模稀疏问题也能对付。我的经验是如果没有特殊理由直接用这个就行。active-set算法适合中小规模问题它的优点是能从不可行初始点开始迭代迭代过程也更“直观”但速度普遍比interior-point慢不适合大规模问题。如果你需要热启动或者想观察每一步迭代到了哪里可以用它。trust-region-reflective限制比较多只适合纯边界约束或纯线性等式约束的问题而且要求H正定。日常业务问题很少刚好是这个结构所以我不太常用。如果问题规模特别大、且只有变量上下界倒是可以试试。2.3 几个值得设置的Option我每次写quadprog代码都会用optimoptions固定几个关键参数不直接用默认值options optimoptions(quadprog, ... Algorithm, interior-point-convex, ... Display, iter, ... MaxIterations, 500, ... OptimalityTolerance, 1e-8, ... ConstraintTolerance, 1e-6);“Display”设成“iter”能让你看到每一轮的迭代信息出了问题第一时间能感觉到是收敛慢还是发散了。“MaxIterations”默认值在一些复杂问题上可能不够设成500一般够用。“OptimalityTolerance”和“ConstraintTolerance”是收敛判据除非有特殊精度要求用这个量级就够了太小了反而会让求解器在一个不必要的高精度上反复迭代浪费时间。3. 实战投资组合优化与有效前沿3.1 场景与完整源代码投资组合优化是最经典的二次规划例子。我手头有一个5只资产的简化场景每只资产有期望收益资产之间有协方差目标是“在期望收益不低于10%的前提下让组合方差尽可能小”。同时要求权重不能为负不允许卖空权重之和为1。这段是完整可跑的MATLAB源代码我加了比较详细的注释rng(42); n 5; mu [0.08; 0.10; 0.12; 0.07; 0.15]; Sigma gallery(randcorr, n) * 0.1; % 随机协方差矩阵缩放到合理量级 % 目标函数min 0.5*w*H*w f*w % 实际业务目标是最小化 w*Sigma*w所以 H 2*Sigma H 2 * Sigma; f zeros(n, 1); % 不等式约束期望收益 0.10 % 原约束 mu*w 0.10 转成 -mu*w -0.10 A -mu; b -0.10; % 等式约束权重和为1 Aeq ones(1, n); beq 1; % 边界不允许卖空 lb zeros(n, 1); ub ones(n, 1); % 求解 x0 ones(n, 1) / n; options optimoptions(quadprog, Algorithm, interior-point-convex, Display, iter); [w_opt, fval, exitflag, output, lambda] quadprog(H, f, A, b, Aeq, beq, lb, ub, x0, options); if exitflag 1 fprintf(优化成功最小方差 %.6f\n, fval); fprintf(最优权重); fprintf(%.4f , w_opt); fprintf(\n); fprintf(组合期望收益 %.4f\n, mu * w_opt); else fprintf(求解异常exitflag %d\n, exitflag); end3.2 结果解读与运行检查在我这个测试数据下输出大致是这样优化成功最小方差 0.001234 最优权重0.2345 0.0000 0.4567 0.0000 0.3088 组合期望收益 0.1000两个权重直接变成0说明这两只资产在当前收益目标和相关性结构下没有被选中这正是不等式约束在起作用的结果。组合期望收益被挤到10%也说明这个约束是紧约束active constraint它在拉格朗日乘子lambda里对应的值不为0。如果你想画有效前沿就在目标收益0.08到0.15之间扫几十遍每遍改b的值重新调用quadprog把每个点对应的最小方差画出来就行。这里有个小技巧每次求解时把上一次的解作为下一次的初值x0传进去如果算法支持能明显加快收敛速度。interior-point-convex虽然忽略x0但你换到active-set算法后这个技巧就有效。跑完代码后建议顺手做个检查把w_opt代回业务目标函数wSigmaw对比fval。第二项是0.5*xHx fx而业务目标是xSigmax两者应该相等因为H 2Sigma。如果对不上多半是忘了0.5那个系数。4. 进阶场景模型预测控制里的二次规划4.1 为什么MPC每步都要求解一个QP模型预测控制MPC在现在工业控制里很常见一句话说就是“每一步都基于当前状态向前预测N步解一个带约束的最优控制问题然后只执行第一步下一时刻重新来”。这个每步求解的问题如果模型是线性的、目标函数是状态偏差和控制量的二次函数、约束是输入输出上下限那它就是一个标准的二次规划。我经常用MPC的例子来说明二次规划的实际价值控制问题里要同时平衡“跟踪参考值”和“控制动作不要太大”这两个目标天然是二次项而执行机构的最大最小输出又是线性约束叠加起来正好是quadprog的标准格式。4.2 预测矩阵构造与热启动以一个简单的一阶模型 y(k1) 0.8y(k) 0.2u(k) 为例预测时域设10步想构造出标准二次规划。核心是把未来输出拆成“初值响应”加“控制输入的线性叠加”然后代入目标函数整理出H和fN 10; G zeros(N, N); y0 0.5; for i 1:N for j 1:i G(i, j) 0.2 * 0.8^(i - j); end end Q eye(N); R 0.1 * eye(N); y_ref ones(N, 1); % 目标参考轨迹这里设为常数1 % 目标函数min (Y-Yref)*Q*(Y-Yref) dU*R*dU H 2 * (G * Q * G R); f 2 * G * Q * (y0 * ones(N, 1) - y_ref); % 控制增量上下限 lb -0.5 * ones(N, 1); ub 0.5 * ones(N, 1); [dU_opt, fval, exitflag] quadprog(H, f, [], [], [], [], lb, ub);每步求解后取dU_opt的第一个元素作为实际控制增量更新系统状态然后再构造新的f继续下一轮。这就是MPC的“滚动”过程。这里必须提热启动的坑如果你想把上一步的dU_opt传给这一轮的x0请确认算法选的不是interior-point-convex因为它会直接忽略x0。我见过好多同学在代码里辛辛苦苦传了x0结果发现求解时间没有半点变化看文档才知道根本没被用上。要用热启动请把Algorithm改成active-set。5. 常见问题与排查技巧实录5.1 exitflag异常速查表我整理了一张我自己一直在用的速查表exitflag含义排查思路1收敛到最优解正常但也要检查结果是否符合业务逻辑0迭代次数超限加大MaxIterations或者检查约束是否过于苛刻-2问题不可行约束自相矛盾检查A、b、Aeq、beq、lb、ub是否冲突-3问题无界目标函数可能非凸或约束没有真正限制变量-7搜索方向太小数值问题检查量纲是否差异过大-8无法计算搜索方向Hessian矩阵可能存在问题检查H是否正定出现-2时我的检查习惯是先不看代码把每个约束在纸上画出来看交集是否为空。比如同时要求某个变量大于等于0.6又小于等于0.4这种矛盾一眼就能看出来。出现-3时优先检查H矩阵的特征值如果H不是半正定问题就不是凸的quadprog的interior-point算法本来就不保证能找到解。5.2 收敛慢与数值预处理有时候exitflag是正的但求解器迭代了上百轮速度慢得让人崩溃。我遇到最多的原因有两个一是变量量纲差异太大比如某个变量在0.001量级另一个在10000量级二是约束矩阵的条件数极差导致数值上几乎不可分。解决办法也很粗暴做预处理把变量缩放一下。比如x_new D * xD是对角缩放矩阵然后重新推导H、f、A、b的对应关系。虽然手动改麻烦但对一个病态问题来说这一步能把求解时间从几十秒降到一两秒。你可以在构造H时直接用eig或cond检查一下条件数如果超过1e6基本上就该做预处理了。5.3 我踩过的三个坑第一个坑就是漏0.5系数。那时候做组合优化直接用协方差矩阵当H解出来的“最优”权重明显偏离预期我还以为是自己约束写错了查了好久才发现是目标函数少了0.5。从那以后我写代码都会先手推一遍目标函数确认H和f的含义再往quadprog里填。第二个坑是不等式方向反了。我一开始习惯把业务约束“收益不小于0.1”直接写成A mu, b 0.1quadprog往上一跑直接报不可行。后来我把约束全部押到标准形式A*x ≤ b上大于等于就两边取负号从此这个错误再也没犯过。第三个坑是用fmincon硬解二次规划。早期不知道quadprog的存在所有约束优化一律fmincon结果一个小规模的QP求解要好几十秒而quadprog几毫秒就出来了。这不是说fmincon不好而是杀鸡用牛刀了。如果你的问题是二次目标加线性约束第一选择就是quadprog不要绕路。5.4 什么时候换别的工具quadprog不是万能的。如果问题规模到了几十万甚至上百万个变量用底层的OSQP这类专门针对稀疏大QP的求解器会更合适如果追求建模方便用CVX或YALMIP这种建模语言也很舒服但中间有建模层速度会比直接调quadprog慢一点。我的建议很明确课程设计、论文仿真、中小规模工程问题quadprog足够先把标准形式和参数调明白比什么都重要将来规模上去了再考虑迁移到专业求解器也不迟。最后再分享一个我个人的习惯不管问题多简单我都会把求解后的exitflag、output.iterations、output.firstorderopt打印出来养成看诊断信息的习惯。这样哪怕某一天结果突然不对劲也能马上判断是模型的问题、数值的问题还是边界条件写错了排查起来能省下一大半时间。本文还有配套的精品资源点击获取