二次规划--积极集法（active set method）

wsh_

已于 2024-05-02 11:27:17 修改

阅读量8.7k

点赞数 6

分类专栏：最优化文章标签：二次规划不等式约束

于 2021-01-08 23:07:41 首次发布

本文链接：https://blog.csdn.net/hxudhdjuf/article/details/112374967

版权

最优化专栏收录该内容

1 篇文章 0 订阅

订阅专栏

本文介绍了如何使用有效集法解决带有等式和不等式约束的二次规划问题。首先阐述了拉格朗日法解决等式约束问题，然后详细解释了有效集法的原理和步骤，包括四种迭代情形。最后，给出了一个具体的MATLAB代码示例来演示有效集法的实现过程。

摘要由CSDN通过智能技术生成

1.等式约束二次规划问题

$\left\{\begin{array}{ll} \min & f(\boldsymbol{x})=\frac{1}{2} \boldsymbol{x}^{\mathrm{T}} \boldsymbol{G} \boldsymbol{x}+\boldsymbol{d}^{\mathrm{T}} \boldsymbol{x} \\ \text { s.t. } & \boldsymbol{A} \boldsymbol{x}=\boldsymbol{b} \end{array}\right.$ (1)

其中 $\boldsymbol{G}\subseteq \mathbb{R}^{n \times n}$ 且为对称矩阵， $\boldsymbol d\subseteq \mathbb{R}^n$ , $\boldsymbol{A}\subseteq \mathbb{R}^{m \times n}$ ， $\boldsymbol b\subseteq \mathbb{R}^n$ . 求解等式约束二次规划问题一般有直接消元法，正交分解法和拉格朗日法。

拉格朗日法

构建Eq.(1)的拉格朗日乘子函数，得

$L(x,\lambda)=\frac{1}{2}x^{\rm T}Gx+d^{\rm T}x + \lambda^{\rm T}(Ax-b)$

由KKT条件可得：

$\left\{ {\begin{array}{*{20}{c}} {{\nabla _{\boldsymbol{x}}}L({\boldsymbol{x}},{\boldsymbol{\lambda }}) = {\boldsymbol{0}}} \\ {{\nabla _{\boldsymbol{\lambda }}}L({\boldsymbol{x}},{\boldsymbol{\lambda }}) = {\boldsymbol{0}}} \end{array}} \right. \Rightarrow \left\{ {\begin{array}{*{20}{l}} {{\boldsymbol{Gx}} + {\boldsymbol{d}} + {{\boldsymbol{A}}^{\text{T}}}{\boldsymbol{\lambda }} = {\boldsymbol{0}}} \\ {{{\boldsymbol{A}}}{\boldsymbol{x}} - {\boldsymbol{b}} = {\boldsymbol{0}}} \end{array}} \right.$

写成矩阵的形式得

$\left[\begin{array}{cc} \boldsymbol{G} & \boldsymbol{A}^{\mathrm{T}} \\ \boldsymbol{A} & \boldsymbol{0} \end{array}\right]\left[\begin{array}{l} \boldsymbol{x} \\ \lambda \end{array}\right]=\left[\begin{array}{l} -\boldsymbol{d} \\ \boldsymbol{b} \end{array}\right]$

系数矩阵 $\left[\begin{array}{cc} \boldsymbol{G} & \boldsymbol{A}^{\mathrm{T}} \\ \boldsymbol{A} & \boldsymbol{0} \end{array}\right]$ 称为Lagrange矩阵，它是对称不定的。如果Lagrange矩阵的逆存在

可得最优解为 $\left\{\begin{array}{l} x^{*}=-Q d+R b \\ \lambda^{*}=-R^{T} d+U b \end{array}\right.$ .如果 $\boldsymbol{G}^{-1}$ 存在，可得

$\left\{\begin{array}{l} \boldsymbol{Q}=\boldsymbol{G}^{-1}-\boldsymbol{G}^{-1} \boldsymbol{A}^{\mathrm{T}}\left(\boldsymbol{A} \boldsymbol{G}^{-1} \boldsymbol{A}^{\mathrm{T}}\right)^{-1} \boldsymbol{A} \boldsymbol{G}^{-1} \\ \boldsymbol{R}=\boldsymbol{G}^{-1} \boldsymbol{A}^{\mathrm{T}}\left(\boldsymbol{A} \boldsymbol{G}^{-1} \boldsymbol{A}^{\mathrm{T}}\right)^{-1} \\ \boldsymbol{U}=-\left(\boldsymbol{A} \boldsymbol{G}^{-1} \boldsymbol{A}^{\mathrm{T}}\right)^{-1} \end{array}\right.$

如果Lagrange矩阵的逆不存在，可以求伪逆。

2.等式和不等式约束二次规划问题

$\left\{\begin{array}{ll} \min & f(\boldsymbol{x})=\frac{1}{2} \boldsymbol{x}^{\mathrm{T}} \boldsymbol{G} \boldsymbol{x}+\boldsymbol{d}^{\mathrm{T}} \boldsymbol{x} \\ \text { s.t. } & \boldsymbol{A_e} \boldsymbol{x}=\boldsymbol{b_e} \\ &\boldsymbol{A_i} \boldsymbol{x}\leqslant \boldsymbol{b_i} \end{array}\right.$

由于带有等式约束的二次规划问题已经有成熟的解法，所以要求解带有不等式约束的二次规划问题，一种很自然的想法就是将不等式约束转换成等式约束。有效集法（Active Set Method）是求解带不等式约束的二次规划问题的一种经典算法。首先看下问题的几何表述

preview

可行域：目标问题的可行域是空间中的一个凸多面体（包括边界和内部）。例如上图中的绿色多边形区域。
等值面：如果矩阵G正定，那么目标函数的等值面是空间中的一族同心超椭球面，所有椭球面相似且同向（各轴同向重合）。例如上图中的同心椭圆曲线族。越内层的椭球面，目标函数在其上的取值越优（低）。
优化问题：在几何表述下，我们的目标就是在可行域内找到达到同心椭圆曲线族最内层的点。例如上图中，P是无约束最优解，Q是约束最优解。

显然加入约束之后，最优解必定出现在某个或某几个不等式约束的边界上。

有效集：任给一个可行解，其有效集就是满足等号成立的约束条件的集合。显然，有效集必然包含所有等式约束，同时包含不等式约束的一个子集。注意这里是等号成立，也就是说，当前解让一个不等式约束的等号成立时，这个不等式约束才会被纳入有效集。
最优有效集：最优解的有效集称为最优有效集。只需把最优有效集中的约束全部改写为等式约束，并扔掉其它约束，然后用Lagrange乘子法直接求解即可。

由于约束条件数目是有限的，而最优有效集是其子集，进而最优有效集只有有限多种可能。于是很自然地就想到了下面大名鼎鼎的暴力破解法——穷举法。

每次，我们从全体约束中选出一个子集（包含全部等式约束和部分不等式约束），称之为本次尝试的工作集，并求解该工作集对应的子优化问题。一个工作集对应的子优化问题是这样定义的：在原问题中，将工作集中的约束全部改为等式约束，同时扔掉工作集之外的约束。由于这是一个等式约束问题，所以可以用Lagrange法轻松求解。之后检查该解是否是原问题的可行解（即它是否也满足工作集之外的约束）。如果是，则记录下来，否则就扔掉。遍历全部可能的子集之后，在记录下来的所有解中，选一个使目标函数值最优的即可。

穷举法虽然已经可以确保在有限步内获得目标问题的解，但其运算量常常超乎想象（尤其是控制标量和不等式约束条件个数较多时），所以我们在以下两方面寻求改进：

改进最优解的识别方式：利用凸二次规划的特点，建立一组识别规则（实际上就是KKT条件），可直接判断一个可行解是否是最优解，因此算法一旦试出最优解就立即停止，不必遍历剩余可能性
改进尝试各种可能性的顺序：不再随机选取下一次尝试的工作集，而是通过求解子优化问题找到能使目标函数值有效下降的新工作集和迭代点。

经过以上两项改进，就形成了效率大幅提高的有效集法。

有效集法

初始可行解：有效集法需要一个初始可行解作为迭代的起点。由于可行解与目标函数无关，所以获得初始可行解的方法与线性规划完全一致，典型方法有两阶段法和大M法，，因此我们假设已经获得一个初始可行解 x0 ，并得到有效集为W0
迭代求解。迭代的过程主要包括两个部分，一个寻找当前有效集中的最优解，另一个是更新有效集以得到最优集。假设迭代k（k=1,2,3,...）步后，当前解为xk，工作集为wk,根据xk落在可行域上的位置，可以划分出四种情形。

情形1： 若xk 满足KKT条件，则它就是最优解，算法停止。
情形2：若 xk 在当前工作集（等式约束集）下已达最优，即 $s_{k}=x^{*}-x_k=0$ ,但工作集中存在至少一个约束，它原来是不等式约束，而且当 xk 朝该不等式约束的可行一侧移动时，可以使目标函数值进一步下降，(可以通过约束对应的 $\lambda$ 是否为负判断，存在 $\lambda<0$ 说明目标函数值可进一步减小)，那么我们令迭代点不动（即 xk+1=xk并从工作集中去掉那个约束（即将它还原放松成不等式约束）。
情形3：若 xk在当前工作集（等式约束集）下还可以改进，但步长不能走满（否则会走出原始问题的可行域），那么就取 xk 延改进方向恰走到原始可行域边界处的位置为下一个迭代点 xk+1 ，同时把碰到的边界处对应的不等式约束加入工作集。情形3的例子请参见下图中x3
情形4：若 xk在当前工作集（等式约束集）下还可以改进，并且步长能够走满（即改进到最优仍未出原始可行域），那么就取这个最优点作为下一个迭代点 xk+1，并保持工作集不变。情形4的例子请参见下图中 x1和 x4。

下面是使用有效集求解带不等式约束优化问题的经典例子。

G = eye(2,2)*2;
d = [-2 -5]';
A = [-1 2; 1 2; 1 -2; -1 0; 0 -1];
b = [2;6;2;0;0];
Aeq = [-1 1];
beq = [0.03];
MaxIterations = 100;
Tolerance = 1.0e-6;
x0 = [0;1]; workset0= [3;5];
x = active_set_method(x0,G,d, A, b, Aeq, beq, MaxIterations, Tolerance)
%
options = optimoptions(@quadprog,'Algorithm','active-set','MaxIterations',500);
y = quadprog(G, d, A, b, Aeq, beq, [], [], x0, options)

function x = active_set_method(x0,G,d, Ai, bi, Aeq, beq, MaxIterations, Tolerance)
%min 1/2*x'Gx+d'x
%st. A*x <= b
%    Aeq*x = b

        ne = size(Aeq,1);
        ni = size(Ai,1);
        workset = zeros(ni,1);
        x = x0;%给定初始解
        for i = 1:ni
            %将满足初始解的不等式都加入有效集
           if abs(Ai(i,:)*x0 -bi(i)) <= Tolerance
               workset(i) = 1;
           end
        end
        
    for k = 1:MaxIterations
     %迭代求解   
        a = [Aeq; Ai(workset == 1,:)];
        b_ = [beq; bi(workset == 1,:)];
        
        [x1, lambda]=kkt(G, d, a, b_);
        s = x1 - x;
        if norm(s,2) < Tolerance
            %得到的解是本次有效集中的最优解，对应情形1,2
            x = x1; 
            min_lambda = 0;
            if length(lambda) > ne
                [min_lambda,index] = min(lambda(ne+1:end));
            end
           if min_lambda >= 0 %满足KKT条件，得到最优解，对应情形1
               break;
           else
               %寻找新的有效集，对应情形2
              for i=1:ni
                if workset(i)&&sum(workset(1:i))==index
                    workset(i)=0;
                    break;
                end
              end 
           end
        else
            %得到的解不是本次有效集中的最优解，继续优化，对应情形3,4
            [alpha, index] = Alpha(workset, x, s, Ai, bi);
            x = x + alpha*s;
            if alpha <1
                workset(index) = 1;
            end
        end
    end

    function [s, lambda]=kkt(G, g, Ae, be)
        ws = size(Ae,1);
        kkt_A = [G, Ae'; Ae, zeros(ws,ws)];
        kkt_B = [-g;be];
        
        slambda = pinv(kkt_A)*kkt_B; %kkt_A不一定可逆，需要求伪逆
        dim = size(G,2);
        s = slambda(1:dim);
        lambda = slambda(dim+1:end);
    end

    function [alpha, index] = Alpha(workset, x, s, Ai, bi)
        outset = find(workset == 0);
        Atp = Ai(outset,:)*s;
        Atp_index = Atp>0;
        out_plus_ind = outset(Atp_index);
        alphatset = (bi(out_plus_ind) - Ai(out_plus_ind,:)*x)./Atp(Atp_index);
        [alpha, index] = min([alphatset;1]);
        if alpha < 1
            index = out_plus_ind(index);
        end
    end
end

参考：

有效集法 - 知乎

浅谈最优化问题的KKT条件 - 知乎

有效集法介绍（Active Set Method）_在线有效集策略-CSDN博客

Active set method介绍 - 知乎