大M单纯形算法的MATLAB实现

一、理论

我的上一篇博文介绍了最基础的单纯形算法的MATLAB实现，但普通的单纯形算法都是从线性规划的规范式开始迭代，并且要求初始的规范式能显易地给出一个初始基可行解，虽说对于任意一个系数矩阵行满秩的标准型线性规划问题都可以转化为规范式，但如果对于一个每一个非规范式的线性规划问题，都需要手动转化成规范式才能求解，这样的算法似乎有点low了，于是就有了本文将介绍的一种初始基可行解的求解方法，大M单纯形算法。

首先来简述一下算法的思想，对于一个线性规划问题，第一步将其转化为标准型，此时若该标准型的系数矩阵中含有一个单位矩阵，则该标准型是一个规范式，否则引入人工变量，并在改人工变量的价格系数前加上一个充分大的正数M，作为惩罚项，这样可以使得在极小化目标函数的过程中，人工变量逐步转化为非基变量，这样最后就可以排除人工变量的影响，引入人工变量后就可以得到一个转化的标准型线性规划问题，该标准型为规范式，就可以继续用普通标准型算法进行迭代，最终求出最优解。

二、代码

算法具体的思想及一些定理的证明，随便找一本最优化理论的书就有，在此不过多赘述，上代码

function [x,y,ResultFlag]=BigMSimplexAlgorithm(A,B,C,varargin)
% 2020-3-31 臻orz
%inputs:
%   A:系数矩阵 m*n
%   B:右端向量 m*1
%   C:价格系数向量 n*1
%alternative inputs:
%   target:优化目标 0 ~ min; 1 ~ max;
%   sign:约束条件符号向量 m*1 其中
%       -1 ~ '<=';0 ~ '=';1 ~ '>'; 默认为0,即为等式
%   M:惩罚系数，为一充分大正数,默认为10*sum(sum(abs(A)))
%outputs:
%   x:最优解 n*1
%   y:最优值 num
%   ResultFlag:是否找到最优解%check inputs
ip = inputParser;
ip.addRequired('A',@(x)validateattributes(x,{'double'},...{'finite','nonnan'},'BigMSimplexAlgorithm','A',1));
ip.addRequired('B',@(x)validateattributes(x,{'double'},...{'size',[size(A,1),1]},'BigMSimplexAlgorithm','B',2));
ip.addRequired('C',@(x)validateattributes(x,{'double'},...{'size',[size(A,2),1]},'BigMSimplexAlgorithm','C',3));
ip.addParameter('target',0,@(x)validateattributes(x,...{'double'},{'scalar'},'BigMSimplexAlgorithm','target'));
ip.addParameter('sign',zeros(size(A,1),1),@(x)validateattributes(x,...{'double'},{'size',[size(A,1),1]},'BigMSimplexAlgorithm','sign'));
ip.addParameter('M',10*sum(sum(abs(A))),@(x)validateattributes(x,...{'double'},{'scalar','positive'},'BigMSimplexAlgorithm','M'));
ip.parse(A,B,C,varargin{:});%initialize
target = ip.Results.target;
sign = ip.Results.sign;
M = ip.Results.M;
[m,n] = size(A);
P = [];
x = zeros(n,1);
y = 0;
ResultFlag = 0;
j = 0;%standardization
if targetC = -C;%目标函数的转化
end
A(B<0,:) = -A(B<0,:);
sign(B<0,:) = -sign(B<0,:);
B = abs(B);%约束条件的转化
for i = sign'j = j+1;switch icase -1%引入松弛变量a = zeros(m,1);a(j) = 1;A = [A a];C = [C;0];case 1%引入剩余变量a = zeros(m,1);a(j) = -1;A = [A a];C = [C;0];end
end
n1 = size(A,2);%记录转化标准型的未知量的个数%找寻单位矩阵
for i = 1:ma = 0;for j = find(A(i,:)==1)if sum(A(:,j)==0) == m-1P = [P j];a = 1;endendif ~a%若该行无基解,引入人工变量j = zeros(m,1);j(i) = 1;A = [A j];P = [P size(A,2)];C = [C;M];end
endP = P(1:m);
CB = C(P);%基变量对应的价值系数
sigma = C'-CB'*inv(A(:,P))*A;
sigma(P) = 0;while 1if ~sum(sigma<0)if sum(P>n1)%如果基变量含有人工变量return;endx = zeros(n1,1);x(P) = B;x = x(1:n);%舍去引入的松弛变量与剩余变量if targety = -CB'*B;elsey = CB'*B;endResultFlag = 1;return;endpivot_y = find(sigma==min(sigma));pivot_y = pivot_y(1);if sum(A(:,pivot_y)<0) == mreturn;endtheta_index = find(A(:,pivot_y)>0);theta = B(theta_index)./A(theta_index,pivot_y);pivot_x = theta_index(theta==min(theta));%确定主元pivot_x = pivot_x(1);P(pivot_x) = pivot_y;%更新PCB(pivot_x) = C(pivot_y);%更新CB%更新系数矩阵B(pivot_x) = B(pivot_x)/A(pivot_x,pivot_y);A(pivot_x,:) = A(pivot_x,:)./A(pivot_x,pivot_y);a = 1:m;a(pivot_x) = [];for i = aB(i) = B(i)-A(i,pivot_y)*B(pivot_x);A(i,:) = A(i,:)-A(i,pivot_y)*A(pivot_x,:);    endsigma = sigma-sigma(pivot_y)*A(pivot_x,:);%更新sigma
end
end

然后来讲一下输入输出，代码中已经有过注释，调用方式为

>> [x,y,result]=BigMSimplexAlgorithm(A,B,C,'target',0,'sign',D,'M',10);

其中A，B，C分别为系数矩阵，右端向量，价格系数向量，注意B，C为列向量即可，后面的三个均为可选输入
‘target’：为优化的目标为求最小值还是最大值，分别对应0和1，可以不设置，默认为0，即求最小值
‘sign’：为约束条件的符号，小于、等于、大于分别对应-1，0，1，也是一列向量，可以不设置，默认为0，即全为相等
‘M’：为惩罚项，是一个充分大的正数，但也不可取得过大，会导致计算量大倒是事小，主要是因为有可能会使得规划问题变得病态，即矩阵中有些元素过小，有些过大，此时在计算时由于计算机存储的精度问题，导致计算出现丢失现象，可以不设置，默认取值为10*sum(sum(abs(A)))

举个例子

对于线性规划问题

就可以有如下输入

A=[1 -2 1 0;-4 1 2 -1;-2 0 1 0];%系数矩阵
B=[11 3 1]';%右端向量
C=[3 -1 -1 0]';%价格系数向量
D=[-1 0 0]';%符号向量[x,y,result]=BigMSimplexAlgorithm(A,B,C,'target',1,'sign',D,'M',10);
%或者也可以
[x,y,result]=BigMSimplexAlgorithm(A,B,C,'sign',D,'target',1);
%后三个可选参数顺序无所谓，M可以自己设置也可以使用默认
%由于这里要求解目标函数为最大值，所以设置'target'为1
%由于约束条件有小于等于，故'sign'不可以使用默认值，需设置

最后得出的结果为

故原线性规划问题的最优解为(4,1,9,0)^T，最优值为2