[最优化方法/运筹学] 大M法解决线性规划代码 单纯形法

问题：

代码：

M=1e6;
f=[-3,2,1,0,0,M,M]; %加入松弛变量和人工变量 
A=[  1,-2,1,1,0,0,0;
    -4,1,2,0,-1,1,0;
    -2,0,1,0,0,0,1]; %x1-2x2+x3<=11
B=[11; 3; 1];           %4x1-x2-2x3<=-3
% Aeq=[-2,0,1];        %-2x1+x3=1
% Beq=1;

% 调用大M方法求解
[x_opt, f_opt] = big_m_method(f, A,B);

% 输出结果
disp('最优解:');
disp(x_opt);
disp('最优目标函数值:');
disp(f_opt);

A f B是根据问题写的 松弛变量人工变量是手动加的

函数：

function [x_opt, f_opt] = big_m_method(c, A, b)
    % 输入参数:
    % c: 目标函数系数向量
    % A: 约束条件系数矩阵
    % b: 约束条件右侧常数向量
    % 输出参数:
    % x_opt: 最优解向量
    % f_opt: 最优目标函数值
    cnt=1;
    Xb=get_Xb(A); %这个有问题 为什么第一次得到的都是零呢 
    while cnt<5
        % disp(A);
        lamda =get_lamda(c,A,Xb); %得到检验数 
        % disp('c');
        % disp(c);
        % disp('lamda');
        % disp(lamda);
        % disp('Xb');
        % disp(Xb);
        cnt=cnt+1;
        if all(lamda>=0) %如果检验数都大于等于0 
            break
        end
        %循环结束之后 就要用lamda里面找到最小的作为入基 
        [~,inx_index]=min(lamda); %得到入基
        %得到入基之后要算出sita 找到最小 得到出基
        sita=get_sita(A,b,inx_index);
        % disp('sita');
        % disp(sita);
        [~,outx_index]=min(sita); %得到出基 出基变量 Xb(outx_index) 使用就是用outx
        %得到入基和出基之后 通过初等行变换把入基的列向量变为只有0和1且和为1的那种形式 A和b都要变
        [A,b]=change_matrix(A,b,inx_index,outx_index);
        %在原来的出基的位置直接变成入基
        Xb=change_Xb(Xb,outx_index,inx_index);
        %进入循环 重新计算lamda
    end
    %循环结束之后 根据xb的位置和b得到最后的结果
    x_opt=zeros(size(A,2),1);
    for i=1:size(Xb,1)
        x_opt(Xb(i))=b(i);
    end
    f_opt=c*x_opt;
end

function [A,b]=change_matrix(A,b,inx,outx)
    num=A(outx,inx);
    A(outx,:)=A(outx,:)/num;b(outx)=b(outx)/num; %这个位置变成1
    for i=1:size(A,1)
        if i==outx
            continue;
        else
            q=A(i,inx)/A(outx,inx);
            A(i,:)=A(i,:)-q*A(outx,:);
            b(i,:)=b(i,:)-q*b(outx,:); %初等行变换 其他行的这个位置变成0
        end
    end
end

function lamda=get_lamda(c,A,Xb)
    lamda=c; %检验数
    for i=1:size(A,1)
        % disp(lamda);
        % disp(A(i,:));
        % disp(c(Xb(i)));
        lamda=lamda-A(i,:)*c(Xb(i));
    end
end

function Xb=get_Xb(A)
    Xb=zeros(3,1);
    index=1;
    for v=1:size(A,2)
        isvalid=judge(A,v);
        if isvalid
            Xb(index)=v;
            disp(Xb(index));
            index=index+1;
        end
    end
end

function Xb=change_Xb(Xb,outx,inx)
    Xb(outx)=inx;
end

function sita=get_sita(A,b,min_index)
    %得到入基之后要算出sita 找到最小 得到出基
    p=A(:,min_index);
    % disp('p');
    % disp(p);
    cnt=size(p,1);
    sita=zeros(cnt,1);
    for i=1:size(b,1)
        if p(i)>0
            sita(i)=b(i)/p(i);
        else 
            sita(i)=1e6; %还是需要保持形状一致 
        end
    end %得到sita
end
    
function isValid=judge(A,v)
    validElements = all(ismember(A(:,v), [0, 1]));
    % 检查元素和是否等于 1
    sumElements = sum(A(:,v)) == 1;
    % 判断是否同时满足两个条件
    isValid = validElements && sumElements;
end

结果：

化为标准型目标函数取了负 最后结果没问题

很多disp做debug用 但是最后怎么改对的我也不知道