求LP問題的基本(可行)解Matlab代碼
基本原理
Matlab代碼
初始化
是否隨機生成LP問題還是自己指定LP問題
初始化係數矩陣A和b
找到初始的基本可行解
展示線性規劃問題
求解線性規劃的單純型法
function script_LP()
考慮如下線性規劃問題的標準型: max Z= CX s.t. {AX=b,X>=0}
假定A是m*n階矩陣,並且A的秩r(A)=m。那麼,求解該問題的最優解的的步驟:
找到初始的基本可行解;
判斷是否最優;
非最優時,從一個基本可行解變換到更好的基本可行解(進基、出基、旋轉);
單純形法的本質上是每次去找更好的頂點(基本可行解),一直找到最優的頂點(基本可行解)為止。
換句話說,單純形法的要點包括以下三個部分:
第一部分,可以用之前的求解基本可行解的代碼來實現。第二和第三部分,則需要用到LP問題的典式的概念。
clc
clear all
close all
format rat
choose = 0;
if choose
m = 3;
n = 8;
is_opt_flag = 0;
NumMax = 10;
while is_opt_flag==0
A = randi(NumMax,m,n);
if rank(A) ==m
is_opt_flag=1;
end
end
b = randi(NumMax,m,1);
C = randi(NumMax,n,1);
[m,n] = size(A);
idxs = nchoosek(1:n,m);
count = size(idxs,1);
have_no_basic_feasible = 1;
k = 1;
while k <= count && have_no_basic_feasible
BIdx = idxs(k,:);
B = A(:,BIdx);
str = [repmat('P%d,',1,m-1) 'P%d'];
if det(B) ~= 0
X = zeros(1,n);
XB = B\b;
XN = 0;
X(BIdx)=XB;
if sum(XB<0) == 0
fprintf(['\n=========================================\n1. 列向量' str '能構成基陣.\n'],BIdx);
fprintf('2. 該基陣對應的初始基本可行解:\n X= [%s].\n\n=========================================\n',rats(X));
have_no_basic_feasible = 0;
end
end
k=k+1;
end
if have_no_basic_feasible
fprintf('\n=========================================\n1. 未找到初始的基本可行解.\n');
return;
end
else
A = [1 2 1 0 0;3 2 0 1 0; 0 2 0 0 1];
b = [30;60;24];
C = [40;50;0;0;0];
[m,n] = size(A);
X = [0;0;30;60;24];
BIdx = [3,4,5];
end
disp('係數矩陣A是:')
disp(A)
disp('列向量b是:')
disp(b)
disp('係數向量c是:')
disp(C)
fprintf('\n=========================================\n')
disp('目標函數f=CX是:')
str = ['f=' repmat('%dx%d+',1,n-1) '%dx%d.\n'];
v =[];
for j =1:n
v =[v,C(j),j];
end
fprintf(str,v)
disp('約束組AX=b是:')
str = [repmat('%dx%d+',1,n-1) '%dx%d=%d.\n'];
for i = 1:m
v =[];
for j =1:n
v =[v, A(i,j),j];
end
fprintf(str,[v,b(i)])
end
fprintf('\n=========================================\n')
[X,Z] = simplex(A,b,C,X,BIdx);
係數矩陣A是:
1 2 1 0 0
3 2 0 1 0
0 2 0 0 1
列向量b是:
30
60
24
係數向量c是:
40
50
0
0
0
=========================================
目標函數f=CX是:
f=40x1+50x2+0x3+0x4+0x5.
約束組AX=b是:
1x1+2x2+1x3+0x4+0x5=30.
3x1+2x2+0x3+1x4+0x5=60.
0x1+2x2+0x3+0x4+1x5=24.
=========================================
function [X,Z] = simplex(A,b,C,X,BIdx)
單純形法的基本原理,可以從不同角度去理解:
1)從基本解的角度來理解,就是從一個基本可行解出發,去找另一個更好的基本可行解;
2)從典式的角度來理解,就是從一個典式變換到另外一個典式;
3)從矩陣運算的角度來理解,就是找初始基本可行解,判斷最優(檢驗數小於0),換基運算(進基、出基、旋轉運算)的過程。
[m,n] = size(A);
is_opt_flag = 0;
NIdx = setdiff(1:n,BIdx);
iter_num = 0;
while ~is_opt_flag
B = A(:,BIdx);
N = A(:,NIdx);
CB = C(BIdx);
CN = C(NIdx);
X = zeros(1,n);
XB = B\b;
XN = 0;
X(BIdx)=XB;
Z = CB'*XB;
iter_num = iter_num+1;
fprintf('第%d次單純形迭代,當前的基本可行解:\n X= [%s].\n',iter_num,rats(X));
sigma = CN'-CB'/B*N;
[max_sigma,enter_basic_idx] = max(sigma);
if max_sigma <=0
is_opt_flag = 1;
fprintf('\n=========================================\n找到該問題的最優解:\n X= [%s].\n',rats(X));
else
enter_basic_idx = NIdx(enter_basic_idx);
dk = (B\b)./(B\A(:,enter_basic_idx));
if sum(dk<=0) == m
fprintf('3. 該問題無有界解.\n');
return;
else
[d_idx,~] = find(dk>0);
[~, leave_basic_idx] = min(dk(dk>0));
leave_basic_idx = d_idx(leave_basic_idx);
leave_basic_idx = BIdx(leave_basic_idx);
BIdx = union(setdiff(BIdx,leave_basic_idx),enter_basic_idx);
NIdx = union(setdiff(NIdx,enter_basic_idx),leave_basic_idx);
end
end
end
第2次單純形迭代,當前的基本可行解:
X= [0 0 30 60 24].
第2次單純形迭代,當前的基本可行解:
X= [0 12 6 36 0].
第3次單純形迭代,當前的基本可行解:
X= [6 12 0 18 0].
第4次單純形迭代,當前的基本可行解:
X= [15 15/2 0 0 9].
=========================================
找到該問題的最優解:
X= [15 15/2 0 0 9].
Published with MATLAB® R2015b