【最优化笔记3】线性规划--求解方法（单纯形法及Matlab实现）

单纯形法是求解线性规划问题的一种通用的有效算法（必考点）。
此篇博客以我的好友Cizeron总结为基础完成，特此表示感谢。

1.前置概念

（1）约束方程的规范形式：
$$\{\begin{array}{ll}minf=c^{T}x& \\ s.t.Ax=b,& \text{x≥0 }\end{array}\text{\hspace{1em}(1)}$$
线性规划（1）的约束条件系数矩阵A通过初等行变换，总可以化为
$$[I_m,N]$$
则约束条件可以写为
$$[I_m,N]\left[\begin{array}{c}{x}_B\\ x_N\end{array}\right]=b^{\prime }\text{\hspace{1em}(2)}$$
其中，显然$B=I_m$即为基，$x=(x_B=b^{\prime },0{)}^T$为(1)的关于B基本解。将线性规划问题化为此规范形式是下文编程实现的第一步。
(2)基变换：
若初始基本可行解$x^{(0)}$不是最优解，那么就还要找一个新的基本可行解：从原可行基中换出一个列向量（离基变量），再换入一个新的列向量(进基变量)（要保证线性无关）,从而得到一个新的可行基，这就是基变换。

2.基本思想

单纯形法的基本思想是：给出一种规则，使由 LP问题一个基本可行解（极点）转移到另一个基本可行解，目标函数值是减小的，而且两个基本可行解之间的转换是容易实现的，经过有限次迭代，即可求得所需的最优基本可行解。

再具体一点来说：
对于一个优化问题$$\{\begin{array}{ll}minf=c^{T}x& \\ s.t.Ax=b,& \text{x≥0 }\end{array}\text{\hspace{1em}(1)}$$ 这里$b\ge 0,rank(A)=m$,记$A=[B,N]=\left[\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1n}\\ a_{21}& a_{22}& ...& a_{2n}\\ ...& ...& ...& ...\\ a_{m1}& a_{m2}& ...& a_{mn}\end{array}\right]=(p_1,p_2,...,p_n)$,其中$rank(B)=m$，即B为初始可行基。
设$x^{(0)}=\left[\begin{array}{c}{B}^{-1}b\\ 0\end{array}\right]$为问题（1）的一个初始可行解，则在$x^{(0)}$处的目标函数值
$f_0=c^T{x}^{(0)}=[c_B,c_N]\left[\begin{array}{c}{B}^{-1}b\\ 0\end{array}\right]=c_B{B}^{-1}b$。
若初始可行解$x^{(0)}$不是最优解，则需要找到一个新的基本可行解，即进行基变换。
设$x=\left[\begin{array}{c}{x}_B\\ x_N\end{array}\right]$为任意一个可行解，则由$Ax=b$得：$x_B=B^{-1}b-B^{-1}N{x}_N$,在该$x$点的目标函数值为：$f=c^{T}x=c_B{x}_B+c_N{x}_N=c_B(B^{-1}b-B^{-1}N{x}_N)+c_N{x}_N=c_B{B}^{-1}b-(c_B{B}^{-1}N-c_N)x_N=c_B(B^{-1}b-B^{-1}N{x}_N)+c_N{x}_N=c_B{B}^{-1}b-{\sum }_{j\in R_N}(c_B{B}^{-1}{p}_j-c_j)x_j=f_0-{\sum }_{j\in R_N}(z_j-c_j)x_j$，其中$R_N$为非基变量的下标集以及$z=c_B{B}^{-1}{p}_j$。最后的可行解函数值为：$$f=f_0-\sum _{j\in R_N}(z_j-c_j)x_j$$因而，适当选取非基变量$x_j(j\in R_N)$的数值就有可能使得${\sum }_{j\in R_N}(z_j-c_j)x_j>0$从而达到目标函数下降的目的。为此，我们在原来的$n-m$个非基变量中，令$n-m-1$个变量仍取0值，而那个“天选之子”（比如$x_k$）取正值，只要$(z_k-c_k)>0$就能达到优化的目的。
至于这个$k$怎么选，怎么生成新的基本可行解，可以看看上一节的线性规划基本定理的证明。下面要具体说说怎么算。

3.算法步骤

Step1： （当然要先化为标准型）然后找到初始可行基B，解$B{x}_B=b$,求得$x_B=B^{-1}b=\tilde{b}$,令$x_N=0$,确定初始可行解。并计算初始目标函数值$f=c_B{x}_B，$其中$c_B$代表$x_B$对应的目标函数系数向量。

Step2： 求单纯形乘子$w=c_B{B}^{-1}$。对于所有非基变量，计算判别（检验)数$(z_j-c_j)=w{p}_j-c_j$ 。令$z_k-c_k=ma{x}_{j\in R_N}${$z_j-c_j$}，$（R_N$是非基变量的下标集）若$z_k-c_k\le 0$，则对于所有非基变量$z_j-c_j\le 0$,而对于基变量的判别数总是零，因此停止计算，现行基本可行解是最优解。否则，进行下一步。

Step3： 解$B_k{y}_k=p_k$，得到$y_k=B_k^{-1}{p}_k$,若$y_k\le 0$,即$y_k$的每一个分量均为非正数，则停止计算：问题不存在有限最优解。否则，进行第四步。

Step4： 确定下标$r：\frac{\widetilde{b_r}}{y_{rk}}=min${$\frac{\widetilde{b_i}}{y_{ik}}|y_{ik}>0$}。$x_{Br}$为离基变量，$x_k$为进基变量。用$p_k$替换$p_{Br}$,得到新的矩阵B，返回第一步。

考试计算中，我们经常使用单纯形表完成上述计算。
不懂？
看栗子就完事了！

4.算例

ONE: 化为标准型，并找到初始可行基B确定对应的基本可行解$x_B$
标准型如下：

从标准型容易得到一个以$B=I_2$为基底的基本可行解：$x_B=B^{-1}b=[80,90{]}^T,x=[x_N,x_B]=[0,0,80,90{]}^T,c_B=[0,0],$非基变量的下标集$R_N=${$1,2$}。

TWO：计算初始检验数$q_j$及目标函数值$z_0$,建立初始单纯形表
计算检验数$q_j=z_j-c_j=w{p}_j-c_j=c_B{B}^{-1}{p}_j-c_j$和目标函数值$z_0=f=c_B^T{x}_B=c_B{B}^{-1}b$,建立如下图所示的单纯形表：

具体对于本题来说这里$w=c_B{B}^{-1}=0$，因此初始检验数为$q_j=-c_j$，目标函数值为$z_0=f_0=c_B^T{x}_B=0,$因此初始单纯形表如下：

THREE：决定进基矢量$a_k$
取最大检验数$q_2=16>0$所对应的变量作为进基矢（变）量，即$k=2$。

FOUR：决定离基矢量$a_r$和主元素$y_{rk}$
$y_{ij},j\ne 0$即单纯形表第i行的第j个变量对应列的对应元素，$y_{i0}$即单纯形表最后一列($B^{-1}b$)的第i行元素。
首先计算比值$y_{i0}/y_{i}k$,这里$k=2,i=1,2$。结果为：
$\frac{y_{10}}{y_{12}}=80/4=20；$$\frac{y_{20}}{y_{22}}=90/3=30$
故由$r：\frac{\widetilde{b_r}}{y_{rk}}=min${$\frac{\widetilde{b_i}}{y_{ik}}|y_{ik}>0$}。得$r=1$。因此主元素为$y_{12}$，离基变量为$x_3$。

FIVE：以$y_{rk}$为主元素更新（进行Gauss消元）单纯形表
对单纯形表进行初等行变换使得主元素$y_{rk}$所在位置变为1，且对应列的其他行都变为0。
对于本题将第一行除以4，然后分别乘以(-3)和(-16)加到第2、3行，完成第一次迭代，函数值减小为为-320。

然后回到THREE进行下次迭代。
那么什么时候停止呢？下面来说

5.算法收敛性

以极小化问题为例，最大检验数$q_r=ma{x}_{j\in R^N}(z_j-c_j)=ma{x}_{j\in R^N}(w{p}_j-c_j)=ma{x}_{j\in R^N}(c_B{B}^{-1}{p}_j-c_j)$每次迭代必然出现下面三种情形：
(1)$q_r\le 0$：这时现行基本可行解就是最优解。
(2)$q_r>0,y_{rk}\le 0$：当$x_k$无线增大时，目标函数趋于负无穷，因此解无解。
(3)$q_r>0,y_{rk}>0$：这时求出的新的基本可行解，经迭代使目标函数下降。
收敛定理： 对于非退化问题（即存在非退化基本可行解），单纯形方法经有限次迭代后达到最优基本可行解，或得出无界的结论。
所以上面的那个例子在迭代了2次后最大检验数$q_r\le 0$,迭代停止，此时的解即为最优解。

6.Matlab实现

这里实现已经化为标准型且约束方程是规范形式。其实后面学习我们知道其他线性优化问题均可以通过M法或两阶段法转化至此。
通过上述算例的步骤一步一步实现即可，不算困难。大致可分四步：1.输入问题；2.建立初始单纯形表；3.迭代寻找最优解；4.输出结果。
因为这次期末时间过于紧张，这里先把代码总表给出，之后再按下面的步骤走一遍。

1.输入问题

2.建立初始单纯形表

3.迭代寻找最优解

4.输出结果

5.附录（代码总表）

% 输入一标准型，实现程序求求解
% 若有m个方程，则初始基本可行解为E，且位于最后mxm矩阵
% matlab代码
% Made by Fei

% 返回最优解（变量名和对应取值）和最优函数值
%c,A,b分别为目标函数系数向量，约束函数系数矩阵，约束函数右端常数

function [Y,T] = Mysolve(c,A,b)

% 预备工作
% 求出变量个数以及方程个数
n_numx = size(A,2);
n_eq = rank(A);
% 记录初始基本可行解及相应的变量所在位置
XB = b;
XB_index=[];
for i=1:n_eq
    XB_index(i)=n_numx-n_eq+i;
end
% 求出单纯形乘子
cB = c(:,(n_numx-n_eq+1):n_numx);
w = cB;
% 求出初始函数值
f0 = cB*XB;
% 求出n_numx个检验数
z=[];
for i=1:n_numx
    if(i>n_numx-n_eq)
        z(i)=0;
    else
        temp = w*A(:,(n_numx-n_eq-i+1))-c(:,i);
        z(i)=temp;
    end
end

% 开始迭代求解
flag=0;
while(flag==0)
    % 判断检验数是否全部全部为非正数，如是则停止迭代平输出结果
    flag=1;
    for i=1:n_numx
        temp=z(i);
        if(temp>0)
            flag=0;
        end
    end
    % 如果全负则输出结果
    if(flag==1)
        Y=[XB_index;XB']; % 变量名及所对应的取值 
        T=f0; % 最优目标函数值
        break
    end
    
    %求解过程
    % 找出最大检验数及并确定进基变量
    Max=z(1);
    index_in=1;
     for i=2:n_numx
        if(z(i)>Max)
            Max=z(i);
            index_in=i;
        end
     end
     % 找出yr以及出基变量
     y=b./A(:,index_in);
     index_out_H=0;
     for i=1:n_eq
         if(y(i)<=0)
            continue;
         else
             min=y(i);
             % 分别记下对应的行以及出基变量
             index_out_H=i;
             index_out=XB_index(i);
             break;
         end
     end
     % 如果yi均为负数则无最优解，结束迭代
     % 否则寻找yr及出基变量
     if(index_out_H==0)
         fprintf("该问题无最优值解");
         break;
     else
         for i=index_out_H:n_eq
             if(y(i)<min)
                 min=y(i);
                 index_out_H=i;
                 index_out=XB_index(i);
             end
         end
     end
     % 到这里我们已经找到了进基变量与出基变量的位置，及相应的y和检验数z
     
     % 接下来我们尝试进行转轴运算更新单纯形表
     b(index_out_H)=b(index_out_H)./A(index_out_H,index_in);
     A(index_out_H,:)=A(index_out_H,:)./A(index_out_H,index_in);
     % 更新A和B-1b
     for i=1:n_eq
         if(i==index_out_H)
             continue;
         else
             temp=A(i,index_in);
             A(i,:)=A(i,:)-temp.*A(index_out_H,:);
             b(i)=b(i)-temp.*b(index_out_H);
         end
     end
     % 更新检验数z和目标函数值f0
     temp=z(index_in);
     for i=1:n_numx
         z(i)=z(i)-temp*A(index_out_H,i);
     end
     f0=f0-temp*b(index_out_H);
     % 更新XB（基本可行解）以及对应的取值
     for i=1:n_eq
         if(XB_index(i)==index_out)
             XB_index(i)=index_in;
         else
             continue;
         end
     end
     XB=b;
end

以上文的例子输入演示结果：

本节完。