【最优化理论】序列二次规划SQP理解和例题代码

前言

最近最优化课程要求完成序列二次规划SQP算法，完成后为了加深理解写了这篇博客。

序列二次规划SQP

序列二次规划主要用于约束优化问题，其是目前公认的求解约束非线性优化问题最有效的方法之一。与其他算法相比，序列二次规划法的优点是收敛性好、计算效率高、边界搜索能力强，但是其在问题规模比较大时，计算复杂度比较高。

KKT条件

为了解决带不等式的约束优化问题，我们先要了解KKT条件。
$$\begin{gathered} minf(x) \\ s.t.g_i(x)\le 0i=1,...,m \\ {h}_j(x)=0j=1,...,n \\ ⟹\text{KKT条件}\{\begin{array}{l}\nabla f(x^{\ast })+{\lambda }_i{\sum }_{i=1}^m\nabla g_i(x^{\ast })+{\mu }_j{\sum }_{j=1}^n\nabla h_j(x^{\ast })=0(1)\\ {\lambda }_i{g}_i(x^{\ast })=0(2)\\ h_j(x^{\ast })=0(3)\\ {\lambda }_i\ge 0(4)\\ g_i(x^{\ast })\le 0(5)\end{array} \end{gathered}$$
定义的拉格朗日函数为：
L ( x , { λ i } , { μ i } ) = f ( x ) + λ i ∑ i = 1 m g i ( x ) + μ j ∑ j = 1 n h j ( x ) L(x,\{\lambda_i\},\{\mu_i\})=f(x)+\lambda_i \sum_{i=1}^{m} g_i(x)+\mu_j \sum_{j=1}^{n} h_j(x) L(x,{λi​},{μi​})=f(x)+λi​i=1∑m​gi​(x)+μj​j=1∑n​hj​(x)
其中，$x^{\ast }$是我们要找的最优解，我们可以利用式子(1)(2)(3)帮助我们找到最优解，找到最优解后再代入式子(4)(5)进行验证，不成立就排除。

对于普遍问题，KKT条件是判断某点是极值点的必要条件，也就是说，最优解$x^{\ast }$一定满足KKT条件，但是KKT条件的解不一定是最优解。

对于凸规划，KKT条件才是判断某点是极值点的充要条件。
详细的解释请看这里

最优化问题

假设我们的最优化问题为：

$$\begin{gathered} minf(x) \\ s.t.h(x)=0 \\ g(x)\le 0 \end{gathered}$$

对于带约束的问题，我们可以使用拉格朗日乘数法，构造一个拉格朗日函数：
$$L(x,\lambda ,\mu )=f(x)+\lambda g(x)+\mu h(x)$$

利用KKT条件

利用我们之前说到的KKT条件，我们可以得到
$$⟹\text{KKT条件}\{\begin{array}{l}\nabla L(x,\lambda ,\mu )=\nabla f(x)+\lambda \nabla g(x)+\mu \nabla h(x)=0\\ \lambda g(x)=0\\ h(x)=0\\ \lambda \ge 0\\ g(x)\le 0\end{array}$$
我们之前说到的，求该含有不等式约束问题的最优解可以转化为求KKT条件的最优解。
那么实际上我们就是要解方程组
$F(x,\lambda ,\mu )=\left[\begin{array}{c}\nabla L(x,\lambda ,\mu )\\ \lambda g(x)\\ h(x)\end{array}\right]=0$

我们设它的雅克比矩阵为
$J_F(x,\lambda ,\mu )$
要解这个方程组，我们可以使用牛顿法
$x_{k+1}=x_k+\frac{f(x_k)}{f(x_k^{{}^{\prime }})}$
在多元函数中，我们不能直接这样解，于是我们有
$$\left(\begin{array}{c}{x}_{k+1}\\ {\lambda }_{k+1}\\ {\mu }_{k+1}\end{array}\right)=\left(\begin{array}{c}{x}_k\\ {\lambda }_k\\ {\mu }_k\end{array}\right)-J_F^{-1}(x,\lambda ,\mu )F(x,\lambda ,\mu )$$
这里注意：$A^{-1}B$实际上就是$Ax=B$的解，我们就不用求逆，可以直接解$Ax=B$这个方程组
也就是$J_F(x,\lambda ,\mu )d_k$$=F(x,\lambda ,\mu )$
解出这个$d_k$，不断迭代，我们就能够求出最优解。
$$\left(\begin{array}{c}{x}_{k+1}\\ {\lambda }_{k+1}\\ {\mu }_{k+1}\end{array}\right)=\left(\begin{array}{c}{x}_k\\ {\lambda }_k\\ {\mu }_k\end{array}\right)-d_k$$
但是可以看到这个方程，要解出它是比较困难的。

线性化

我们说到，SQP算法适用于解决带有非线性约束的问题，上述的h(x)和g(x)都有可能是非线性的，为了更方便求解，我们需要将上述KKT条件中的非线性方程线性化。
$$\begin{gathered} \nabla L(x,\lambda ,\mu )=\nabla f(x)+\left[\begin{array}{cc}\nabla g(x)& \nabla h(x)\end{array}\right]\left[\begin{array}{c}\lambda \\ \mu \end{array}\right] \\ \lambda g(x)=0 \\ h(x)=0 \\ g(x)\le 0 \end{gathered}$$
对上面的$\nabla L(x,\lambda ,\mu )$、$g(x)$和$h(x)$式子泰勒展开得到
$$\begin{gathered} \nabla L(x,\lambda ,\mu )\approx \nabla L(x_k,{\lambda }_k,{\mu }_k)+\left[\begin{array}{ccc}{\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k)& {\nabla }_{x\lambda }^{2}L(x_k,{\lambda }_k,{\mu }_k)& {\nabla }_{x\mu }^{2}L(x_k,{\lambda }_k,{\mu }_k)\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta \mu \end{array}\right] \\ h(x)\approx h(x_k)+{\nabla }_{x}h(x_k)\Delta x \\ g(x)\approx g(x_k)+{\nabla }_{x}g(x_k)\Delta x \end{gathered}$$
所以原问题的KKT条件就会变成
$$\begin{gathered} \nabla L(x_k,{\lambda }_k,{\mu }_k)+\left[\begin{array}{ccc}{\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k)& {\nabla }_{x\lambda }^{2}L(x_k,{\lambda }_k,{\mu }_k)& {\nabla }_{x\mu }^{2}L(x_k,{\lambda }_k,{\mu }_k)\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta \mu \end{array}\right]=0 \\ h(x_k)+{\nabla }_{x}h(x_k)\Delta x=0 \\ g(x_k)+{\nabla }_{x}g(x_k)\Delta x\le 0 \\ \lambda (g(x_k)+{\nabla }_{x}g(x_k)\Delta x)=0 \\ \lambda \ge 0 \end{gathered}$$
把第一个式子中的各个部分列一下：
$$\begin{gathered} \nabla L(x,\lambda ,\mu )=\nabla f(x)+\left[\begin{array}{cc}\nabla g(x)& \nabla h(x)\end{array}\right]\left[\begin{array}{c}\lambda \\ \mu \end{array}\right] \\ {\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k)={\nabla }_{xx}^{2}f(x)+\left[\begin{array}{cc}{\lambda }_k& {\mu }_k\end{array}\right]\left[\begin{array}{c}{\nabla }_{xx}^{2}h(x)\\ {\nabla }_{xx}^{2}g(x)\end{array}\right] \\ {\nabla }_{x\lambda }^{2}L(x_k,{\lambda }_k,{\mu }_k)=\nabla g(x) \\ {\nabla }_{x\mu }^{2}L(x_k,{\lambda }_k,{\mu }_k)=\nabla h(x) \end{gathered}$$

但是我们希望得到KKT条件的求解迭代方向，我们将这些代入KKT条件的第一个式子得到
$\nabla L(x_k,{\lambda }_k,{\mu }_k)+\left[\begin{array}{ccc}{\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k)& {\nabla }_{x\lambda }^{2}L(x_k,{\lambda }_k,{\mu }_k)& {\nabla }_{x\mu }^{2}L(x_k,{\lambda }_k,{\mu }_k)\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta \mu \end{array}\right]$
$={\nabla }_{x}f(x)+\left[\begin{array}{cc}{\nabla }_{x}g(x)& {\nabla }_{x}h(x)\end{array}\right]\left[\begin{array}{c}{\lambda }_k\\ {\mu }_k\end{array}\right]+\left[\begin{array}{c}{\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k),{\nabla }_{x}g(x),{\nabla }_{x}h(x)\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta \mu \end{array}\right]$
$={\nabla }_{x}f(x)+\left[\begin{array}{cc}{\nabla }_{x}g(x)& {\nabla }_{x}h(x)\end{array}\right]\left[\begin{array}{c}{\lambda }_k\\ {\mu }_k\end{array}\right]+{\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k)\Delta x+\left[\begin{array}{c}{\nabla }_{x}g(x),{\nabla }_{x}h(x)\end{array}\right]\left[\begin{array}{c}\Delta \lambda \\ \Delta \mu \end{array}\right]$
$={\nabla }_{x}f(x)+{\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k)\Delta x+\left[\begin{array}{cc}{\nabla }_{x}g(x)& {\nabla }_{x}h(x)\end{array}\right]\left[\begin{array}{c}{\lambda }_k\\ {\mu }_k\end{array}\right]+\left[\begin{array}{c}{\nabla }_{x}g(x),{\nabla }_{x}h(x)\end{array}\right]\left[\begin{array}{c}\Delta \lambda \\ \Delta \mu \end{array}\right]$
$={\nabla }_{x}f(x)+{\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k)\Delta x+\left[\begin{array}{c}{\nabla }_{x}g(x),{\nabla }_{x}h(x)\end{array}\right]\left[\begin{array}{c}{\lambda }_{k+1}\\ {\mu }_{k+1}\end{array}\right]$
这样，我们就找出了参数迭代的关系式。

QP二次规划子问题

说了那么多，为什么SQP叫序列二次规划，这里我们就能理解。
我们发现上面这个被线性化的KKT条件，刚好是下面这个最优化问题的KKT条件
$$\begin{gathered} min\frac{1}{2}\Delta x^T{\nabla }_{xx}^{2}L(x_k,{\lambda }_k,{\mu }_k)\Delta x+\Delta x^T{\nabla }_{x}f(x_k) \\ s.t.h(x_k)+{\nabla }_{x}h(x_k)\Delta x=0 \\ g(x_k)+{\nabla }_{x}g(x_k)\Delta x\le 0 \end{gathered}$$
我们需要优化的变量为$\Delta x$，而上面这个问题刚好就是一个QP问题，这意味着，我们每次求迭代的方向，实际上就是解一个QP问题，现在明白为什么叫SQP问题了吧？因为每找到最优解的过程就是一连串地解QP问题的过程！
这里解QP问题的方法就不详细介绍了(因为来不及了，后面有时间我再补上)，可以用内点法等等。
matlab有QP问题求解的函数库叫quadprog，后面的代码我就要成为调库侠了！
经过QP问题求解，我们就能完成一轮迭代，得到
$x=x_k+\Delta x$
$\lambda ={\lambda }_{k+1}$

代码实现

本文的最优化问题为：
$$\begin{gathered} minf(x)=(x_1-1{)}^2+(x_2-2{)}^2 \\ s.t.x_1^2-x_2=0 \\ {x}_1^2+x_2^2\le 25 \end{gathered}$$
本人没怎么学过matlab，matlab代码真的不是很规范。。。

clear;clc;
syms x1 x2 lambda mu

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});

% 等式约束
Aeq = [1, -1];
beq = 0;


% 不等式约束
A = [1, 1];
b = 25;

% 约束函数
[ceq,c] = nonlcon([x1;x2], Aeq, beq, A, b);
disp(ceq)
disp(c)
ceq_fun = matlabFunction(ceq, 'Vars', {x1, x2});
c_fun = matlabFunction(c, 'Vars', {x1, x2});
% 等式约束求导
d_ceq=jacobian(ceq, [x1, x2])';
dd_ceq= hessian(ceq, [x1, x2])';

% 不等式约束求导
d_c=jacobian(c, [x1, x2])';
dd_c= hessian(c, [x1, x2])';

d_ceq_fun = matlabFunction(d_ceq, 'Vars', {x1, x2});
dd_ceq_fun = matlabFunction(dd_ceq, 'Vars', {x1, x2});
d_c_fun = matlabFunction(d_c, 'Vars', {x1, x2});
dd_c_fun = matlabFunction(dd_c, 'Vars', {x1, x2});

% 拉格朗日函数
L = f_sym + mu*ceq+lambda*c;
% 计算拉格朗日函数二阶偏导数
dd_L = dd_f+mu*dd_ceq+lambda*dd_c;
L_fun = matlabFunction(L, 'Vars', {x1, x2,lambda, mu});
dd_L_fun = matlabFunction(dd_L, 'Vars', {x1, x2,lambda, mu});


%初始猜想
max_iter = 50;
x0 = [5; 25];
x_k=x0;
lambda_k =6;
mu_k =1;
x1=x0(1);
x2=x0(2);

tol=0.001;
k=0;

path = []; 

for k = 1:max_iter
    % 把数值代入所有符号表达式
    f_val = f(x1,x2);
    d_f_val=d_f_fun(x1,x2);
    dd_f_val=dd_f_fun(x1,x2);
    
    % 不等式约束
    c_val=c_fun(x1,x2);
    d_c_val=d_c_fun(x1,x2);
    dd_c_val=dd_c_fun(x1,x2);
    % 等式约束
    ceq_val=ceq_fun(x1,x2);
    d_ceq_val=d_ceq_fun(x1,x2);
    dd_ceq_val=dd_ceq_fun(x1,x2);
    
    
    L_val=L_fun(x1,x2,lambda_k, mu_k);
    dd_L_val=dd_L_fun(x1,x2,lambda_k, mu_k);
    
    % 求二次规划问题
    % 公式为= ?x'?f(x)+(1/2)?x'??L?x
    % 等式约束为hj(x) + ?hj(x)'?x = 0
    % 不等式约束为gi(x) + ?gi(x)'?x = 0
    [delta_x,fval,exitflag,output,lambda] = quadprog(dd_L_val,d_f_val,d_c_val',-c_val,d_ceq_val', -ceq_val, [], [],x0);
    x_k = x_k + delta_x;
    
    x1=x_k(1);
    x2=x_k(2);
    disp("当前迭代点为:")
    disp(x_k);
    lambda_k = lambda.ineqlin;
    mu_k = lambda.eqlin;
    % 检查收敛
    if norm(delta_x) < tol
        break;
    end
    
    path = [path; x1, x2];
end

% 输出最终结果
disp('最优解:');
disp(x_k);
disp('计算x_k:');
disp(f(2.1272,4.5249));
disp('最优值:');
disp(f_val);


% 创建网格点以生成图像
x1 = linspace(-5, 5, 20);
x2 = linspace(-5, 5, 20);
[X1, X2] = meshgrid(x1, x2);
Z = f(X1, X2);
 
% 绘制3D图像
figure;
surf(X1, X2, Z, 'FaceAlpha', 0.7); % 函数图像
hold on;

xlabel('x1');
ylabel('x2');
zlabel('f(x)');

% 定义 x3 的范围
x3 = linspace(-30, 30, 20);
% 创建网格
[X1, X3] = meshgrid(x1, x3);
% 等式约束 x1^2 - x2 = 0
X2_eq = min(max(X1.^2,-5),5);
surf(X1, X2_eq, X3, 'FaceAlpha', 0.5, 'EdgeColor', 'none'); % 绘制平面

% 不等式约束就是个球体
radius = 5;
[X, Y, Z] = sphere;
X = radius * X;
Y = radius * Y;
Z = radius * Z;
% 绘制球体
surf(X, Y, Z,'FaceAlpha', 0.25,'EdgeColor', 'none');


% 绘制SQP算法下降路径
plot3(path(:,1), path(:,2), f(path(:,1), path(:,2)), 'r-o', 'LineWidth', 2, 'MarkerSize', 8, 'MarkerFaceColor', 'r');

function [ceq, c] = nonlcon(x, Aeq, beq, A, b)
    % 等式约束 h1(x) = x1^2 - x2 = 0
    ceq = Aeq(1)*x(1)^2 + Aeq(2)*x(2) - beq;
    % 不等式约束 g1(x) = x1 <= 5
    % c = A*x - b;
    % 不等式约束 g1(x) = x1^2 + x2^2 <= 25
    c = A(1)*x(1)^2+A(2)*x(2)^2 - b;
end

图像结果：

要在等式条件构成的抛物面上，同时也要在不等式约束球体与目标函数交线（白色的线）内。

总结

看了这么久的SQP，我个人的理解是SQP问题本质上就是将带非线性约束的问题线性化，然后将每次迭代转化为求一个QP问题，经过求解一系列的QP问题，最终我们能够得到最优解。
这篇博客只是个人的一次记录，matlab代码写得比较差，内容写得比较急可能也会有错，如果有错欢迎指出。

参考

[1]KKT条件，原来如此简单 | 理论+算例实践
[2]最优化抄书笔记：序列二次规划
[3]【最优化】序列（逐步）二次规划法（SQP）
[4]手动构建SQP求解NMPC问题（一）