破圈法求最小支撑树matlab_【运筹学】Vogel法求运输问题初始基解的Matlab实现

function B = Vogel(A,S,D)
%实现按Vogel法生成运输问题的初始基可行解
%A为运价矩阵, S为供应向量, D为需求向量
%返回B为初始可行调运方案
M = 10000;       %足够大数
[m,n] = size(A);
B = zeros(m,n) * NaN;  %存放结果
while true
    rA = sort(A);
    cA = sort(A,2);
    %最小差值向量
    minD = [(cA(:,2) - cA(:,1))', rA(2,:) - rA(1,:)];  %前m个是最小行差, 后n个是最小列差
    [~,I] = max(minD);      %选出最大差值索引
    %选出供应位置(Ir,Ic)
    if I > m
        Ic = I - m;
        [~,Ir] = min(A(:,Ic));
    else
        Ir = I;
        [~,Ic] = min(A(Ir,:));
    end
    %进行调运
    if S(Ir) > D(Ic)             %若供应大于需求
        B(Ir,Ic) = D(Ic);
        S(Ir) = S(Ir) - D(Ic);
        D(Ic) = 0;
        A(:,Ic) = M;
    elseif S(Ir) < D(Ic)         %若供应小于需求
        B(Ir,Ic) = S(Ir);
        D(Ic) = D(Ic) - S(Ir); 
        S(Ir) = 0;
        A(Ir,:) = M;               
    else                         %若供应等于需求
        B(Ir,Ic) = S(Ir);
        S(Ir) = 0;
        D(Ic) = 0;
        A(Ir,:) = M;
        A(:,Ic) = M;
        if min(min(A)) == M
            break
        else
            ind = find(isnan(B(:,Ic))); 
            B(ind(1),Ic) = 0;           %补0  
        end
    end
end

A = [3 11 3 10; 1 9 2 8; 7 4 10 5];  %运价矩阵
S = [7 4 9];    %供应
D = [3 6 5 6];  %需求

Bv = Vogel(A,S,D)
cv = sum(sum(A .* Bv, 'omitnan'))  %总运费