先来介绍一下单纯形法，下面解释是从国科大算法最优化课程林姝老师的课件中截取的。
接下来写代码，单纯形法函数：
%% SimplexMax.m function [x, c, z, pt, ind_B, ind_N] = SimplexMax(c, A, b, ind_B, iter_tag) % 单纯形法求解标准形线性规划问题: max cx s.t. Ax=b x>=0 % 输入参数: c为目标函数系数, A为等式约束方程组系数矩阵, b为等式约束方程组常数项, ind_B为松弛变量索引 % 输出参数: % x: 最优解, x中前len(ind_b)个是x_b, 后面的是x_N; % c: 最优系数, c中从前往后分别是x1的系数、x2的系数、x3的系数......; % z: 最优目标函数值 z_0； % ST存储单纯形表数据； % res_case=0表示有最优解，res_case=1表示有无界解 [m,n] = size(A); %m约束条件个数, n决策变量数 ind_N = setdiff(1:n, ind_B); %非松弛变量的索引 ST = []; format rat step_iter = 0; % 循环求解 while true x0 = zeros(n,1); x0(ind_B) = b; %初始基可行解 cB = c(ind_B); %计算cB Sigma = zeros(1,n); Sigma(ind_N) = c(ind_N) - cB*A(:,ind_N); % 计算检验数 [~, k] = max(Sigma); % 选出最大检验数, 确定松弛变量索引k Theta = b ./ A(:,k); % 计算θ Theta(Theta<=0) = 10000; [~, q] = min(Theta); % 选出最小θ el = ind_B(q); % 确定出松弛变量索引el, 主元为A(q,k) vals = [cB', ind_B', b, A, Theta]; vals = [vals; NaN, NaN, NaN, Sigma, NaN]; ST = [ST; vals]; if ~any(Sigma > 0) %此基可行解为最优解, any存在某个>0 disp('当前基可行解为最优解'); x = x0; z = c * x; pt = Sigma(ind_N); res_case = 0; return else disp('当前基可行解不为最优解'); x = x0; z = c * x; pt = Sigma(ind_N); res_case = 0; end % if ~any(Sigma > 0) if all(A(:,k) <= 0) %有无界解 disp('有无界解'); x = []; res_case = 1; break end % if all(A(:,k) <= 0) step_iter = step_iter + 1; if step_iter == iter_tag break end % if step_iter % 换基 ind_B(ind_B == el) = k; % 新的松弛变量索引 ind_N = setdiff(1:n, ind_B); % 自变量索引 % 更新A和b A(:,ind_N) = A(:,ind_B) \ A(:,ind_N); b = A(:,ind_B) \ b; A(:,ind_B) = eye(m,m); end end % function [x, c, z,ST,res_case] = SimplexMax(c,A,b,ind_B)
输入函数，调用单纯形法，然后输出：
%% 单纯形法求线性规划 A = [4 1 1 0; -1 1 0 1]; % 等式约束方程组（包括自变量与松弛变量）系数矩阵 b = [16; 6]; % 等式约束方程组常数项（等号右边的部分） c = [2 3 0 0]; % 目标函数的系数项（包括自变量与松弛变量） ind_B = [3 4]; % 松弛变量的变量索引，一般松弛变量设置地自变量更大，并紧截止自变量定义 iter_tag = 3; % 迭代次数 从1开始，n表示单纯形法运算了几轮 % SimplexMax求的是max，输出改为min的输出就行了 [x, c, z, pt, ind_B, ind_N] = SimplexMax(c, A, b, ind_B, iter_tag); %% 求min时把c、z、pt转负数 c = -c; z = -z; for i=1:length(pt) pt(i)=-pt(i); end %% 输出 disp('B='); disp(A(:,ind_B)); disp('N='); disp(A(:,ind_N)); disp('x_B='); disp(x(ind_B)); disp('x_B的x索引='); disp(ind_B); disp('x_N='); disp(x(ind_N)); disp('x_N的x索引='); disp(ind_N); disp('c_B的转置='); disp(c(ind_B)); disp('c_N的转置='); disp(c(ind_N)); disp('z_0='); disp(z); disp('p的转置 pT='); disp(pt);
输出示例：
当前基可行解不为最优解 当前基可行解不为最优解 当前基可行解为最优解 B= 4 1 -1 1 N= 1 0 0 1 x_B= 2 8 x_B的x索引= 1 2 x_N= 0 0 x_N的x索引= 3 4 c_B的转置= -2 -3 c_N的转置= 0 0 z_0= -28 p的转置 pT= 1 2
注意出现了三行文字：
当前基可行解不为最优解
当前基可行解不为最优解
当前基可行解为最优解
这个意思是说，前两次迭代都不是最优解，可以通过控制iter_tag 来设定最大迭代次数，最后一次是最优解。
