线性规划 Part 1：修正单纯形法

单纯形是 N 维空间中的 N+1 个顶点的凸包，是一个多胞体：直线上的一个线段，平面上的一个三角形，三维空间中的一个四面体等等，都是单纯形。而线性规划的可行域正是由多个超平面围成的单纯形。其最优解一定是在单纯形的顶点上，因此单纯形法就是从一个顶点出发，走向另一个更优的顶点。

1. 问题描述
2. 原理推导
3. 算法流程
4. 实验
5. 参考

1. 问题描述

线性规划的原始问题如下：

\[\begin{aligned} \min_{x\in \mathbb{R}^n} &\quad c^{\top} x \\ \text { s.t. } &\quad Ax=b \\ &\quad x \ge 0 \end{aligned}\]

其中 $A \in \mathbb{R}^{m \times n}$，且行满秩 $\mathrm{Rank}(A)=m$.

2. 原理推导

首先我们可以思考，$A$ 虽然是行满秩，但并不是列满秩，我们可以取出其中的 $m$ 列来构成 $A$ 的一组基. 假设这些列的指标为 $B=\{b_1,b_2,\ldots,b_m\}$，非基指标记为 $N=\{1,2,\ldots,n\}-B$.

那么有

\[Ax=A_Bx_B+A_Nx_N=b\]

由于 $A_B$ 可逆，于是可通过 $x_N$ 唯一确定 $x_B$

\[x_B=A_B^{-1}(b-A_Nx_N)\]

对于每个固定的 $B$，都有唯一 $\hat{x}$ 使得

\[A_B\hat{x}_B=b,\quad \hat{x}_N=0\]

$x$ 与 $\hat{x}$ 的关系为:

\[\begin{aligned} x_B&=A_B^{-1}(b-A_Nx_N)\\ &=A_B^{-1}(A_B\hat{x}_B-A_Nx_N)\\ &=\hat{x}_B-A_B^{-1}A_Nx_N \end{aligned}\]

那么优化函数可以写成：

\[\begin{aligned} c^{\top}x&=c^{\top}_Bx_B+c^{\top}_Nx_N\\ &=c^{\top}_B\hat{x}_B+(c^{\top}_N-c^{\top}_BA_B^{-1}A_N)x_N\\ & = c^{\top}_B\hat{x}_B+r^{\top}x_N \end{aligned}\]

这里 $r=c_N-A_N^{\top}A_B^{-\top}c_B$. 显然，当 $r\ge0$ 时，$\hat{x}$ 即为该线性规划的最优解.

如果 $r$ 中存在分量 $r_s<0$（注意这里的 $s$ 是相对于 $N$ 而言的，我们设其在整体的指标中为 $j$）.

那么我们就可以进行转轴操作(Pivot operations). 转轴操作是单纯形法中的核心操作，其作用是将一个基变量与一个非基变量进行互换。可以将转轴操作理解为从单纯形上的一个顶点走向另一个顶点。

我们将 $j$ 提到基指标中来，称为入基(enter the basis).

一个指标入基了，那么就要有一个指标出基. 对于当前可行解 $\hat{x}$, 我们假设有一个可行方向 $d$, 满足

\[Ad=A_Bd_B+A_Nd_N=0\]

由于 $j$ 入基了，所以 $d_N=e_s\in \mathbb{R}^{n-m}$. 于是可得:

\[d_B=-A_B^{-1}A_Ne_s=-A_B^{-1}a_j\]

这里的 $a_j$ 为$A$ 的第 $j$ 列.

我们设步长为 $t$，于是新的迭代点为 $\hat{x}+td$，新的目标函数为：

\[\begin{aligned} c^{\top}(\hat{x}+td)&=c^{\top}\hat{x}+t(c^{\top}_N-c^{\top}_BA_B^{-1}A_N)e_s\\ &=c^{\top}\hat{x}+tr_s\\ &< c^{\top}\hat{x} \end{aligned}\]

可见目标函数是严格下降的，而且 $t$ 越大，目标函数下降得越多. 但是当 $t$ 过大时，$\hat{x}+td$ 可能会超出可行区域，所以我们要找到最早抵达边界处的分量 $i\in B$

\[\left |\frac{\hat{x}_i}{d_i} \right | = \min_{d_k<0}\left |\frac{\hat{x}_k}{d_k} \right |\]

于是就可以把 $i$ 出基，$j$ 入基:

\[B'=B-\{i\}+\{j\}\]

然后重复以上步骤即可.

3. 算法流程

对于 $B$，计算 $\hat{x}_B=A_B^{-1}b$.
计算 $r=c_N-A_N^{\top}A_B^{-\top}c_B$.
对于 $r$:
- 如果 $r\ge 0$, $\hat{x}$ 即为最优解.
- 否则选出最小的 $r_s<0$，跳转第4步.
计算更新方向 $d_B=-A_B^{-1}a_j$.
如果 $d_B\ge 0$ 则无下界.
寻找 $d_B<0$ 中使得 ${\hat{x}_k}/{d_k}$ 最大的基指标 $i$.
更新 $B' \leftarrow B-\{i\}+\{j\}$.
跳转第1步.

4. 实验

\[\begin{aligned} \min_{x} &\quad -3 x_{1}-x_{2}-3 x_{3} \\ \text { s.t. } &\quad 2 x_{1}+x_{2}+x_{3} \leq 2 \\ &\quad x_{1}+2 x_{2}+3 x_{3} \leq 5 \\ &\quad 2 x_{1}+2 x_{2}+x_{3} \leq 6 \\ &\quad x_{\nu} \geq 0 \quad(1 \leq \nu \leq 3) \end{aligned}\]

这个问题并不是标准的等式约束，但是我们可以通过引入松弛变量(slack variables) $x_{4}, x_{5}, x_{6} \geq 0$, 使之变成：

\[\begin{aligned} \min_{x} &\quad -3 x_{1}-x_{2}-3 x_{3} \\ \text { s.t. } &\quad 2 x_{1}+x_{2}+x_{3}+x_{4}=2 \\ &\quad x_{1}+2 x_{2}+3 x_{3} +x_5= 5 \\ &\quad 2 x_{1}+2 x_{2}+x_{3} +x_{6}=6 \\ &\quad x_{\nu} \geq 0 \quad(1 \leq \nu \leq 6)\\ \end{aligned}\]

那么将以下数据输入程序:

\[A=\left(\begin{array}{llllll}2 & 1 & 1 & 1 & 0 & 0 \\ 1 & 2 & 3 & 0 & 1 & 0 \\ 2 & 2 & 1 & 0 & 0 & 1\end{array}\right), \quad b=\left(\begin{array}{l}2 \\ 5 \\ 6\end{array}\right), \quad c=(-3,-1,-3,0,0,0)^{T}\]

可行基(feasible basis) 为 $(4,5,6)$.

clc,clear

A = [2 1 1 1 0 0;1 2 3 0 1 0;2 2 1 0 0 1];
b = [2;5;6];
c = [-3 -1 -3 0 0 0]';
ind = [4,5,6];
[x, z] = LP_simplex(c, A, b, ind)


function [x, z] = LP_simplex(c,A,b,B)
% 单纯形法求解标准形线性规划问题: 
% max cx s.t. Ax=b x>=0
% 输入参数: c为目标函数系数, 
% A为约束方程组系数矩阵, 
% b为约束方程组常数项, 
% B为基变量索引
[m,n] = size(A);    %m约束条件个数, n决策变量数
N = setdiff(1:n, B);  %非基变量的索引
format rat
% 循环求解
while true
    x0 = zeros(n,1);
    A_B = A(:,B);
    A_N = A(:,N);
    A_Binv = inv(A_B);
    x0(B) = A_Binv*b;   %初始基可行解
    r = c(N) - A_N'*A_Binv'*c(B);   %计算检验数
    if ~any(r < 0)           %此基可行解为最优解, any存在某个>0        
        x = x0;
        z = c' * x;
        fprintf("Find opt solution!")
        return
    end
    [~, s] = min(r);         %进基变量索引 s
    j = N(s);
    d_B = -A_Binv*A(:,j);
    if ~any(d_B<0)
        fprintf("No opt solution!")
        x = x0;
        z = c * x;
        return
    end
    Theta = x0(B) ./ d_B;         %计算θ
    Theta(Theta>=0) = -1e+9;
    [~, q] = max(Theta);
    i = B(q);        
    % 换基
    B(B == i) = j;      %新的基变量索引
    N = setdiff(1:n, B); %非基变量索引
end
end

输出为:

Find opt solution!
x =

       1/5     
       0       
       8/5     
       0       
       0       
       4       


z =

     -27/5     

注意，如果这里的 $b_1=-2$ ，那么我们的初始解 $x_4=-2<0$，是不在可行域内的. 这时可以先将 $A$ 的第一行和 $b_1$ 同时乘一个 $-1$.

5. 参考

W. Forst and D. Hoffmann, Optimization—Theory and Practice. New York, NY: Springer New York, 2010. https://doi.org/10.1007/978-0-387-78977-4.
单纯形法求解线性规划问题的Matlab实现 - 知乎
简单理解线性规划的单纯形算法 - 知乎
数值优化笔记整理（A）——线性规划中的单纯形法与内点法 - 知乎