用Matlab编写的单纯形法（两阶段）

• 因为这个学期学习了《运筹学》这门课，所以尝试编写了里面介绍的第一类算法
• 用于解决LP（linear program）问题的单纯形法
  这是我的GitHub仓库，你可以去fork
• 里面整体没用问题，但是运行显示还是有问题，看有空的时候再改吧！

下面就贴上代码

建议先看一下README.md文件

1. seekUnitMatrix.m文件（用来寻找单位矩阵）

function [ indexUnitMatrix ] = seekUnitMatrix( inputMatrix )
%seekUnitMatrix 从矩阵（可以不是方阵）中找到可以构成单位矩阵的列
%parameter:inputMatrix---->输入矩阵
%return:indexUnitMatrix--->可以顺序构成单位矩阵的列，若构不成单位矩阵，则返回0

indexUnitMatrix = [];
[row_inputMatrix,colomn_inputMatrix] = size(inputMatrix);
if row_inputMatrix <= colomn_inputMatrix
    for i = 1:row_inputMatrix
        E = zeros(row_inputMatrix,1);
        E(i) = 1;
        for j = 1:colomn_inputMatrix
            if all(inputMatrix(:,j)==E)
                indexUnitMatrix = [indexUnitMatrix,j];
                break; % 可能存在多列是相同的，只要一列的索引就行
            end
        end
    end
    
elseif row_inputMatrix > colomn_inputMatrix
    inputMatrix = inputMatrix';
    indexUnitMatrix = seekUnitMatrix( inputMatrix ); % 递归调用

end
if isempty(indexUnitMatrix)
    indexUnitMatrix = 0;
end
end

2. simplexMethod.m文件（普通的单纯形法，只是还不能判断无解的情况）

function [S,V,rA,rb,rXB] = simplexMethod(C,A,b)
%simplexMethod 单纯形法解线性规划,利用单纯形表，基本思路：
%   step1：找出一初始基及基可行解（一般都是找单位矩阵）
%   step2：判断该可行解是否为最优解，如果是，则停止运算；否就转step3
%   step3：换入换出转换基，转step2
%现在似乎只能解增添的松弛变量构成单位矩阵作为初始基（或者本身存在单位矩阵）的情况
%C--->目标系数向量,价值向量,按行向量输入
%A--->系数矩阵,按普通输入就行
%b--->资源向量，按列向量格式输入
%S--->返回的最优解，在没用写输入参数时会作为ans返回
%V--->返回的最优目标值
%rA--->返回运算结束时的A矩阵，供后续调用
%rb--->返回运算结束时的b向量，供后续调用
%rXB--->返回运算结束时的XB向量，供后续调用
%XB-->定位基变量
%theta--->换出准则向量
%Q--->为检验向量Q=CN-CB * N
lb = length(b); %b向量的长度
lc_list = 1:length(C);
flag = 0;
% 对于单位基矩阵不在矩阵的最后几列，可以人工交换列的次序
N = A;
CN = C;
XN = lc_list;
% 这里的操作以单位阵出现在系数矩阵右端为前提
[row_A, colomn_A] = size(A);
index_A = colomn_A - row_A;
XB = index_A+1:length(C);
XN(XB) = [];
B = A(:,XB);
N(:,XB) = [];
CB = C(XB);
CN(XB) = []; % C中剔除CB即为CN

while 1
    
    Z = CB * b;
    Q = CN - CB * N; % Q为检验向量
    index_Q_Positive = find(Q>0); % Q中可能存在多个正的元素
    
    for i = index_Q_Positive
        if max(N(:,i))<=0
            fprintf('此线性规划问题解无界！\n');
            S=[];V=[];
            flag = 1;
            break;
        end
    end
    if flag
        break; % 如果解无界的话，就跳出大循环，不用再求解了！
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%CTRL + I为批量缩进%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    X = zeros(1, length(C)); % 先把解赋成全0，然后把基向量对应的位置覆盖掉
    n = 1;
    for i = XB
        X(i) = b(n);
        n = n + 1;
    end
    
    if max(Q) < 0 %此时有唯一的最优解，其实最后显示的时候不用显示松弛变量的值，它对函数值也无影响
        S = X;
        V = Z;
        rA = A;
        rb = b;
        rXB = XB;
        disp('最优解为：'),disp(S);
        disp('最优目标值为：'),disp(V);
        break;
    end 
    
    if max(Q) == 0
        %为了后面调用过程中不显示中间过程，把输出注释掉了
        fprintf('此线性规划问题有无穷多最优解\n');
        S = X;
        V = Z;
        rA = A;
        rb = b;
        rXB = XB;
        disp('其中一个最优解为：'),disp(S);
        disp('最优目标值为：'),disp(V);
        break;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 确定换入向量（索引）
    % 要小心这里Q中可能存在多个最大值相同（还不是退化）
    % 当存在多个相同最大值时，就默认取第一个
    index_Q_Max = find(Q == max(Q));
    index_Q_Max = index_Q_Max(1);
    in_A = N(:,index_Q_Max); % in_A为换入向量
    % 确定换出向量（索引）
    theta = zeros(1,lb); %先初始化theta，防止下面有些未赋值的情况
    for i = 1:length(theta) %这里不用计算分母为0或为负数的情况
        if in_A(i) < 1e-8
            continue;
        end
        theta(i) = b(i) / in_A(i);
    end
    % 因为上面theta初始化为零向量，所以这里找的是最小的正theta
    index_theta_Min = find(theta == min(theta(theta~=0))); % 这里至少会有一项正的theta，因为没有的话已经被选中为上面的解无界的情况了
    pivot = N(index_theta_Min, index_Q_Max); % pivot即为系数矩阵中初等行变换的主元
    
    XB(index_theta_Min) = XN(index_Q_Max); % 对XB进行更新，并更新对应的系数
    CB(index_theta_Min) = CN(index_Q_Max);
    XN = lc_list;
    XN(XB) = []; % X中剔除XB即为XN
    % 下面其实也可以这样写,如果先求出XN的话
    %CN = C(XN);
    CN = C;
    CN(XB) = []; % C中剔除CB即为CN
    
    
    %%%%%%%%%接下来要做初等行变换了%%%%%%%%%
    A(index_theta_Min, :) = A(index_theta_Min, :) / pivot;
    b(index_theta_Min) = b(index_theta_Min) / pivot;
  
    for i = 1:row_A
        if i==index_theta_Min
            continue;
        end
        times = N(i,index_Q_Max);
        A(i,:) = A(i,:) - times * A(index_theta_Min, :);
        b(i) = b(i) - times * b(index_theta_Min);
    end
    
    B = A(:,XB);
    % 下面也可以这样写
    % N = A(:,XN);
    N = A;
    N(:,XB) = [];
end
end

3. twoPhaseMethod.m文件（两阶段法，只要用这个就行，会调用前两个函数）

function [ optimal_solution,optimal_value ] = twoPhaseMethod( C,A,b,num_AV )
%twoPhaseMethod 运用两阶段法，调用之前的单纯形法进行LP问题的求解
%C,A,b都是化为标准型之后的系数
%parameter：C--->目标变量系数向量,按行向量输入，其中人工变量系数为-1
%parameter：A--->系数矩阵,按普通逐行输入就行,经分析可知A几乎都是行数小于列数
%parameter：b--->资源向量，按列向量格式输入
%parameter：num_AV--->化成标准型中人工变量的个数（0<=num_AV<=row_A）
%return：optimal_solution--->最优解
%return：optimal_value--->最优值

%[row_A,colomn_A] = size(A);
length_C = length(C);
% flag = 1;
% 计算人工变量之外变量的个数
num_V = length_C - num_AV;
% if num_AV == 0
%     flag = 0;
% end

% 第一阶段：调用之前的simplexMethod.m来检验原问题是否有可行解
% 若返回的最优值为0，则进入第二阶段；否则，原问题无可行解
% 若是只通过添加松弛变量就能构成单位矩阵的，那样默认进入第二阶段（肯定有解）
% 因为其实第一阶段也只是给第二阶段找一个初始解，如果本来就有就不需要第一阶段了
C11 = zeros(1,num_V);
C12 = -1 * ones(1,num_AV);
C1 = [C11 C12];
% 因不想改之前的单纯形法的代码输入格式（要改也是可以的，就是比较大改）
% 所以这里得符合之前输入系数矩阵最后是单位矩阵的格式
% 以下就是把系数矩阵单位阵后移，并移动相应的CB和XB，这样交换一下对最优值没用影响，
% 只是最优解某些顺序换了一下（其实并不会移动人工变量）
index_A_UnitMatrix = seekUnitMatrix(A);
move_index_A = index_A_UnitMatrix(1:length(index_A_UnitMatrix)-num_AV);
length_move_index_A = length(move_index_A);
% 这里的交换不会出现把换好的又换出去的问题，下面就会出现了
for i=1:length_move_index_A
    % 交换A，C1，C
    A(:,[move_index_A(i),num_V-length_move_index_A+i]) = A(:,[num_V-length_move_index_A+i,move_index_A(i)]);
    C1([move_index_A(i),num_V-length_move_index_A+i]) = C1([num_V-length_move_index_A+i,move_index_A(i)]); %下面其实要用的是C1
    C([move_index_A(i),num_V-length_move_index_A+i]) = C([num_V-length_move_index_A+i,move_index_A(i)]); % C在这里也交换了为了保持一致
end
[rS,rV,rA,rb,rXB] = simplexMethod(C1,A,b);
for i=1:length_move_index_A
    % 交换rA
    rA(:,[move_index_A(i),num_V-length_move_index_A+i]) = rA(:,[num_V-length_move_index_A+i,move_index_A(i)]);
end
if rV ~= 0
    disp('原问题无可行解！');
    optimal_solution = [];
    optimal_value = [];
    
end

if rV == 0
    disp('原问题有可行解！');

% 第二阶段，去掉人工变量，还原目标函数系数，利用第一阶段得到初始单纯形表
AA = rA(:,1:num_V);
[row_AA,colomn_AA] = size(AA);
C = C(1:num_V);
% 此时，为了适合simplexMethod.m的输入格式，又得交换列
index_AA_UnitMatrix = seekUnitMatrix(AA);
W = 1:num_V;
W(index_AA_UnitMatrix) = [];
not_index_AA_UnitMatrix = W;
length_move_index_AA = length(index_AA_UnitMatrix);
AAA=zeros(row_AA,colomn_AA);
CCC=zeros(1,length(C));
for i=1:length_move_index_AA
    % 在凑单纯形法的输入格式
    AAA(:,colomn_AA-length_move_index_AA+i) = AA(:,index_AA_UnitMatrix(i));
    CCC(colomn_AA-length_move_index_AA+i) = C(index_AA_UnitMatrix(i));
end
for j = 1:length(not_index_AA_UnitMatrix)
    AAA(:,j) = AA(:,not_index_AA_UnitMatrix(j));
    CCC(j) = C(not_index_AA_UnitMatrix(j));
end
[finalS,finalV,~,~,~] = simplexMethod(CCC,AAA,rb);
finalSS = zeros(1,length(finalS));
% 现在得到的解还得交换顺序，好与最之前的解的顺序对应
% 交换1
for j = 1:length(not_index_AA_UnitMatrix)
    finalSS(j) = finalS(not_index_AA_UnitMatrix(j));
end
% 交换2
for i=1:length_move_index_AA
    finalSS(index_AA_UnitMatrix(i)) = finalS(colomn_AA-length_move_index_AA+i);
end
% 交换3
disp('原问题最优解为：'),disp(finalSS);
disp('原问题最优值为：'),disp(finalV);
end

end