線性規劃 單純形法演算法 MATLAB實現

例如： 

問題一：假設產生Q1、Q2、Q3的量分別為X1、X2、X3，那麼

1. 約束條件：
   • P1量<1500，即 2X1+3X2+0X3 ≤1500
   • P2量<800， 即 0X1+2X2+4X3 ≤800
   • P3量<2000，即 3X1+2X2+5X3 ≤2000
   • 產量非負， 即X1,X2,X3 > 0
2. 目標函式最大化  max(3X1+5X2+4X3)

加入鬆弛變數，化為標準形式：

A=[2 3 0 1 0 0;0 2 4 0 1 0;3 2 5 0 0 1];

B=[1500;800;2000];

C=[3 5 4 0 0 0];

然後執行function [ X, Z] = ssimplex( A, b, C ,D)；

% Z: 目標函式的極小值 % A: 約束函式的係數矩陣 % B: 約束函式的常數列向量 % C: 目標函式的係數向量 % D: 求最大值為1，求最小值為0

函式程式碼如下所示：

function [ X, Z] = myssimplex( A, B, C ,D)
% 單純形法的實現
% X: 目標函式的最優解
% Z: 目標函式的極小值
% A: 約束函式的係數矩陣
% B: 約束函式的常數列向量
% C: 目標函式的係數向量
% D: 求最大值為1，求最小值為0
% E: 將B，A拼接成新的矩陣
flag = 1;
S=0;
Z=0;
[m, n] = size(A);
BIndex = n - m  + 1 : n;                % 基向量下標集合
    if D==0                             %如果求最小值將C取反
        C=-C;
    end
    if (n < m)
    disp('係數矩陣不符合要求！')
    else
        while flag
            Cb = C(BIndex);             % 基矩陣對應的目標值b
            Zj = (Cb)*A;
            Rj =C-Zj;                   %判別數
            X = zeros(1, n);
            if max(Rj)<=0
            for i=1:m
                   X(BIndex(i))=B(i);
            end
            for i=1:n
                Z=Z+(C(i)*X(i));
            end
               flag = 0;
               fprintf('迭代次數為:%d\n',S);
               disp('已找到最優解：')
               S=0;
               break;
            else
                S=S+1;
                 [~, k1] = max(Rj); %獲得換入基變數的位置
                 for i=1:m
                     if A(i,k1)>0
                        P(i)=B(i)/A(i,k1);
                     else 
                        P(i)=inf;
                     end   
                 end
                [~, k2]=min(P);    %獲得換出基變數的位置
                BIndex(k2)=k1;     %新的基變數
                E=[B,A];           %B,A合成一個新的矩陣
                E(k2,:)=E(k2,:)/E(k2,k1+1); %將A中k2行，除於A中k2行，k1列；
                for i=1:m                   %更新E
                  if(i==k2)
                      continue;
                  end
                    while(1)
                        E(i,:)= E(i,:)-E(i,k1+1)*E(k2,:);%將E進行初等行換
                        if(E(i,k1+1)==0)
                        break;
                        end
                    end
                    B=E(1:m,1); %E拆開
                    A=E(1:m,2:n+1);
               end
            end
        end
    end
end
 

執行程式結果如下所示：

A=[2 3 0 1 0 0;0 2 4 0 1 0;3 2 5 0 0 1];
B=[1500;800;2000];
C=[3 5 4 0 0 0];
D=1;
[ X, Z] = myssimplex( A, B, C ,D)
迭代次數為:3
已找到最優解：

X =

   375   250    75     0     0     0


Z =

        2675

注：在約束函式的係數矩陣中，單位向量一定放在矩陣的最後，如下所示：

A =

     2     3     0     1     0     0      0     2     4     0     1     0      3     2     5     0     0     1