最优化-线性规划-单纯形法2

单纯形法转化成标准式后，便可解决通过切换基可行解即凸集顶点进行优化，依据当前判断状态来更迭计算，此处直接给出两阶段法的单纯形法；

clc;clear;close
%% 问题描述-线性规划
%目标函数的系数-目标是最大还是最小-约束条件<=\<不同需要区分-决策变量的取值范围
% 适用面太窄
% 线性规划解决的闭合集合，非闭合则确定极限
%%  输入数据
amb=1;%1 :目标为Max；   -1 :目标函数为Min
cz=rands(1,50).*100;%目标函数的决策变量系数
ca=rand(20,50).*20;%约束条件决策变量的系数
cb=rand(1,20).*100;%约束条件所对应的右项值
lrel=round(rands(1,20));%约束条件左项式与右项式的关系
% <=小于等于 -1 ;== 等于 0; >= 大于等于1
ul=[[-inf,-inf,round(rands(1,48).*100)];[20,inf,120,ones(1,47).*inf]];%决策变量的最小值\最大值
save 'cz.mat' cz;save 'ca.mat' ca;save 'cb.mat' cb;save 'lrel.mat' lrel;
save 'ul.mat' ul;
% load('cz.mat');load('ca.mat');load('cb.mat');load('lrel.mat');load('ul.mat');
%% 以上都需要缺一不可
tic
%STEP 1 --转化成标准形式 
[zlast,xlast,newcz,newca,newcb,addn,deln]=ori2sta(amb,cz,ca,cb,lrel,ul);
%STEP 2--大M法转化成新的标准形式-形成初始解
[mcz,mca,mcb,X0,manualn,xp]=bigm(newcz,newca,newcb,addn);
%STEP 3--大M求解验证是否有可行解
% 1 检验解的情况 是否可能有更好的解 2 如果有转移基可行解，转移方式，
sum(newcb);
[f,texta,mca,mcb,X0,x1,xp,lastf]=simpl(mcz,mca,mcb,X0,xp);
%STEP 4-- 转变成需要求的方程
judge=0;
if lastf(end)==0 %如果存在最佳值
    % 首先检验 根是否还有大M法中的人工变量
    limb=numel(mcb);%当前的约束条件的变化
    n=0;
    dela=[];
    for i=1:limb
        if sum(manualn(:,2)'==xp(i,2))>0 %如果是之前的
            n=n+1;
            dela(n)=i;
        end
    end
    if n>0 %存在人工变量，但此约束无用
        mca(dela,:)=[];
        mcb(dela)=[];
        xp(dela,:)=[];
    end
    mca(:,manualn(:,2)')=[];
    dema=numel(newcz);%当前的决策变量个数
    X0=zeros(1,dema);
    X0(xp(:,2)')=1;
    limb=numel(mcb);%当前的约束条件的变化
    xp(:,1)=(1:limb)';
    [f,texta,mca,mcb,X0,x1,xp,lastf]=simpl(newcz,mca,mcb,X0,xp);
else
    judge=1;
    '无解'
end
% STEP 5 --转化成需要的值
if judge==0 % 有解
    limb=size(xlast,2);
    lastx=zeros(1,limb);
    for i=1:limb
        buf=xlast{i};
        for j=1:size(buf,2)
            lastx(i)=lastx(i)+x1(buf(1,j))*buf(2,j)+buf(3,j);
        end
    end
    lastz=(lastf(end)+zlast(2))*zlast(1);
end

toc

function [zlast,xlast,cz,ca,cb,addn,deln]=ori2sta(amb,cz,ca,cb,lrel,ul)
% author:major-岩寐
% decisionnum 决策变量数量 limitnum约束条件限制个数
% coef_A 约束条件系数 coef_B右项式 coef_C决策变量系数 
% limitrelation 约束条件左项式和右项式的关系 -1：=《 ， 0 ：= = ， 1：= 》
% dicisionrelation 约束条件的取值-1=《0 ， 0 ：无界 ， 1=
% 》0（若决策变量大于或小于某值，根据取值范围来判断，譬如，x1>=-1,则将该约束写入约束条件，x1取值为无界
%默认 使用两阶段求解，来产生初始解，大M法不可靠在计算机里
zlast=[1,0];
dema=numel(cz);%当前的决策变量个数
xlast=cell(1,dema);
for i=1:dema
    xlast{i}=[i;1;0];
end %转化后与之前的关系明确
%% 1 转化决策变量--
dema=numel(cz);%当前的决策变量个数
limb=numel(cb);%当前的约束条件的变化
n=0;%新的决策变量数量
m=0;%约束条件的变化
for i=1:dema %依据取值范围进行求解
    if ul(1,i)>-inf % 如果最小值大于无限最小 x0>=ul(1,i),x1+ul(1,i)>=ul(1,i)
        % 关系变化
        xlast{i}=[i;1;ul(1,i)];
        zlast=[1,zlast(2)+ul(1,i)*cz(i)];%目标函数
        % 更新其他 cz ca cb  lrel
        cb=cb-ul(1,i).*ca(:,i)';
        if ~(ul(2,i)<inf)%表明取值范围合理
        else% 增加一个约束
            m=m+1;%增加一个约束
            buf=zeros(1,dema+n);%决策变量
            buf(i)=1;
            ca(limb+m,:)=buf;
            cb(limb+m)=ul(2,i)-ul(1,i);
            lrel(limb+m)=-1;%小于等于
        end
    else %如果可以取值无穷小
        if ~(ul(2,i)<inf)%表明取值为实数 则需要添加新变量 x0=x1-x2>-inf<inf
            n=n+1;%添加新变量
            % 关系变化
            xlast{i}=[i,dema+n;1,-1;0,0];%替换关系
            % 更新其他 cz ca cb 
            cz(dema+n)=-cz(i);%决策系数
            ca(:,dema+n)=-1.*ca(:,i);
        else %最大值有限制，num>x0>-inf,num>(-x1+num)>-inf,0<x1-num<inf,
            % 关系变化
            xlast{i}=[i;-1;ul(2,i)];
            zlast=[1,zlast(2)+ul(2,i)*cz(i)];%目标函数
            % 更新其他 cz ca cb 
            cz(i)=cz(i)*-1;
            cb=cb-ul(2,i).*ca(:,i)';
            ca(:,i)=ca(:,i).*-1;
        end
    end
end
%% 2 转化目标函数 min
if amb~=-1 %如果不是min
    % 关系变化
    zlast=-1.*zlast;%目标函数
    % 更新其他 cz ca cb 
    cz=cz.*-1;
end
%% 3 转化右项式 以及 左右项式关系 先变动右项式再变动项式之间关系，不然松弛变量的系数为-1
n=0;%新添加的+项变量 
addn=[];
m=0;%新添加的-项变量
deln=[];
dema=numel(cz);%当前的决策变量个数
limb=numel(cb);%当前的约束条件的变化
for i=1:limb
    if cb(i)<0% 如果右项式小于0 小于 后面的
        % 更新其他 ca cb lrel
        ca(i,:)=ca(i,:).*-1;cb(i)=cb(i).*-1;lrel(i)=lrel(i)*-1;
    end
    if lrel(i)==-1
        n=n+1;
        % 更新其他 cz ca cb lrel
        cz(dema+m+n)=0;
        ca(:,dema+m+n)=zeros(limb,1);
        ca(i,dema+m+n)=1;
        addn(n,:)=[i,dema+m+n];
    elseif lrel(i)==1
        m=m+1;
        % 更新其他 cz ca cb lrel
        cz(dema+m+n)=0;
        ca(:,dema+m+n)=zeros(limb,1);
        ca(i,dema+m+n)=-1;
        deln(m,:)=[i,dema+m+n];
    end
end
% relation.store=store;relation.stylend=stylend;relation.intialcase=intialcase;
% standard.decisionnum=n;standard.limitnum=limitnum;standard.coef_A=newcoef_A;standard.coef_B=coef_B;standard.coef_C=newcoef_C;
end

function [f,texta,mca,mcb,X0,x1,xp,lastf]=simpl(mcz,mca,mcb,X0,xp)
f=0;
n=0;
lastf=[];
while f==0
    [f,texta,mca,mcb,X0,x1,xp,zhi]=checkx(mcz,mca,mcb,X0,xp);
    n=n+1;
    lastf(n)=zhi;
end
end

function [f,texta,mca,mcb,X0,x1,xp,zhi]=checkx(mcz,mca,mcb,X0,xp)
%% 检验值的原理-弱确定了基可行解，则将基可行解带入到目标函数中来判断
%% 所以验证的前提在于约束矩阵已经形成了基可行解的事态
% 计算检验值--存在退化或者说循环
f=0;
texta='';
limb=numel(mcb);%当前的约束条件的变化
testn=mcz-mcz(xp(:,2)')*mca;
nok=find(X0==0);%非基项
[num,~]=min(testn(nok));%最小值
if num<0 %如果存在小于0 则可能还有更小的值
    % 判断是否是无界解
    if min(sum(mca(:,nok(testn(nok)<0))>0,1))==0 %如果存在增加该值大于0，但约束系数均为0
        f=1;%
        texta='无界';
    else
        [ink,our,ouk,num1]=bio(nok,testn,mca,mcb,xp);%进出明确
        %% 更新
        buf=1:limb;buf(buf==our)=[];
        mcb(our)=mcb(our)/mca(our,ink);
        mcb(buf)=mcb(buf)- mcb(our).*mca(buf,ink)';
        mca(our,:)=mca(our,:)./mca(our,ink);
        mca(buf,:)=mca(buf,:)-mca(buf,ink)*mca(our,:);
        xp(our,2)=ink;
        X0(ouk)=0;X0(ink)=1;
    end
else
    if num==0
        f=2;
        texta='无穷最优解';
    else
        f=3;
        texta='唯一最优解';
    end
end
x1=X0;
x1(xp(:,2)')=mcb;
zhi=sum(mcz.*x1);
end
function [newcz,newca,newcb,X0,manualn,xp]=bigm(cz,ca,cb,addn)
%% 大M法 首先添加人工变量 形成真正的标准型
dema=numel(cz);%当前的决策变量个数
limb=numel(cb);%当前的约束条件的变化
maln=limb-size(addn,1);%需要添加的人工数
manualn=zeros(maln,2);
newcz=zeros(1,dema+maln);
X0=zeros(1,dema+maln);
newca=ca;newcb=cb;
n=0;%添加的人工变量
for i=1:limb
    if sum(addn(:,1)==i)==0
        n=n+1;
        newcz(dema+n)=1;
        newca(:,dema+n)=zeros(limb,1);
        newca(i,dema+n)=1;
        manualn(n,:)=[i,dema+n];
    end
end
X0(addn(:,2)')=1;X0(manualn(:,2)')=1;
xp=[addn;manualn];
[~,buf]=sort(xp(:,1)');
xp=xp(buf,:);
end
 function [ink,our,ouk,num]=bio(nok,testn,mca,mcb,xp)
%选择入基 出基
%% 最大入基
[num1,index]=min(testn(nok));
ink=nok(index);
%% 确定出基
buf=find(mca(:,ink)>0);
buf1=mcb(buf)./mca(buf,ink)';
[num,index1]=min(buf1);
our=buf(index1);
ouk=xp(our,2);
% num1*num    
end