网络流之最大流问题

نُشر في 2023-02-11 عُدل في 2023-02-13

网络流之最大流问题 - 徐子睿

1 算法的引入

UNIX插座

有 n 个插座，m 个设备和 k (n,m,k≤100) 种转换器，每种转换器有无限多。已知每个插座的类型，每个设备的插头类型，以及每种转换器的插座类型和插头类型。插头和插座类型都用不超过 24 个字母表示，插座只能插到类型名称相同的插座中。例如，有 4 个插座，类型分别为 A，B，C，D；有 5 个设备，插头类型分别为 B，C，B，B，X；还有三种转换器，分别是 B->X，X->A 和 X->D。这里用 B -> X 表示插座类型为 B，插头类型为 X，因此一个插座类型为 B 的设备插上这种转换器之后就 “变成” 了一个插头类型为 X 的设备。转换器可以级联使用，例如插头类型为 A 的设备依次接上 A->B, B->C, C->D 这 3 个转换器之后会 “变成” 插头类型为 D 的设备。要求插的设备尽量多。问最少剩几个不匹配的设备。

题目分析

题目中提到每种转换器都有无限多个，所以我们可以预处理出每个插头可以转换出的插头类型。

这部分的预处理，考虑到 n,m,k≤100 ，即使所有的插头、设备、转换器都不同最多也只可能有400种不同的名字（而这种情况显然是无解的），我们可以用Floyd算法算出转换器之间的传递闭包来处理转换器之间的可转换关系

接下来的问题就是如何分配各个插头（可能是转换后的）

而这就需要我们下面的模型——网络流。它就是专门来处理分配问题的模型

2.1 网络流

在一个有向图上选择一个源点，一个汇点，每一条边上都有一个流量上限（以下称为容量），即经过这条边的流量不能超过这个上界，同时，除源点和汇点外，所有点的入流和出流都相等，而源点只有流出的流，汇点只有汇入的流。这样的图叫做网络流。

所谓网络或容量网络指的是一个连通的赋权有向图 D = （V、E、C），其中V 是该图的顶点集，E是有向边(即弧)集，C是弧上的容量。此外顶点集中包括一个起点和一个终点。网络上的流就是由起点流向终点的可行流，这是定义在网络上的非负函数，它一方面受到容量的限制，另一方面除去起点和终点以外，在所有中途点要求保持流入量和流出量是平衡的。（引自百度百科）

上图中1为源点，6为汇点

为了下面讲解方便，我们引入如下几个概念：

源点$S$：只有流出去的点
汇点$T$：只有流进来的点
流量$Flow$：一条边上流过的流量
容量$Capacity$：一条边上可供流过的最大流量
残量$Rest$：一条边上的容量 - 流量

网络流图中有这样三个重要的性质值得我们注意：

对任意一条边，总有 $Flow < Capacity$
对任意一个非S/T点u总有 $\sum^{}_{p \in V} flow[p][u] \equiv \sum^{}_{v \in V} flow[u][v]$ (即这个点的入流与出流相等)
对任意一个有向边(u, v)总有$flow[u][v] \equiv -flow[v][u]$

2.2 最大流问题

最大流问题就是求一个网络中S-T流量的最大值。

2.2.1Edmonds-Karp算法

Edmonds-Karp算法基于增广路定理，下面我们介绍该定理：

我们计算每条边上的残量，得到原网络的残量网络

残量网络中任意一条s-t通路都能够对应原图中的一条增广路。而这条道路上的最小残量就是这条通路上可以增大的流量。

自然的，若残量网络中不存在增广路时，原网络即达到最大流

基于这个定理，我们便可以设计出相应的算法算出网络中的最大流，即Edmonds-Karp算法。

在EK算法中我们在加入正向边的同时，又加入了一条反向的、capacity为0的反向边，这样做的原因是如果我们在增广的时候找到的增广路并不是最优解的话，反向边的存在可以使得从其他路流经这条边时，可以通过流过他的反向边来撤销把这条边作为增广路的决策。

我们给出一个例子：

若是一开始我们增广的决策是 u -> v

而实际上的最优解是u -> q;p -> v;

那么我们就可以通过增广到p时走u -> v的反向边来撤销u -> v上原有的流。

当然，这种撤销也可以是不完全的，如果我们的最优解需要u,p同时对v输入流量的话，也可以只撤销一部分的u -> v流。

下面给出完整的代码并配有细节讲解：

class EdmondsKarp{
    private:
        
    vector <edge> edges;           //存储边
    vector <int> g[maxn];          //存储以某节点为起点的边序号
    const int inf = 0x3f3f3f3f;     
    int pre[maxn], adv[maxn];      //pre存储某节点增广路上的入弧
                                   //adv存储由s到该点的最小余量（可改进量）
    public:

    void addedge(int u, int v, int c){
        edges.push_back((edge){u, v, c, 0});
        edges.push_back((edge){v, u, 0, 0});
        int i = edges.size();
        g[u].push_back(i - 2);      //正向边
        g[v].push_back(i - 1);      //反向边
    }

    int maxFlow(int s, int t){      //BFS搜索增广路
        int res = 0;
        while(1){
            memset(adv, 0, sizeof(adv));
            memset(pre, 0, sizeof(pre));
            queue <int> q;
            q.push(s);
            adv[s] = inf;
            while(!q.empty()){
                int u = q.front(); q.pop();
                for(auto i : g[u]){
                        edge &e = edges[i];         //用引用减少代码量，增加可读性
                    if(!adv[e.v] && e.c > e.f){     
                        //若capacity > flow则说明这条边在残量网络中可以是增广路，将其记录
                        //adv可以同时充当vis数组的作用
                        adv[e.v] = min(adv[e.u], e.c - e.f); //计算可改进量
                        pre[e.v] = i;               //记录入弧
                        q.push(e.v);
                    }
                }
                if(adv[t]) break;   //一旦搜索到t点，则代表增广路已经找到
            }
            if(!adv[t]) break;      //若已经没有增广路，则已找到最大流
            else{
                for(int v = t;v != s;v = edges[pre[v]].u){
                    edge &e = edges[pre[v]], &einv = edges[pre[v] ^ 1];
                    e.f += adv[v]; 
                    einv.f -= adv[v];
                }
                res += adv[t];
            }
        }
        return res;
    }

};

时间复杂度分析：

每增广一次，层次网络中必定有一条边会被删除。层次网络中最多有m条边，所以认为最多可以增广m次。在最短增广路算法中，用BFS来增广，一次增广的复杂度为$O(n+m)$，其中$O(m)$为BFS的花费，$O(n)$为修改流量的花费。所以在每一阶段寻找增广路的复杂度为$O(m(m+n)) = O(m^2)$。因此n个阶段寻找增广路的算法总复杂度为$O(nm^2)$。

EK算法虽然可以计算出网络中的最大流，但其有一些有待改进的地方：

在搜索过程中，一旦搜索到t点的增广路就退出循环，浪费了对未搜索到的增广路的计算。有大量的计算冗余。

下面我们来看一个例子来更好地理解这一点

我们有一个图，如图一。

按照套路，我们先BFS，找S−T最短路。所有的距离标号都画在了图二上（EK算法可能用不到，但Dinic用得到）。

假设我们选的是S−3−T这条路，增广（如图三，绿色）

然后我们再来一遍BFS。。。等等！

细心的你可能也发现了，S−1−T也是一条S−T最短路。

那就增广吧！（如图四）

可以检查一下，这时候没有长度为2的最短路了。

但EK算法不会这样。它会再笨拙地BFS一遍，这就浪费了不少时间。

所以说，多路增广是很重要的。

故我们需要找到一个新的算法，让其在一次循环内多找出几条增广路，避免己经计算出的数据的浪费。

2.2.2 Dinic算法

在Dinic算法中，我们引入了分层图这一概念，更具体地，我们先计算出各节点到s点的无权最短路并将其记录为该点的层次，且最短路上的边必须可以增广，这一步相当于把网络中的节点分成不同的层。
而当t点的层次不存在时则意味着没有从s到t的增广路，此时最大流已经找到

与EK算法不同的是，当我们用BFS计算出网络的分层图后，对于增广路的搜索我们用的是DFS。DFS时，可以继续搜索的条件就是v的层次比u大1.

为什么只能大1呢，回到我们BFS时，我们在分层图中计算的是各节点到s点的无权最短路，所以层次差超过2的一对点（u， v）直接相连的情况是不存在的，因为如果这样的点对真的存在，那么在距离为1的u，v中，u、v到s的最短路长度竟然差了不止1，这样显然是错误的。

下面给出代码并进行细节讲解：

class Dinic{
    private:
    
    vector <edge> edges;
    vector <int> g[maxn];
    int level[maxn];                //存储节点层次
    int res = 0;
    const int inf = 0x3f3f3f3f;

    public:

    void addedge(int u, int v, int c){
        edges.push_back((edge){u, v, c, 0});
        edges.push_back((edge){v, u, 0, 0});
        int i = edges.size();
        g[u].push_back(i - 2);
        g[v].push_back(i - 1);
    }

    bool BFS(int s, int t){
        memset(level, 0, sizeof(level));
        queue <int> q;
        q.push(s);
        level[s] = 1;
        while(!q.empty()){
            int u = q.front(); q.pop();
            for(auto i : g[u]){
                edge &e = edges[i];
                if(e.c > e.f && level[e.v] == 0){
                    //只有增广路存在时才能分配层次
                    level[e.v] = level[e.u] + 1;
                    q.push(e.v);
                }
            }
        }
        if(!level[t]) return false; //若t点层次不存在则说明没有s->t增广路，返回false退出循环
        else return true;
    } 

    int DFS(int u, int adv, int t){
        if(u == t) return adv; //搜索到t点代表已经找到增广路，返回可改进量
        for(auto i : g[u]){
            edge &e = edges[i], &einv = edges[i ^ 1];
            if(level[e.v] == level[e.u] + 1 && e.c > e.f){
                int res = DFS(e.v, min(adv, e.c - e.f), t);
                if(res > 0){
                    e.f += res;
                    einv.f -= res;
                    return adv;
                }
            }
        }
        return 0;
    }

    int maxFlow(int s, int t){
        int res = 0;
        while(BFS(s, t)){                       //只要t节点的层次存在则还有增广路
            while(int adv = DFS(s, inf, t)){    //DFS搜索出增广路返回总可改进量
                res += adv;
            }
        }
        return res;
    }
};

复杂度分析：

Dinic算法最多被分为n个阶段，这是因为每次建立层次图后t的层次是严格递增的。

现在来分析DFS过程的总复杂度。在每一阶段，将DFS分成两部分分析。

修改增广路的流量并后退的花费。在每一阶段，最多增广m次，每次修改流量的费用为$O(n)$。而一次增广后在增广路中后退的费用也为$O(n)$。所以在每一阶段中，修改增广路以及后退的复杂度为$O(nm)$。
DFS遍历时的前进与后退。在DFS遍历时，如果当前路径的最后一个顶点能够继续扩展，则一定是沿着第i层的顶点指向第i+1层顶点的边向汇点前进了一步。因为增广路经长度最长为n，所以最坏的情况下前进n步就会遇到汇点。在前进的过程中，可能会遇到没有边能够沿着继续前进的情况，这时将路径中的最后一个点在层次图中删除。

注意到每后退一次必定会删除一个顶点，所以后退的次数最多为n次。在每一阶段中，后退的复杂度为$O(n)$。

假设在最坏的情况下，所有的点最后均被退了回来，一共共后退了n次，这也就意味着，有n次的前进被“无情”地退了回来，这n次前进操作都没有起到“寻找增广路”的作用。除去这n次前进和n次后退，其余的前进都对最后找到增广路做了贡献。增广路最多找到m次。每次最多前进n个点。所以所有前进操作最多为n+m*n次，复杂度为$O(nm)$。

于是得到，在每一阶段中，DFS遍历时前进与后退的花费为$O(nm)$。

综合以上两点，一次DFS的复杂度为$O(nm)$。因此，Dinic算法的总复杂度即$O(n ^ 2 m)$。

当然，Dinic算法可以再被优化：

注意到当我们在一次DFS是如果搜索到了一条增广路，我们要回退到s点再重新寻找增广路，而这时我们仍可能搜索到之前已经增广过的路线，所以我们引入当前弧这一概念，用一个数组cur记录点u之前循环到了哪一条边，以此来加速避免重复的寻找同一条边。

优化后代码如下：

class Dinic{
    private:
    
    vector <edge> edges;
    vector <int> g[maxn];
    int level[maxn], cur[maxn];
    int res = 0, n = 0;
    const int inf = 0x3f3f3f3f;

    public:

    void addedge(int u, int v, int c){
        edges.push_back((edge){u, v, c, 0});
        edges.push_back((edge){v, u, 0, 0});
        int i = edges.size();
        g[u].push_back(i - 2);
        g[v].push_back(i - 1);
    }

    inline void setN(int N){
        n = N;
    }

    bool BFS(int s, int t){
        memset(level, 0, sizeof(level));
        queue <int> q;
        q.push(s);
        level[s] = 1;
        while(!q.empty()){
            int u = q.front(); q.pop();
            for(auto i : g[u]){
                edge &e = edges[i];
                if(e.c > e.f && level[e.v] == 0){
                    level[e.v] = level[e.u] + 1;
                    q.push(e.v);
                }
            }
        }
        if(!level[t]) return false;
        else return true;
    } 

    int DFSAdv(int u, int adv, int t){
        if(u == t) return adv;
        for(int &i = cur[u];i < g[u].size();i++){  //引用使得i变化时同时改变u点的当前弧
            edge &e = edges[g[u][i]], &einv = edges[g[u][i] ^ 1];
            if(level[e.v] == level[e.u] + 1 && e.c > e.f){
                int res = DFS(e.v, min(adv, e.c - e.f), t);
                if(res > 0){
                    e.f += res;
                    einv.f -= res;
                    return res;
                }
            }
        }
        return 0;
    }

    int maxFlowAdv(int s, int t){
        int res = 0;
        while(BFS(s, t)){
            memset(cur, 0, sizeof(cur)); // 在新层次图上搜索时要重置当前弧
            while(int adv = DFSAdv(s, inf, t)){
                res += adv;
            }
        }
        return res;
    }
};

3 引入题的解答

有了上面算法的介绍，我们便可以设计出UNIX插头问题的代码了

首先用Floyd计算出转换器的传递闭包
接下来构造网络，令第i个设备的插头为plug[i]，第i个插座为in[i]。
构造s点，对于每一个设备都连接一条s到plug[i]，capacity为1的边；构造t点，对于每一个插座都连接一条in[i]到t，capacity为1的边；对于每个plug[i]，连接一条capacity为inf的边到其可以转换的插座上
以上网络的最大流即为最多的匹配插头，用总数减去最大流即为最小剩余插头。

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;

int KASE;
const int maxn = 1000;
const int inf = 0x3f3f3f3f;


struct edge{
    int u, v, c, f;
};

class Dinic{
    private:
    
    vector <edge> edges;
    vector <int> g[maxn];
    int level[maxn], cur[maxn];
    int res = 0, n = 0;
    const int inf = 0x3f3f3f3f;

    public:

    void addedge(int u, int v, int c){
        edges.push_back((edge){u, v, c, 0});
        edges.push_back((edge){v, u, 0, 0});
        int i = edges.size();
        g[u].push_back(i - 2);
        g[v].push_back(i - 1);
    }

    bool BFS(int s, int t){
        memset(level, 0, sizeof(level));
        queue <int> q;
        q.push(s);
        level[s] = 1;
        while(!q.empty()){
            int u = q.front(); q.pop();
            for(auto i : g[u]){
                edge &e = edges[i];
                if(e.c > e.f && level[e.v] == 0){
                    level[e.v] = level[e.u] + 1;
                    q.push(e.v);
                }
            }
        }
        if(!level[t]) return false;
        else return true;
    } 

    int DFS(int u, int adv, int t){
        if(u == t) return adv;
        for(auto i : g[u]){
            edge &e = edges[i], &einv = edges[i ^ 1];
            if(level[e.v] == level[e.u] + 1 && e.c > e.f){
                int res = DFS(e.v, min(adv, e.c - e.f), t);
                if(res > 0){
                    e.f += res;
                    einv.f -= res;
                    return res;
                }
            }
        }
        return 0;
    }

    int maxFlow(int s, int t){
        int res = 0;
        while(BFS(s, t)){
            while(int adv = DFS(s, inf, t)){
                res += adv;
            }
        }
        return res;
    }
};

int d[maxn][maxn];

void solve(){
    int plugCnt, inCnt, convertCnt;
    vector <string> name;
    vector <int> plugs, ins;
    memset(d, 0, sizeof(d));
    map <string, int> idx;
    cin>>inCnt;
    for(int i = 0;i < inCnt;i++){
        string s;
        cin>>s;
        if(!idx[s]){
            name.push_back(s);
            idx[s] = name.size() + 1;
        }
        ins.push_back(idx[s]);
    }
    cin>>plugCnt;
    for(int i = 0;i < plugCnt;i++){
        string temp, s;
        cin>>temp>>s;
        if(!idx[s]){
            name.push_back(s);
            idx[s] = name.size() + 1;
        }
        plugs.push_back(idx[s]);
    }
    cin>>convertCnt;
    for(int i = 0;i < convertCnt;i++){
        string u, v;
        cin>>u>>v;
        if(!idx[u]){
            name.push_back(u);
            idx[u] = name.size() + 1;
        }
        if(!idx[v]){
            name.push_back(v);
            idx[v] = name.size() + 1;
        }
        d[idx[u]][idx[v]] = 1;
    }
    for(int k = 2;k <= name.size() + 1;k++) d[k][k] = 1;
    for(int k = 2;k <= name.size() + 1;k++){
        for(int i = 2;i <= name.size() + 1;i++){
            for(int j = 2;j <= name.size() + 1;j++){
                d[i][j] = d[i][j] || (d[i][k] && d[k][j]);
            }
        }
    }
    int S = name.size();
    int s = 1, t = S + 2;
    Dinic D;
    for(auto i : plugs){
        D.addedge(s, i, 1);
        for(auto j : ins){
            if(d[i][j]){
                D.addedge(i, j + S + 1, inf);
            }
        }
    }
    for(auto i :ins){
        D.addedge(i + S + 1, t, 1);
    }
    cout<<plugs.size() - D.maxFlow(s, t);
}

int main(){
    ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    cin>>KASE;
    while(KASE){
        solve();
        cout<<"\n";
        if(--KASE) cout<<"\n";
    }
    return 0;
}

4 例题

通过这道题我们可以更好地理解网络流是如何解决分配问题的。

UVA11082 矩阵解压 Matrix Decompressing

已知一矩阵的行数$R$与列数$C$，以及前$i$行的前缀和$a_i$、前$j$列的前缀和$b_j$，矩阵元素的取值为$1-20$求一合法矩阵。

题目分析

由前缀和我们可以反向得到每一行、列的元素和。

值得注意的是，每一列的元素都分别属于不同的行，更抽象的说，这一列的元素在列上汇聚到了一起又分配给了各行，如下图所示

这启示了我们也许可以用网络流来进行元素的分配

还有一个小问题，网络流中允许的最小流量是0，但题目要求最少是1

对于这个问题，我们可以把每个colsum 减去 R，rowsum 减去 C。这样处理之后元素的取值范围就变成了0-19，可以构造网络流。

那么本题的思路就已经基本出来了

计算每行、列的元素和$rowsum[R]$,$colsum[C]$（元素取值范围为0-19）;
构建网络，构造源点$s$，对于每个行索引$i$连接一条$s$到$i$、$capacity$为$colsum[i]$的边；构造汇点t，对于每个列索引j连接一条$j$到$t$、$capacity$为$rowsum[j]$的边；对于每个$i$，对所有$j$连接一条$i$到$j$，$capacity$为$19$的边
计算网络的最大流，如果最大流等于所有元素之和则有解，且第$i$行第$j$列元素即为从$i$到$j$边上的$flow$。

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;

int KASE, T;
const int maxn = 100;

struct edge{
    int u, v, c, f;
};

class Dinic{
    public:
    
    vector <edge> edges;
    vector <int> g[maxn];
    int level[maxn] = {0}, cur[maxn] = {0};
    int res = 0, n = 0;
    const int inf = 0x3f3f3f3f;

    

    void addedge(int u, int v, int c){
        edges.push_back((edge){u, v, c, 0});
        edges.push_back((edge){v, u, 0, 0});
        int i = edges.size();
        g[u].push_back(i - 2);
        g[v].push_back(i - 1);
    }

    bool BFS(int s, int t){
        memset(level, 0, sizeof(level));
        queue <int> q;
        q.push(s);
        level[s] = 1;
        while(!q.empty()){
            int u = q.front(); q.pop();
            for(auto i : g[u]){
                edge &e = edges[i];
                if(e.c > e.f && level[e.v] == 0){
                    level[e.v] = level[e.u] + 1;
                    q.push(e.v);
                }
            }
        }
        if(!level[t]) return false;
        else return true;
    } 

    int DFS(int u, int adv, int t){
        if(u == t) return adv;
        for(auto i : g[u]){
            edge &e = edges[i], &einv = edges[i ^ 1];
            if(level[e.v] == level[e.u] + 1 && e.c > e.f){
                int res = DFS(e.v, min(adv, e.c - e.f), t);
                if(res > 0){
                    e.f += res;
                    einv.f -= res;
                    return res;
                }
            }
        }
        return 0;
    }

    int maxFlow(int s, int t){
        int res = 0;
        while(BFS(s, t)){
            while(int adv = DFS(s, inf, t)){
                res += adv;
            }
        }
        return res;
    }
};

void solve(){
    int R, C;
    cin>>R>>C;
    int rowsum[maxn] = {0}, colsum[maxn] = {0};
    for(int i = 1;i <= R;i++) cin>>rowsum[i];
    for(int i = R;i >= 1;i--) rowsum[i] -= (rowsum[i - 1] + C);
    for(int i = 1;i <= C;i++) cin>>colsum[i];
    for(int i = C;i >= 1;i--) colsum[i] -= (colsum[i - 1] + R);
    Dinic D;
    int s = R + C + 1, t = R + C + 2;
    for(int i = 1;i <= R;i++){
        D.addedge(s, i, rowsum[i]);
    }
    for(int i = 1;i <= C;i++){
        D.addedge(i + R, t, colsum[i]);
    }
    for(int i = 1;i <= R;i++){
        for(int j = 1;j <= C;j++){
            D.addedge(i, j + R, 19);
        }
    }
    D.maxFlow(s, t);
    cout<<"Matrix "<<T - KASE + 1<<"\n";
    for(int i = 1;i <= R;i++){
        for(int j = 1;j <= C;j++){
            for(auto k : D.g[i]){
                if(D.edges[k].v == j + R){
                    cout<<D.edges[k].f + 1<<" ";
                    break;
                }
            }
        }
        cout<<"\n";
    }
}

int main(){
    ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    cin>>KASE;
    T = KASE;
    while(KASE){
        solve();
        if(KASE--) cout<<"\n";
    }
    return 0;
}