前面一篇讲的单纯形方法的实现,但程序输入的必须是已经有初始基本可行解的单纯形表。
但实际问题中很少有现成的基本可行解，比如以下这个问题：
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)
计算结果：
