线性规划 单纯形法算法 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