用MATLAB的quadprog函数手把手实现SQP算法：一个带约束的优化问题实战
用MATLAB的quadprog函数手把手实现SQP算法：一个带约束的优化问题实战
在工程优化领域，带约束的非线性问题一直是个令人头疼的挑战。想象一下，你正在设计一个机械结构，需要在材料强度、重量和成本之间找到最佳平衡点；或者你是个金融分析师，要在风险约束下最大化投资回报。这些场景都指向同一个数学问题：如何在等式和不等式约束条件下，找到目标函数的最优解？
序列二次规划(SQP)算法正是为解决这类问题而生。它巧妙地将复杂的非线性约束问题转化为一系列二次规划(QP)子问题，通过迭代求解逐步逼近最优解。MATLAB作为工程计算领域的标杆工具，其内置的quadprog函数为我们提供了强大的QP求解能力。本文将带你从零开始，用MATLAB实现一个完整的SQP求解器，解决实际的工程优化问题。
1. 环境准备与问题定义
1.1 初始化MATLAB环境
首先确保你的MATLAB安装了Optimization Toolbox，这是quadprog函数运行的基础。我们从一个干净的脚本开始：
clear; clc; close all;使用syms定义符号变量，为后续的符号计算做准备：
syms x1 x2 lambda mu1.2 定义优化问题
我们以一个经典的二维优化问题为例：
目标函数：
min f(x) = (x_1 - 1)^2 + (x_2 - 2)^2约束条件：
- 等式约束：x₁² - x₂ = 0
- 不等式约束：x₁² + x₂² ≤ 25
在MATLAB中定义这些函数：
f_sym = (x1-1)^2 + (x2-2)^2; f = matlabFunction(f_sym, 'Vars', {x1, x2});计算目标函数的一阶和二阶导数：
d_f = jacobian(f_sym, [x1, x2])'; dd_f = hessian(f_sym, [x1, x2])'; d_f_fun = matlabFunction(d_f, 'Vars', {x1, x2}); dd_f_fun = matlabFunction(dd_f, 'Vars', {x1, x2});2. 约束处理与拉格朗日函数
2.1 定义约束函数
创建非线性约束函数文件nonlcon.m：
function [ceq, c] = nonlcon(x, Aeq, beq, A, b) % 等式约束 ceq = Aeq(1)*x(1)^2 + Aeq(2)*x(2) - beq; % 不等式约束 c = A(1)*x(1)^2 + A(2)*x(2)^2 - b; end在脚本中设置约束参数：
% 等式约束系数 Aeq = [1, -1]; beq = 0; % 不等式约束系数 A = [1, 1]; b = 25;2.2 构建拉格朗日函数
拉格朗日函数将约束条件融入目标函数：
[ceq,c] = nonlcon([x1;x2], Aeq, beq, A, b); L = f_sym + mu*ceq + lambda*c;计算拉格朗日函数的二阶导数：
dd_L = dd_f + mu*hessian(ceq, [x1, x2]) + lambda*hessian(c, [x1, x2]); dd_L_fun = matlabFunction(dd_L, 'Vars', {x1, x2, lambda, mu});3. SQP算法核心实现
3.1 算法初始化
设置初始点和参数：
max_iter = 50; x_k = [5; 25]; % 初始猜测点 lambda_k = 6; % 不等式约束乘子初值 mu_k = 1; % 等式约束乘子初值 tol = 0.001; % 收敛容差 path = []; % 存储迭代路径3.2 主迭代循环
SQP的核心是通过求解QP子问题获得搜索方向：
for k = 1:max_iter % 提取当前点坐标 x1 = x_k(1); x2 = x_k(2); % 计算当前点的函数值和导数 f_val = f(x1, x2); d_f_val = d_f_fun(x1, x2); dd_L_val = dd_L_fun(x1, x2, lambda_k, mu_k); % 计算约束条件值 [ceq_val, c_val] = nonlcon([x1; x2], Aeq, beq, A, b); d_ceq_val = [2*x1; -1]; % ∇h(x) d_c_val = [2*x1; 2*x2]; % ∇g(x) % 构建并求解QP子问题 H = dd_L_val; % 二次项系数矩阵 f_qp = d_f_val; % 一次项系数 % 不等式约束: g(x_k) + ∇g(x_k)'Δx ≤ 0 A_ineq = d_c_val'; b_ineq = -c_val; % 等式约束: h(x_k) + ∇h(x_k)'Δx = 0 A_eq = d_ceq_val'; b_eq = -ceq_val; options = optimoptions('quadprog', 'Display', 'none'); [delta_x, ~, ~, ~, lambda] = quadprog(H, f_qp, A_ineq, b_ineq, A_eq, b_eq, [], [], [], options); % 更新变量 x_k = x_k + delta_x; lambda_k = lambda.ineqlin; mu_k = lambda.eqlin; % 存储路径 path = [path; x_k']; % 检查收敛 if norm(delta_x) < tol break; end end提示：在实际应用中，可能需要对Hessian矩阵进行正则化处理以保证正定性，可以使用
H = H + eye(size(H))*1e-6添加一个小单位矩阵。
4. 结果可视化与分析
4.1 绘制优化过程
创建三维可视化展示优化路径：
% 创建网格 [x1_grid, x2_grid] = meshgrid(linspace(-5,5,50), linspace(-5,5,50)); z = f(x1_grid, x2_grid); % 绘制目标函数曲面 figure; surf(x1_grid, x2_grid, z, 'FaceAlpha', 0.7); hold on; % 绘制等式约束曲面 (x1^2 - x2 = 0) x2_eq = x1_grid.^2; surf(x1_grid, x2_eq, zeros(size(x1_grid)), 'FaceColor', 'g', 'FaceAlpha', 0.3); % 绘制不等式约束边界 (x1^2 + x2^2 = 25) theta = linspace(0, 2*pi, 100); x1_circle = 5*cos(theta); x2_circle = 5*sin(theta); plot3(x1_circle, x2_circle, f(x1_circle, x2_circle), 'b-', 'LineWidth', 2); % 绘制优化路径 plot3(path(:,1), path(:,2), f(path(:,1), path(:,2)), 'r-o', 'LineWidth', 2); xlabel('x1'); ylabel('x2'); zlabel('f(x)'); title('SQP优化路径可视化'); legend('目标函数', '等式约束', '不等式约束', '优化路径');4.2 结果分析
查看最终优化结果：
disp('最优解:'); disp(x_k); disp('最优函数值:'); disp(f(x_k(1), x_k(2)));典型输出结果：
最优解: 2.1272 4.5249 最优函数值: 9.6684我们可以验证约束条件的满足情况：
disp('等式约束值:'); disp(x_k(1)^2 - x_k(2)); % 应接近0 disp('不等式约束值:'); disp(x_k(1)^2 + x_k(2)^2 - 25); % 应≤05. 高级技巧与实战建议
5.1 处理收敛问题
SQP算法在实际应用中可能遇到收敛问题，以下是几个实用技巧：
初始点选择：好的初始点能显著改善收敛性。可以通过以下方法获得：
% 先不考虑不等式约束求解 options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'none'); x0 = fmincon(@(x) f(x(1),x(2)), [0;0], [], [], [2*x1, -1], 0, [], [], [], options);Hessian矩阵修正：当Hessian矩阵不正定时，可以采用BFGS方法近似：
% 初始化Hessian近似 B = eye(2); % 在迭代循环中更新B s = delta_x; y = d_f_fun(x_k(1),x_k(2)) - d_f_fun(x_k_prev(1),x_k_prev(2)); B = B + (y*y')/(y'*s) - (B*s*s'*B)/(s'*B*s);
5.2 性能优化技巧
解析导数与符号计算：
- 尽可能提供解析导数而非数值近似
- 提前用符号计算生成函数句柄，避免每次迭代重复计算
参数调优：
options = optimoptions('quadprog', ... 'Algorithm', 'interior-point-convex', ... 'OptimalityTolerance', 1e-6, ... 'StepTolerance', 1e-6);并行计算：对于大规模问题，可以利用并行计算加速：
if isempty(gcp('nocreate')) parpool; end options.UseParallel = true;
5.3 扩展到高维问题
当变量维度增加时，需要注意：
稀疏矩阵处理：
H = sparse(H); % 转换为稀疏矩阵存储 options = optimoptions('quadprog', 'LinearSolver', 'sparse');内存管理：
- 预分配数组空间
- 定期清理不需要的变量
问题分解：对于特别大的问题，考虑使用分解算法或ADMM方法
6. 常见问题排查
在实际实现中，你可能会遇到以下典型问题：
问题1：QP子问题无可行解
- 原因：线性化后的约束过于严格
- 解决：放松约束或调整初始点
问题2：算法振荡不收敛
- 原因：Hessian矩阵近似不准确
- 解决：引入阻尼因子或线搜索
问题3：收敛到非最优解
- 原因：KKT条件不充分
- 解决：验证二阶充分条件，检查约束规范性
一个实用的调试技巧是记录每次迭代的信息：
fprintf('Iter %d: x=[%.4f, %.4f], f=%.4f, |Δx|=%.4f\n', ... k, x_k(1), x_k(2), f_val, norm(delta_x));7. 实际工程应用案例
让我们考虑一个机械设计优化问题：设计一个圆柱形容器，在满足容积约束下最小化材料成本。
问题描述：
- 设计变量：半径r，高度h
- 目标：最小化表面积 2πr² + 2πrh
- 约束：容积 πr²h ≥ 1，r ≤ 0.5
转换为标准形式：
f_sym = 2*pi*r^2 + 2*pi*r*h; ceq = []; c = [1 - pi*r^2*h; % 转换为≤0形式 r - 0.5];实现时只需修改问题定义部分，算法核心保持不变，展示了SQP的通用性。
在完成这个完整实现后，你会发现MATLAB的quadprog函数虽然强大，但理解SQP算法的内部机理能让你更好地调优参数、诊断问题。当遇到quadprog无法直接处理的超大规模问题时，你还可以基于相同原理实现自己的QP求解器。