下面的最小費用最大流算法采用的是“基于Floyd最短路算法的Ford和Fulkerson疊加算法”，其基本思路為：把各條弧上機關流量的費用看成某種長度，用Floyd求最短路的方法确定一條自V1至Vn的最短路；再将這條最短路作為可擴充路，用求解最大流問題的方法将其上的流量增至最大可能值；而這條最短路上的流量增加後，其上各條弧的機關流量的費用要重新确定，如此多次疊代，最終得到最小費用最大流。本源碼由GreenSim團隊原創，轉載請注明
function [f,MinCost,MaxFlow]=MinimumCostFlow(a,c,V,s,t)
%% MinimumCostFlow.m
% 最小費用最大流算法通用matlab函數
%% 基于Floyd最短路算法的Ford和Fulkerson疊加算法
% GreenSim團隊原創作品，轉載請注明
%% 輸入參數清單
% a 機關流量的費用矩陣
% c 鍊路容量矩陣
% V 最大流的預設值，可為無窮大
% s 源節點
% t 目的節點
%% 輸出參數清單
% f 鍊路流量矩陣
% MinCost 最小費用
% MaxFlow 最大流量
%% 第一步：初始化
N=size(a,1);%節點數目
f=zeros(N,N);%流量矩陣，初始時為零流
MaxFlow=sum(f(s,:));%最大流量，初始時也為零
flag=zeros(N,N);%真實的前向邊應該被記住
for i=1:N
for j=1:N
if i~=j&&c(i,j)~=0
flag(i,j)=1;%前向邊标記
flag(j,i)=-1;%反向邊标記
end
if a(i,j)==inf
a(i,j)=BV;
w(i,j)=BV;%為提高程式的穩健性，以一個有限大數取代無窮大
end
end
end
if L(end)
RE=1;%如果路徑長度小于大數，說明路徑存在
else
RE=0;
end
%% 第二步：疊代過程
while RE==1&&MaxFlow<=V%停止條件為達到最大流的預設值或者沒有從s到t的最短路
%以下為更新網絡結構
MinCost1=sum(sum(f.*a));
MaxFlow1=sum(f(s,:));
f1=f;
TS=length(R)-1;%路徑經過的跳數
LY=zeros(1,TS);%流量裕度
for i=1:TS
LY(i)=c(R(i),R(i+1));
end
maxLY=min(LY);%流量裕度的最小值，也即最大能夠增加的流量
for i=1:TS
u=R(i);
v=R(i+1);
if flag(u,v)==1&&maxLY
f(u,v)=f(u,v)+maxLY;%記錄流量值
w(u,v)=a(u,v);%更新權重值
c(v,u)=c(v,u)+maxLY;%反向鍊路的流量裕度更新
elseif flag(u,v)==1&&maxLY==c(u,v)%當這條邊為前向邊且是飽和邊時
w(u,v)=BV;%更新權重值
c(u,v)=c(u,v)-maxLY;%更新流量裕度值
w(v,u)=-a(u,v);%反向鍊路權重更新
elseif flag(u,v)==-1&&maxLY
w(v,u)=a(v,u);
c(v,u)=c(v,u)+maxLY;
w(u,v)=-a(v,u);
elseif flag(u,v)==-1&&maxLY==c(u,v)%當這條邊為反向邊且是飽和邊時
w(v,u)=a(v,u);
c(u,v)=c(u,v)-maxLY;
w(u,v)=BV;
else
end
end
MaxFlow2=sum(f(s,:));
MinCost2=sum(sum(f.*a));
if MaxFlow2<=V
MaxFlow=MaxFlow2;
MinCost=MinCost2;
[L,R]=FLOYD(w,s,t);
else
f=f1+prop*(f-f1);
MaxFlow=V;
MinCost=MinCost1+prop*(MinCost2-MinCost1);
return
end
if L(end)
RE=1;%如果路徑長度小于大數，說明路徑存在
else
RE=0;
end
end
function [L,R]=FLOYD(w,s,t)
n=size(w,1);
D=w;
path=zeros(n,n);
%以下是标準floyd算法
for i=1:n
for j=1:n
if D(i,j)~=inf
path(i,j)=j;
end
end
end
for k=1:n
for i=1:n
for j=1:n
if D(i,k)+D(k,j)
D(i,j)=D(i,k)+D(k,j);
path(i,j)=path(i,k);
end
end
end
end
L=zeros(0,0);
R=s;
while 1
if s==t
L=fliplr(L);
L=[0,L];
return
end
L=[L,D(s,t)];
R=[R,path(s,t)];
s=path(s,t);
end
