MATLAB 單純形法演算法
阿新 • • 發佈:2018-12-18
線性規劃 單純形法演算法 MATLAB實現
例如:
問題一:假設產生Q1、Q2、Q3的量分別為X1、X2、X3,那麼
- 約束條件:
- P1量<1500,即 2X1+3X2+0X3 ≤1500
- P2量<800, 即 0X1+2X2+4X3 ≤800
- P3量<2000,即 3X1+2X2+5X3 ≤2000
- 產量非負, 即X1,X2,X3 > 0
- 目標函式最大化 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