前面一篇講的單純形方法的實現,但程序輸入的必須是已經有初始基本可行解的單純形表。
但實際問題中很少有現成的基本可行解，比如以下這個問題：
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)
計算結果：
