远航休息栈
首页 > 笔记 > 最小费用最大流算法总结

最小费用最大流算法总结

考前水图论

我们已经学过了求最大流的算法，现在我们想象假如我们有一个流量网络，现在每个边除了流量，现在还有一个单位费用，这条边的费用相当于它的单位费用乘上它的流量，我们要保持最大流的同时，还要保持边权最小，这就是最小费用最大流问题。

我们该如何求最小费用最大流呢？

其实是可以用最短路更新的。

我们每次要更新的时候就可以在费用网络上跑最短路。因为有负权的问题，我们可以先用spfa跑。

核心代码如下

while(!que.empty())
    {
        int now=que.front();que.pop();vis[now]=0;
        for (int i=st[now];~i;i=e[i].next)
            if (e[i].flow>0 && dis[e[i].to]>dis[now]+e[i].cost)
            {
                dis[e[i].to]=dis[now]+e[i].cost;
                pe[e[i].to]=i,pv[e[i].to]=now;
                if (!vis[e[i].to])
                    vis[e[i].to]=1,que.push(e[i].to);
            }
    }

这时spfa更新的条件就是：还有流量可以流且现在所连到的点的dis（费用）比现在自己的dis加上当前边的费用要大。（相当于spfa更新边权）同时还要更新一下能增广的最大流量。

最后的时候我们需要更新一下

核心代码：

while(spfa(S,T))
    {
        flow=INF;
        for (int i=T;i!=S;i=pv[i])
            flow=min(flow,e[pe[i]].flow);
        COST+=flow*dis[T];
        FLOW+=flow;
        for (int i=T;i!=S;i=pv[i])
            e[pe[i]].flow-=flow,e[pe[i]^1].flow+=flow;
    }

因为我们求出了最短路，实际上就是求出了这条增广路上的单位权值的和。那么总增广的费用就是最短路*总流量

在上一段代码中我们记录了一个p的数组，这样的话就可以通过汇点反着找到这条最短路，接着增广流量。

 

以为这种方法在稠密图较慢，但是它可以处理负权，而且可以方便的改成最大费用最大流。

代码

#include <cstdio>
#include <algorithm>
#include <queue>
#include <cstring>
#define M 2001000
#define N 1001000
#define INF 0x33333333
#define min(x,y) ((x<y)?(x):(y))
using namespace std;
typedef pair<int,int> Pair;
struct node{int from,to,next,flow,cost;}e[M];
int tot=-1,st[M];
int n,m,x,y,z;
int pe[N],pv[N],dis[N],vis[N];
void add(int x,int y,int z,int zz){
    e[++tot].to=y;
    e[tot].from=x;
    e[tot].flow=z;
    e[tot].cost=zz;
    e[tot].next=st[x];
    st[x]=tot;
}
void add_edge(int x,int y,int z,int zz){add(x,y,z,zz),add(y,x,0,-zz);}
queue<int>que;
bool spfa(int S,int T)
{
    memset(dis,0x33,sizeof dis);
    memset(vis,0,sizeof vis);
    que.push(S),vis[S]=1,dis[S]=0;
    while(!que.empty())
    {
        int now=que.front();que.pop();vis[now]=0;
        for (int i=st[now];~i;i=e[i].next)
            if (e[i].flow>0 && dis[e[i].to]>dis[now]+e[i].cost)
            {
                dis[e[i].to]=dis[now]+e[i].cost;
                pe[e[i].to]=i,pv[e[i].to]=now;
                if (!vis[e[i].to])
                    vis[e[i].to]=1,que.push(e[i].to);
            }
    }
    return dis[T]<INF;
}
Pair mfmc(int S,int T)
{
    int COST=0,FLOW=0,flow;
    while(spfa(S,T))
    {
        flow=INF;
        for (int i=T;i!=S;i=pv[i])
            flow=min(flow,e[pe[i]].flow);
        COST+=flow*dis[T];
        FLOW+=flow;
        for (int i=T;i!=S;i=pv[i])
            e[pe[i]].flow-=flow,e[pe[i]^1].flow+=flow;
    }
    return make_pair(FLOW,COST);
}
main()
{
    int m,from,to;
    scanf("%d%d%d%d",&n,&m,&from,&to);
    memset(e,-1,sizeof e);
    memset(st,-1,sizeof st);
    int S=from,T=to;
    for(int i=0;i<m;i++)
    {
        int u,v,flow,cost;
        scanf("%d%d%d%d",&u,&v,&flow,&cost);
        add_edge(u,v,flow,cost);
    }
    Pair ans=mfmc(S,T);
    printf("%d %d\n",ans.first,ans.second);
}

 

接下来就是重点，用Dijkstra求解费用流

有人可能会说,Dijkstra只能在正权图上求最短路呀,而费用流的网络构图中\(w[x][y] = -w[y][x]\),势必会产生负权边,怎么能用Dijkstra呢？(下文用\(w[i][j]\)表示费用,\(p[i]\)表示点,方便讲解)

我们求解费用流一般用spfa或zkw流,这两种方法在稀疏图上运行效率还不错,但是在稠密图上的表现却无法令人满意,而且其运行效率与图的形态关系太大,无法令人安心。

反观Dijkstra,它在稠密图\(O(n^2)\)和稀疏图\(O(mlogm)\)上的表现都不错,而且其运行时间与图的形态关系不大,较为稳定,唯一的缺点是不能运行在负权图上。

事实上在对网络进行一些处理后,Dijkstra是可以用来求解费用流的。

我们对网络G中的每个点i设置一个势函数\(h[i]\),在算法运行的每一时刻,对于每条残余网络中的边\((x, y)\),都要满足\(h[x] + w[x][y] – h[y] >= 0\)(三角不等式)。维护h[i]的方法下文再讨论。

构建新网络\(G’,V’ = V,w'[x][y] = h[x] + w[x][y] – h[y]\),这样在新网络中所有边的权值均非负,可以使用Dijkstra求解最短路。

\(G’\)中的一条路径 \(p[1], p[2],\cdots, p[m]\) 的权值为

\begin{align*}
w'[path]&=w'[p[1]][p[2]] + w'[p[2]][p[3]] + \cdots + w'[p[m-1]][p[m]]\\
&=h[p[1]]-h[p[2]] + w[p[1]][p[2]]\\
&+h[p[2]]-h[p[3]] + w[p[2]][p[3]]\\
&+\cdots\\
&+h[p[m-1]]-h[p[m]] + w[p[m-1]][p[m]]\\
&=w[p[1]][p[2]] + w[p[2]][p[3]] + \cdots + w[p[m-1]][p[m]] + h[p[1]]-h[p[m]] \\
&=w[path] + h[p[1]]-h[p[m]]
\end{align*}

而当p[1]和p[m]确定时,最优化\(w[path]\)与最优化\(w'[path]\)是等价的,所以我们可以在G’中求解S到T的最短路\(w'[path]\),G中S到T的最短路\(w[path] = w'[path] – h[S] + h[T]\)。

下面我们来讨论如何维护势函数\(h[]\)。

初始情况下,如果所有的权值都为正,那么可以简单地将所有\(h[i]\)设置为0；如果有负权值,那么我们做一遍spfa,让h[]等于距离函数。如果初始网络为DAG,也可以采用递推\(h[i] = min(h[j] + w[j][i], g[j][i] = true)\)。

让\(d'[i]\)表示G’中S到点i的距离,当某次增广结束后,G’中会新加入某些边\((j, i)\),而这些(j, i)必定满足\(d'[i] + w'[i][j] = d'[j]\)(否则(i, j)就不会在增广路中)。对上式进行一些代数变换：

\begin{align*}
d'[i] + w'[i][j] = d'[j]
\\d'[i] + w[i][j] + h[i]-h[j] = d'[j]
\\(d'[j] + h[j])-(d'[i] + h[i])-w[i][j] = 0
\\(d'[j] + h[j])-(d'[i] + h[i]) + w[j][i] = 0
\end{align*}
(因为是费用流,所以有w[i][j] = -w[j][i])

因此让所有\(h[i] += d'[i]\)后,新加入的边(j, i)也会满足势函数的性质。

同时我们有：

\begin{align*}
d'[i] + w'[i][j] – d'[j] >= 0\\
d'[i] + h[i] – h[j] + w[i][j] – d'[j] >= 0\\
(d'[i] + h[i]) – (d'[j] + h[j]) + w[i][j] >= 0\\
\end{align*}

因此修改\(h[i]\)后,\((i, j)\)依然会满足势函数的性质。

算法过程如下：

>S1 初始化h[]
>S2 在残留网络中做Dijkstra
>S3 若S到T有可行路径,则修改增广路上的边的容量并所有h[i] += d'[i],转S2,否则退出

算法时间复杂度为：\(O(spfa() + K * Dijkstra())\)或\(O(K * Dijkstra())\),这取决于初始化h[]时是否调用spfa,K表示增广次数。

代码如下

#include <cstdio>
#include <algorithm>
#include <queue>
#include <cstring>
#define M 2001000
#define N 1001000
#define INF 0x33333333
#define min(x,y) ((x<y)?(x):(y))
using namespace std;
typedef pair<int,int> Pair;
struct node
{
    int from,to,next,flow,cost;
}e[M];
int tot=-1,st[M];
int n,m,x,y,z;
void add(int x,int y,int z,int zz)
{
    e[++tot].to=y;
    e[tot].from=x;
    e[tot].flow=z;
    e[tot].cost=zz;
    e[tot].next=st[x];
    st[x]=tot;
}
Pair main_pro(int s,int t)
{
    static int h[N];
    int flow=0,cost=0;
    while(1)
    {
        static int dis[N],pv[N],pe[N];
        memset(dis,0x33,sizeof dis);    //距离初始化为inf
        dis[s]=0;
        priority_queue<Pair,vector<Pair>,greater<Pair> >que;
        que.push(Pair(0,s));
        while(!que.empty())
        {
            Pair now=que.top();
            que.pop();
            if (now.first!=dis[now.second]) continue;
            //dis[now.second]是目前的值，now.first是入堆时的值，如果目前的优，直接continue
            if (now.second==t) break;
            for (int i=st[now.second];~i;i=e[i].next)
            {
                int nowcost=e[i].cost+h[now.second]-h[e[i].to];
                //现在的费用为这条路径的费用+起始点的势函数-结尾点的势函数
                if (e[i].flow>0 && dis[e[i].to]>dis[now.second]+nowcost) //松弛
                {
                    dis[e[i].to]=dis[now.second]+nowcost;
                    que.push(Pair(dis[e[i].to],e[i].to));
                    pv[e[i].to]=now.second; //更新节点所对应的路径
                    pe[e[i].to]=i;
                }
            }
        }
        if (dis[t]==INF)  //没有增广到T
            break;
        for (int i=0;i<n;i++)  //更新势函数
            h[i]=min(h[i]+dis[i],INF);
        int newflow=INF;
        for (int x=t;x!=s;x=pv[x])  //找到最小流量
            newflow=min(newflow,e[pe[x]].flow);
        flow+=newflow;
        cost+=newflow*h[t];
        for (int x=t;x!=s;x=pv[x])  //更新流量
            e[pe[x]].flow-=newflow,e[pe[x]^1].flow+=newflow;
    }
    return make_pair(flow,cost);
}

int main()
{
    int m,from,to;
    scanf("%d%d%d%d",&n,&m,&from,&to);
    memset(e,-1,sizeof e);
    memset(st,-1,sizeof st);
    for(int i=0;i<m;i++)
    {
        int u,v,flow,cost;
        scanf("%d%d%d%d",&u,&v,&flow,&cost);
        add(u,v,flow,cost);
        add(v,u,0,-cost);
    }
    Pair ans=main_pro(from,to);
    printf("%d %d\n",ans.first,ans.second);
}