下面的最小费用最大流算法采用的是“基于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

0