本文构建SimplexMax函数,通过构建单纯型表和循环迭代,求解线性规划问题的最优解
clc;clear;
%% 设置变量,调用函数
% 题目参数
A = [0 5 1 0 0;
6 2 0 1 0;
1 1 0 0 1];
b = [15; 24; 5];
c = [2 1 0 0 0];
ind = [3 4 5]
% 无界解
% A = [-2 1 1 0;
% 1 -1 0 1];
% b = [4; 2];
% c = [1 1 0 0];
% ind = [3 4];
% 无穷多最优解
% A = [-1 -1.9 1 0 0 0;
% 1 -1.9 0 1 0 0;
% 1 1.9 0 0 1 0;
% -1 1.9 0 0 0 1];
% b = [-3.8; 3.8; 10.2; 3.8];
% c = [3 5.7 0 0 0 0];
% ind = [3 4 5 6];
[x, z, ST, ca] = SimplexMax(c, A, b, ind) % 调用下述构建的方法进行求解
delete('mysimpleMax.xlsx'); % 删除原有的Excel表(如果存在的话)
fprintf('正在将单纯形表写入Excel...\n');
writematrix('C1','mysimpleMax.xlsx','Sheet',1,'Range','C1');
writecell({'cB','xB', 'b'},'mysimpleMax.xlsx','Sheet',1,'Range','A2');
[~,n] = size(A);
X = strcat('x', string(1:n));
writematrix(X,'mysimpleMax.xlsx','Sheet',1,'Range','D2');
writematrix(c,'mysimpleMax.xlsx','Sheet',1,'Range','D1');
writematrix(ST,'mysimpleMax.xlsx','Sheet',1,'Range','A3');
fprintf('单纯性表写入完成\n');
%% 构造解单纯形法函数
function [x,z,ST,res_case] = SimplexMax(c,A,b,ind_B)
% 单纯形法求解标准形线性规划问题: max cx s.t. Ax=b x>=0
% 输入参数: c为目标函数系数, A为约束方程组系数矩阵, b为约束方程组常数项, ind_B为初始基变量索引
% 输出参数: x最优解, z最优目标函数值, ST存储单纯形表数据,
% res_case=0表示有最优解,res_case=1表示有无界解,res_case=2表示有无穷多解
[m,n] = size(A); % m约束条件个数, n决策变量数
ind_N = setdiff(1:n, ind_B); % 非基变量的索引
ST = [];
format rat
% 循环求解
iteration_round = 0; % 记录迭代次数
while true
x0 = zeros(n,1);
x0(ind_B) = b; % 初始基可行解
cB = c(ind_B); % 计算cB
Sigma = zeros(1,n); % 创建Sigma存储检验数并初始化为0
Sigma(ind_N) = c(ind_N) - cB*A(:,ind_N); % 计算并存入(非基变量)检验数
[~, k] = max(Sigma); % 选出最大检验数, 确定进基变量索引k (注:[a,b]=max(c)中,a为c中的最大值,b为c中最大值的index)
% fprintf('确定进基变量索引:%d\n',k)
Theta = b ./ A(:,k); % 计算θ(取A的k列,使b./(A的k列))
Theta(Theta == 1/0) = NaN;
Theta(Theta <= 0) = NaN;
% Theta(Theta<=0) = 10000;
% [~, q] = min(Theta); % 选出最小θ
[~, q] = min(Theta); % 选出最小正数θ
el = ind_B(q); % 确定出基变量索引el, 主元为A(q,k)
% fprintf('确定出基变量索引:%d\n',el)
vals = [cB',ind_B',b,A,Theta]; % 存储单纯形表
vals = [vals; NaN, NaN, NaN, Sigma, NaN]; % 存储单纯性表Sigma列
[~,vals_width] = size(vals); % 创建一个与vals等宽的空矩阵
vals = [vals; [1:vals_width]* NaN]; % 将该行添加到vals最后作为分隔
ST = [ST; vals]; % 将上述构建的单纯性表vals添加到ST矩阵中
if ~any(Sigma > 0) % 若不存在任何Sigma>0,此基可行解为最优解
fprintf('得到最优解');
x = x0;
z = c * x;
res_case = 0;
return
end
if all(A(:,k) <= 0) %有无界解
fprintf('有无界解');
x = [];
z = NaN;
res_case = 1;
break
end
if ~any(Sigma(ind_B)>0) && any(Sigma(ind_N)==0) % 有无穷多最优解
fprintf('有无穷多最优解')
x = x0;
z = c * x;
res_case = 2;
break
end
iteration_round = iteration_round + 1;
fprintf('迭代第%d次\n',iteration_round)
fprintf('确定进基变量索引:%d 出基变量索引: %d\n',k,el)
% 换基
ind_B(ind_B == el) = k; %新的基变量索引
ind_N = setdiff(1:n, ind_B); %非基变量索引(选取了ind_B和[1:n]的差集)
% 更新A和b
A(:,ind_N) = A(:,ind_B) \ A(:,ind_N); % Q1 = inv(A(:,ind_B))
b = A(:,ind_B) \ b; % b2 = inv(Q1) * b1
A(:,ind_B) = eye(m,m);
end
end
注意:本文采用的对单纯性表进行操作的方法并不是改进单纯型法,如果需要使用改进单纯型法进行求解请参考其他文章,并对下述操作进行改进
% 换基
ind_B(ind_B == el) = k; %新的基变量索引
ind_N = setdiff(1:n, ind_B); %非基变量索引(选取了ind_B和[1:n]的差集)