【程序】求网络的最小费用最大流

宝宝 · 发表于 2004-6-10 08:36:14

<

>求网络的最小费用最大流，弧旁的数字是容量(运费)。
<IMG src="http://vip.6to23.com/yunyan8/shu/images/1.jpg">
一．Ford和Fulkerson迭加算法.
基本思路：把各条弧上单位流量的费用看成某种长度,用求解最短路问题的方法确定一条自V1至Vn的最短路;在将这条最短路作为可扩充路,用求解最大流问题的方法将其上的流量增至最大可能值;而这条最短路上的流量增加后,其上各条弧的单位流量的费用要重新确定,如此多次迭代,最终得到最小费用最大流.
迭加算法：
1) 给定目标流量F或∞,给定最小费用的初始可行流{fij}=0
2) 若V(f)=F,停止,f为最小费用流;否则转(3).
3) 构造{fij} 相应的新的费用有向图W(fij),在W(fij)寻找Vs到Vt的最小费用有向路P(最短路),沿P增加流f的流量直到F,转(2);若不存在从Vs到Vt的最小费用的有向路P,停止.f就是最小费用最大流.
具体解题步骤：
设图中双线所示路径为最短路径，费用有向图为W（fij）．

在图(a)中给出零流 f,在图(b)中找到最小费用有向路,修改图(a)中的可行流,δ=min{4,3,5}=3,得图(c),再做出(c)的调整容量图,再构造相应的新的最小费用有向路得图(d), 修改图(c)中的可行流, δ=min{1,1,2,2}=1,得图(e),以此类推,一直得到图(h),在图(h)中以无最小费用有向路,停止,经计算:
图(h)中最小费用=1*4+3*3+2*4+4*1+1*1+4*2+1*1+3*1=38
图(g)中最大流=5 
<

>C语言程序如下(运行通过):
/*Ford和Fulkerson迭加算法*/
#include"stdio.h"
void main()
{int a,b,i,j,k,p,n,B[10][10],C[10][10],D[10][10],P[10][10][10];
printf("please input n:\n");
scanf("%d",&n);
printf("please input C[%d][%d],B[%d][%d]:\n",n,n,n,n);
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
scanf("%7d,%7d",&C[j],&B[j]); //输入各点(i,j)之间的容量C[j]和运费B[j]
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
{D[j]=B[j];
for(k=0;k<=n;k++)
P[j][k]=0;
if(D[j]<100) //注：100表示Vi到Vj之间无可行路
{P[j]=1[j][j]=1;}
}
for(k=0;k<=n;k++)
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
if(D[k]+D[k][j]<D[j])
{D[j]=D[k]+D[k][j];
for(p=0;p<=n;p++)
P[j][p]=P[k][p]||P[k][j][p];
} //调用FLOYD算法
printf("D[%d][%d]:\n",n,n);
for(i=0;i<=n;i++)
{for(j=0;j<=n;j++)
printf("%7d",D[j]);
printf("\n");
} //由表D中的数值确定V0到V5的最短路
a=C[1][3];b=B[1][3]*D[0][5];
printf("a=%d,b=%d\n",a,b);
B[1][3]=100; //将容量已满的路改为不可行路
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
{D[j]=B[j];
for(k=0;k<=n;k++)
P[j][k]=0;
if(D[j]<100)
{P[j]=1;P[j][j]=1;}
}
for(k=0;k<=n;k++)
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
if(D[k]+D[k][j]<D[j])
{D[j]=D[k]+D[k][j];
for(p=0;p<=n;p++)
P[j][p]=P[k][p]||P[k][j][p];
} //调用FLOYD算法
printf("D[%d][%d]:\n",n,n);
for(i=0;i<=n;i++)
{for(j=0;j<=n;j++)
printf("%7d",D[j]);
printf("\n");
} //由表D中的数值确定V0到V5的新最短路
a=a+C[1][2];b=b+D[0][5];
printf("a=%d,b=%d\n",a,b);
B[0][1]=100;B[1][2]=100;B[4][3]=100; //将容量已满的路改为不可行路
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
{D[j]=B[j];
for(k=0;k<=n;k++)
P[j][k]=0;
if(D[j]<100)
{P[j]=1;P[j][j]=1;}
}
for(k=0;k<=n;k++)
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
if(D[k]+D[k][j]<D[j])
{D[j]=D[k]+D[k][j];
for(p=0;p<=n;p++)
P[j][p]=P[k][p]||P[k][j][p];
} //调用FLOYD算法
printf("D[%d][%d]:\n",n,n);
for(i=0;i<=n;i++)
{for(j=0;j<=n;j++)
printf("%7d",D[j]);
printf("\n");
} //由表D中的数值确定V0到V5的新最短路
a=a+C[2][4]-C[4][3];b=b+D[0][5];
printf("a=%d,b=%d\n",a,b);
B[2][4]=100; //将容量已满的路改为不可行路
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
{D[j]=B[j];
for(k=0;k<=n;k++)
P[j][k]=0;
if(D[j]<100)
{P[j]=1;P[j][j]=1;}
}
for(k=0;k<=n;k++)
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
if(D[k]+D[k][j]<D[j])
{D[j]=D[k]+D[k][j];
for(p=0;p<=n;p++)
P[j][p]=P[k][p]||P[k][j][p];
} //调用FLOYD算法
printf("D[%d][%d]:\n",n,n);
for(i=0;i<=n;i++)
{for(j=0;j<=n;j++)
printf("%7d",D[j]);
printf("\n");
} //检验有无V0到V5的新最短路
}

运行结果:
please input n:
5
please input C[5][5],B[5][5]:
0,0 4,1 5,3 0,100 0,100 0,100
0,100 0,0 1,1 3,3 0,100 0,100
0,100 0,100 0,0 0,100 2,4 0,100
0,100 0,100 0,100 0,0 0,100 5,2
0,100 0,100 0,100 1,1 0,0 2,4
0,100 0,100 0,100 0,100 0,100 0,0
D[5][5]:
0 1 2 4 6 6
100 0 1 3 5 5
100 100 0 5 4 7
100 100 100 0 100 2
100 100 100 1 0 3
100 100 100 100 100 0
a=3,b=18
D[5][5]:
0 1 2 7 6 9
100 0 1 6 5 8
100 100 0 5 4 7
100 100 100 0 100 2
100 100 100 1 0 3
100 100 100 100 100 0
a=4,b=27
D[5][5]:
0 100 3 100 7 11
100 0 100 100 100 100
100 100 0 100 4 8
100 100 100 0 100 2
100 100 100 100 0 4
100 100 100 100 100 0
a=5,b=38
D[5][5]:
0 100 3 100 100 100
100 0 100 100 100 100
100 100 0 100 100 100
100 100 100 0 100 2
100 100 100 100 0 4
100 100 100 100 100 0
//注：100表示Vi到Vj之间无可行路
　
二．圈算法：
1) 利用Ford和Fulkson标号算法找出流量为F(<=最大流)的流f.
2) 构造f对应的调整容量的流网络N'(f).
3) 搜索N'(f)中的负费用有向图C(Floyd算法),若没有则停止,否则转(4).
4) 在C上找出最大的循环流,并加到N上去,同时修改N'(F)中C的容量,转(3).
具体解题步骤:

利用Ford和Fulkson标号算法,得网络的最大流F=5,见图(a),由图(a)构造相应的调整容量的流网络(b),图(b)中不存在负费用有向图,故停止．经计算:
图(b)中最小费用= 4*1+3*1+1*1+3*3+4*2+1*1+4*1+2*4=38
图(a)中最大流为F=5
　
C语言程序如下(运行通过):
/*圈算法*/
#include"stdio.h"
int min(int x,int y)
{if(x<y) return(x);
else return(y);
}
void main()
{int a=0,b=0,i,j,k,p,n,t,A[10][10],P[10][10],B[10][10],C[10][10],D[10][10];
printf("please input n:\n");
scanf("%d",&n);
printf("please input C[%d][%d],B[%d][%d]:\n",n,n,n,n);
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
scanf("%7d,%7d",&C[j],&B[j]); //输入各点(i,j)之间的容量C[j]和运费B[j]
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
{A[j]=C[j];D[j]=0;P[j]=B[j];}
for(i=n;i>1;i--)
for(j=i-1;j>0;j--)
for(k=j-1;k>=0;k--)
if(A[j]!=0&&A[k][j]!=0)
D[k][j]=min(A[j],A[k][j]);
printf("D[%d][%d]:\n",n,n);
for(i=0;i<=n;i++)
{for(j=0;j<=n;j++)
printf("%7d",D[j]);
printf("\n");
}
for(i=0;i<n-1;i++)
for(j=i+1;j<n;j++)
for(k=j+1;k<=n;k++)
if(D[j]!=0&&D[j][k]!=0)
if(D[j]+D[j][k]==C[j])
{A[j]=0;A[j][k]=0;A[j-i][k-j]=0;
P[j]=100;P[j][k]=100;P[j-i][k-j]=0;
a=a+C[j];
b=b+B[j]*C[j]+B[j][k]*C[j][k]+B[j-i][k-j]*C[j-i][k-j];
for(p=k+1;p<=n;p++)
if(C[j]<C[k][p])
{b=b+B[k][p]*C[j];
A[k][p]=C[k][p]-C[j];
}
for(p=k-j+1;p<=n;p++)
if(C[j-i][k-j]<C[k-j][p])
{b=b+B[k-j][p]*C[j-i][k-j]+B[4][3]*C[4][3];
A[k-j][p]=C[k-j][p]-C[j-i][k-j];
}
A[4][3]=0;
}
printf("a=%d,b=%d\n",a,b);
for(i=0;i<=n;i++)
for(j=0;j<=n;j++)
D[j]=0;
for(i=n;i>1;i--)
for(j=i-1;j>0;j--)
for(k=j-1;k>=0;k--)
if(A[j]!=0&&A[k][j]!=0)
D[k][j]=min(A[j],A[k][j]);
printf("D[%d][%d]:\n",n,n);
for(i=0;i<=n;i++)
{for(j=0;j<=n;j++)
printf("%7d",D[j]);
printf("\n");
}
for(i=0;i<n-1;i++)
for(j=i+1;j<n;j++)
for(k=j+1;k<=n;k++)
if(D[j]!=0&&D[j][k]!=0)
if(D[j]==D[j][k])
for(p=k+1;p<=n;p++)
{t=min(min(A[j],A[j][k]),min(A[j][k],A[k][p]));
a=a+t;
b=b+t*(B[j]+B[j][k]+B[k][p]);
}
printf("a=%d,b=%d\n",a,b);
}

运行结果:
please input n:
5
please input C[5][5],B[5][5]:
0,0 4,1 5,3 0,100 0,100 0,100
0,100 0,0 1,1 3,3 0,100 0,100
0,100 0,100 0,0 0,100 2,4 0,100
0,100 0,100 0,100 0,0 0,100 5,2
0,100 0,100 0,100 1,1 0,0 2,4
0,100 0,100 0,100 0,100 0,100 0,0
D[5][5]:
0 1 2 0 0 0
0 0 1 3 0 0
0 0 0 0 2 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
a=4,b=27
D[5][5]:
0 0 1 0 0 0
0 0 0 0 0 0
0 0 0 0 1 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
a=5,b=38
//注：100表示Vi到Vj之间无可行路

knight · 发表于 2004-6-10 21:01:35

<

>//我做的最小费最短距离用最大流，是在原来的最小费用最大流基础上我做的改进算法<

>/*
 Name: NetWorkFlow.h
 Copyright: Knight_Studio
 Author: Ran NI
 Date: 22-05-04 09:40
 Description: The Declaration of the NewWorkFlow Class
*/
#include <vector>
#ifndef _NETWORKFLOW_H_
#define _NETWORKFLOW_H_<

>using namespace std;class NetWorkFlow
{
 public:
 NetWorkFlow::NetWorkFlow(vector<vector<double> > L,vector<vector<int> > BB,vector<vector<double> > cc,vector<vector<double> > ff,vector<double> aa,vector<double> bb):U(L),B(BB),F(ff),a(aa),b(bb),C(cc),x(aa.size())
 {
 if(Weight.size()!=0)
 {
 Weight.clear();
 }
 if(Connect.size()!=0)
 {
 Connect.clear();
 }
 if(Bound.size()!=0)
 {
 Bound.clear();
 }
 if(Length.size()!=0)
 {
 Length.clear();
 }
 for(int i=0;i<(int)(a.size()+b.size()+2);++i)
 {
 vector<double> zz(a.size()+b.size()+2,0);
 vector<int> c(a.size()+b.size()+2,0);
 Connect.push_back(c);
 Weight.push_back(zz);
 Length.push_back(zz);
 Bound.push_back(zz);
 }
 for(int i=0;i<(int)B.size();++i)
 {
 for(int j=0;j<(int)B.size();++j)
 {
 if(B[j]!=0)
 {
 Weight[i+1][j+a.size()+1]=C[j];
 Length[i+1][j+a.size()+1]=F[j];
 Bound[i+1][j+a.size()+1]=U[j];
 Connect[i+1][j+a.size()+1]=1;
 }
 }
 }
 for(int i=0;i<(int)a.size()+1;++i)
 {
 Connect[0]=1;
 if(i>0)
 {
 Bound[0]=a[i-1];
 }
 }
 for(int i=(int)Connect.size()-1;i>(int)Connect.size()-b.size()-2;--i)
 {
 Connect[Connect.size()-1]=1;
 if(i<(int)Connect.size()-1)
 {
 Bound[(int)Connect.size()-1]=b[(int)b.size()-((int)Connect.size()-i-1)];
 }
 }
 }
 vector<int> Find();
 void Augment(vector<int>);
 double LowestCost();
 vector<vector<double> > GetOptimalSolution()
 {
 return Solution;
 }

 protected:
 vector<vector<int> > B;
 vector<vector<double> > U;
 vector<vector<double> > C;
 vector<vector<double> > F;
 vector<vector<double> > Weight;
 vector<vector<int> > Connect;
 vector<vector<double> > Bound;
 vector<vector<double> > Length;
 vector<double> a;
 vector<double> b;
 vector<vector<double> > x;
 vector<vector<double> > Solution;
};
#endif

knight · 发表于 2004-6-10 21:02:35

<

>/*
 Name: NetWorkFlow.cpp
 Copyright: Knight_Studio
 Author: Ran NI
 Date: 22-05-04 09:40
 Description: The Definition of the NewWorkFlow Class
*/
#include <vector>
#include <algorithm>
#include <functional>
#include <iostream>
#include <math.h>
#include "NetWorkFlow.h"<

>using namespace std;<

>bool Judge(vector<double> a,vector<double> b)
{
 bool result=true;
 for(int i=0;i<(int)a.size();++i)
 {
 if(a!=b)
 {
 result=false;
 }
 }
 return result;
}vector<int> NetWorkFlow::Find()
{
 vector<vector<int> > path(Connect.size());
 path[0].push_back(0);
 double inf=0;
 inf=1/inf;
 vector<int> contain;
 contain.push_back(0);
 vector<double> dis(path.size(),inf);
 vector<double> len(Connect.size(),inf);
 len[0]=0;
 dis[0]=0;
 vector<double> dis1(path.size(),inf);
 dis1[0]=0;
 do
 {
 dis=dis1;
 for(int j=1;j<(int)Connect.size();++j)
 {
 int nn=0;
 int it=0;
 double ls=inf;
 double ll=inf;
 for(int k=0;k<(int)contain.size();++k)
 {
 int s=contain[k];
 if(Connect[j]!=0)
 {
 if(dis+Weight[j]<=ll)
 {
 if(len+Length[j]<ls)
 {
 ll=dis+Weight[j];
 nn=s;
 ls=len+Length[j];
 ++it;
 }
 else if(dis+Weight[j]<ll)
 {
 ll=dis+Weight[j];
 nn=s;
 ls=len+Length[j];
 ++it;
 }
 }
 }
 }
 if(it!=0)
 {
 dis1[j]=ll;
 len[j]=ls;
 path[j]=path[nn];
 path[j].push_back(j);
 if(find(contain.begin(),contain.end(),j)==contain.end())
 {
 contain.push_back(j);
 }
 }
 }
 }while(!Judge(dis,dis1));
 return path[path.size()-1];
}void NetWorkFlow::Augment(vector<int> path)
{
 double delta=Bound[path[0]][path[1]]-x[path[0]][path[1]];
 for(int i=1;i<(int)path.size()-1;++i)
 {
 if(path<path[i+1])
 {
 double xx=Bound[path][path[i+1]]-x[path][path[i+1]];
 delta=min(xx,delta);
 }
 else
 {
 double xx=x[path[i+1]][path];
 delta=min(xx,delta);
 }
 }
 for(int i=0;i<(int)path.size()-1;++i)
 {
 if(path<path[i+1])
 {
 x[path][path[i+1]]+=delta;
 if(Bound[path][path[i+1]]-x[path][path[i+1]]==0)
 {
 Connect[path][path[i+1]]=0;
 Weight[path[i+1]][path]=(-Weight[path][path[i+1]]);
 Length[path[i+1]][path]=(-Length[path][path[i+1]]);
 }else if(x[path][path[i+1]]>0)
 {
 Connect[path][path[i+1]]=1;
 Connect[path[i+1]][path]=1;
 Weight[path[i+1]][path]=(-Weight[path][path[i+1]]);
 Length[path[i+1]][path]=(-Length[path][path[i+1]]);
 }
 }
 else
 {
 x[path[i+1]][path]-=delta;
 if(Bound[path[i+1]][path]-x[path[i+1]][path]==0)
 {
 Connect[path[i+1]][path]=0;
 Weight[path][path[i+1]]=(-Weight[path[i+1]][path]);
 Length[path][path[i+1]]=(-Length[path[i+1]][path]);
 }else if(x[path[i+1]][path]>0)
 {
 Connect[path[i+1]][path]=1;
 Connect[path][path[i+1]]=1;
 Weight[path][path[i+1]]=(-Weight[path[i+1]][path]);
 Length[path][path[i+1]]=(-Length[path[i+1]][path]);
 }else if(x[path[i+1]][path]==0)
 {
 Connect[path][path[i+1]]=0;
 Connect[path[i+1]][path]=1;
 }
 }
 }
}double NetWorkFlow:owestCost()
{
 if(x.size()!=0)
 {
 x.clear();
 }
 for(int i=0;i<(int)Connect.size();++i)
 {
 vector<double> d(Connect.size(),0);
 x.push_back(d);
 }
 vector<int> path=this->Find();
 int it=0;
 while(path.size()!=0)
 {
 for(int i=0;i<(int)path.size();++i)
 {
 cout<<path<<"\t";
 }
 cout<<endl;
 Augment(path);
 path=this->Find();
 ++it;
 }
 cout<<it<<endl;
 if(Solution.size()!=0)
 {
 Solution.clear();
 }
 for(int i=0;i<(int)a.size();++i)
 {
 vector<double> ss(b.size(),0);
 Solution.push_back(ss);
 }
 for(int i=0;i<(int)a.size();++i)
 {
 for(int j=0;j<(int)b.size();++j)
 {
 Solution[j]=x[i+1][1+a.size()+j];
 }
 }
 double result=0;
 for(int i=0;i<(int)a.size();++i)
 for(int j=0;j<(int)b.size();++j)
 {
 result+=Solution[j]*C[j];
 }
 return result;
}

lwd1981 · 发表于 2004-8-25 05:22:12

太长了啊！

liurenan · 发表于 2004-8-25 15:58:56

<

>此程序完全可以用matlab语言<

>那样会更简便的啊

knight · 发表于 2004-8-26 04:01:34

我想不一定会简便多少，Matlab的优势在于工具箱，这个算法似乎用到的工具箱不多。

		自动登录	找回密码
密码			注-册-帐-号