对偶单纯形算法的MATLAB实现

最新推荐文章于 2023-12-16 10:29:14 发布

臻orz

最新推荐文章于 2023-12-16 10:29:14 发布

阅读量2.4k

点赞数 4

文章标签：线性代数矩阵算法动态规划

本文链接：https://blog.csdn.net/weixin_43727383/article/details/105276545

版权

一、理论

略

二、代码

function [x,y,ResultFlag]=DualSimplexAlgorithm(A,B,C,varargin)
% 2020-4-2 臻orz
%inputs:
%   A:系数矩阵 m*n
%   B:右端向量 m*1
%   C:价格系数向量 n*1
%alternative inputs:
%   target:优化目标 0 ~ min; 1 ~ max;
%   sign:约束条件符号 -1 ~ <=; 1 ~ >=;
%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',-1,@(x)validateattributes(x,...
    {'double'},{'scalar'},'BigMSimplexAlgorithm','sign'));
ip.parse(A,B,C,varargin{:});

%initialize
target = ip.Results.target;
[m,n] = size(A);
sign = repmat(ip.Results.sign,m,1);
P = [];
x = zeros(n,1);
y = 0;
ResultFlag = 0;
j = 0;

%standardization
if target
    C = -C;%目标函数的转化
end
A(B<0,:) = -A(B<0,:);
sign(B<0,:) = -sign(B<0,:);
B = abs(B);%约束条件的转化
for i = sign'
    j = j+1;
    switch i
        case -1%引入松弛变量
            a = zeros(m,1);a(j) = 1;
            A = [A a];
            C = [C;0];
        case 1%引入剩余变量
            A(j,:) = -A(j,:);
            B(j) = -B(j);
            a = zeros(m,1);a(j) = 1;
            A = [A a];
            C = [C;0];
    end
end

%找寻单位矩阵
for i = 1:m
    for j = find(A(i,:)==1)
        if sum(A(:,j)==0) == m-1
            P = [P j];
        end
    end
end

P = P(1:m);
CB = C(P);%基变量对应的价值系数
sigma = C'-CB'*inv(A(:,P))*A;
sigma(P) = 0;

while 1
    if ~sum(B<0)%有可行解
        x = zeros(size(A,2),1);
        x(P) = B;
        x = x(1:n);%舍去引入的松弛变量与剩余变量
        if target
            y = -CB'*B;
        else
            y = CB'*B;
        end
        ResultFlag = 1;
        return;
    end
    for i = find(B<0)
        if ~sum(A(i,:)<0)
            return;%无可行解
        end
    end
    pivot_x = find(B==min(B));%确定主元
    pivot_x = pivot_x(1);
    theta_index = find(A(pivot_x,:)<0);
    theta = sigma(theta_index)./A(pivot_x,theta_index);
    pivot_y = theta_index(theta==max(theta));
    pivot_y = pivot_y(1);
    P(pivot_x) = pivot_y;%更新P
    CB(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 = a
        B(i) = B(i)-A(i,pivot_y)*B(pivot_x);
        A(i,:) = A(i,:)-A(i,pivot_y)*A(pivot_x,:);    
    end
    sigma = sigma-sigma(pivot_y)*A(pivot_x,:);%更新sigma
end
end

调用方式为

 [x,y,result]=DualSimplexAlgorithm(A,B,C,'target',0,'sign',-1);

其中A，B，C分别为系数矩阵，右端向量，价格系数向量，注意B，C为列向量即可，后面的两个为可选输入
‘target’：为优化的目标为求最小值还是最大值，分别对应0和1，可以不设置，默认为0，即求最小值
‘sign’：为约束条件的符号，小于、大于分别对应-1,1，是一常值，可以不设置，默认为-1，即全为小于