大M法（Big M Method）

大M法（Big M Method）

前面一篇讲的单纯形方法的实现,但程序输入的必须是已经有初始基本可行解的单纯形表。

但实际问题中很少有现成的基本可行解，比如以下这个问题：

min f(x) = –3x1 +x2 + x3

s.t. x1 – 2x2 + x3 + x4=11

      -4x1 + x2 + 2x3 - x5=3

      -2x1+x3=1

      xj>=0 , j=1,2,3,4,5

写成单纯形表就是

  x1 x2 x3 x4 x5 b f 3 -1 -1 0 0 0   1 -2 1 1 0 11   -4 1 2 0 -1 3   -2 0 1 0 0 1

很难找到秩为3的基阵，更不用说直接出现3阶单位阵了。

在实际问题中，尤其是约束条件变多时，找到基阵甚至是判定A是否满秩都十分困难，因此在程序中引入大M法（Big M Method）来获得初始的基本可行解，这样我们能处理的问题就更加多样化了。

上篇已经说过，对于m*n的矩阵A来说，找到一个m*m 的满秩方阵就能得到基本可行解，但是在这么多列向量中怎样挑出m个线性无关的向量来组成一个满秩方阵呢？如果找起来麻烦的话，不如直接添加一个m阶单位阵来的方便！

大M法

大M法又称惩罚法，它是在目标函数中添加m个人工变量M*x（M是一个任意大的正数），同时在A矩阵中添加一个m阶单位矩阵。

这样一来新的A矩阵中就有了一个m*m满秩方阵，满足单纯形法求解的初始要求，但是若要得到最小值f(x),新添加的人工变量的值必然是0的，因为M可以是很大的数，如果Xn+1不为0，f(x)可能会很大，如果无法做到令人工变量取0值，那么原问题就无可行解。

需要注意的是，添加完人工变量之后，人工变量构成一组可行解的基变量，但单纯形初始矩阵要求基变量对应的检验数为0，故需要做行变换把基变量对应的检验数置0。

例如，本文开始引入的问题经过添加人工变量后变为

  x1 x2 x3 x4 x5 x6 x7 x8 b f 3 -1 -1 0 0 -M -M -M 0 x6 1 -2 1 1 0 1 0 0 11 x7 -4 1 2 0 -1 0 1 0 3 x8 -2 0 1 0 0 0 0 1 1

再进行行变换把基变量x6，x7，x8对应的检验数置0，得到：

  x1 x2 x3 x4 x5 x6 x7 x8 b f 3-5M -1-M -1+4M 0 0 0 0 0 0 x6 1 -2 1 1 0 1 0 0 11 x7 -4 1 2 0 -1 0 1 0 3 x8 -2 0 1 0 0 0 0 1 1

进行完这步之后，就回到了单纯形法求解的基本问题，利用原来的算法继续计算就好了。

Matlab实现

BigM.m

function [ x,y ] = BigM( f,A,b )

%输入f是检验数的数组，1*n维

%输入A是约束矩阵， m*n维

%输入b是约束向量， 1*m维

%输出x是解向量

%输出y是最优解

%判断输入维数是否相符

%做初始单纯形表,加入M变量

[n,m]=size(A);%n行m列

M=10000;

S=[f -1*M*ones(1,n) 0;

   A   eye(n)       b'];

format rat %将结果以分数表示

[n,m]=size(S);

%将人工变量的检验数置零

for k=1:n-1

    S(1,:)=S(1,:)+S(k+1,:)*M;

end

%判断检验数 r<=0

r=find(S(1,1:m-1)>0);

len=length(r);

flag=0;

%有大于0的检验数则进入循环

while(len)

    %检查非负检验数所对列向量元素是否都小于等于0

    for k=1:length(r)

        d=find(S(:,r(k))>0);

        if(length(d)+1==2)

        error('无最优解！！！')

        %break;

        end

    end

    %找到检验数中最大值

    [Rk,j]=max(S(1,1:m-1));

    %最大值所在列比值为正数且最小值br/a_rk

    br=S(2:n,m)./S(2:n,j);

    %把比值中的负数都变无穷

    for p=1:length(br)

        if(br(p)<0)br(p)=Inf;

        end

    end

    [h,i]=min(br);%列向量比值最小值

    % i+1为转轴元行号(在S中)，j为转轴元列号

    i=i+1;

    %进行换基,转轴元置1

    S(i,:)=S(i,:)./S(i,j);

    %转轴元所在列其他元素都置0

    for k=1:n

        if(k~=i)

            S(k,:)=S(k,:)-S(i,:)*S(k,j);

        end

    end

    %判断检验数 r<=0

    r=find(S(1,1:m-1)>0);

    len=length(r);

%     %调试用，控制循环步数

%     if(len>0)flag=flag+1;end

%     if(flag==2)break;end

%     S

end

    %检验数全部非正，找到最优解

    %非基变量置0

    x=zeros(1,m-1);

    for i=1:m-1

        %找到基变量

        j=find(S(:,i)==1);%每列中1的个数

        k=find(S(:,i)==0);%每列中0的个数

        if((length(j)+1==2)&&(length(k)+1==n))

            %i为基本元列号，j是行号

            x(i)=S(j,m);

        end

    end

    y=S(1,m);%最优解

    S

end

测试程序：

f=[3 -1 -1 0 0];

A=[1 -2 1 1  0;

  -4  1 2 0 -1;

  -2  0 1 0  0 ];

b=[11 3 1 ];

[x,y]=BigM(f,A,b)

f=[5 2 3 -1];

A=[1 2 3 0 ;

   2 1 5 0 ;

   1 2 4 1 ];

b=[15 20 26];

[x,y]=BigM(f,A,b)

f=[5 10 0 0 0 ];

A=[1/14 1/7 1 0 0;

   1/7 1/12 0 1 0;

   1    1   0 0 1 ];

b=[1 1 8];

[x,y]=BigM(f,A,b)

[x,y]=Simplex(f,A,b)

计算结果：