最短增益路径法求解最大流问题

1. 思路：

 首先我做了一些初始的规定：对于图中的点，规定第一个点（即下标为0）为起点，最后一个点为终点。用3个二维矩阵分别存储初始容量网络、流量网络和残余网络。由于课本上已经对最大流最小割定理进行详细的证明，这里就不再赘述。除了输入输出的操作外，我将整个算法分成3个步骤：

•  步骤1：根据容量网络C和当前流量网络F，计算残余网络R。
•  步骤2：在残余网络R中，利用广度优先搜索算法，计算出增广路径并返回，若没有找到增广路径则返回空。
•  步骤3：利用步骤2返回的增广路径，计算该路径上的流量，将该流量增加到流量网络F中。

 整体大致的思路如图1所示。

图1 整体算法思路

2. 输入：

 按照上面的思路，首先需要对图进行输入，这里我做了两种操作选择：

1.  1. 手动输入边以及权重，输入已知的图来验证程序的正确性。
2.  2. 随机生成边以及权重，用来计算时间复杂度。（随机生成的边数约为节点数的平方）

 输入的操作伪代码如下面所示。

伪代码

ManualInput(C):                // 手动输入边
while(TRUE):
    input i, j, w              // 分别输入边的起点，终点和权重
    if i < 0:                  // 由于起点不小于0，当输入i小于0时表示输入结束
        break
    C[i][j] = w                // 用邻接矩阵存储图
 
RandomInput(C):                // 随机生成边
fori=0 to C.length-1:        
    for j=0 to C.length-1:
        if i == j:             // 规定图没有自环，即权重为0
             C[i][j] = 0
        else:
             C[i][j] = rand()  // 用邻接矩阵存储图

3. 计算残余网络

 

 由于有了容量网络C和流量网络F，计算残余网络的操作非常简单。首先用容量网络C减去流量网络F，接着增加流量的反向边，就能得到残余网络R。 

伪代码
calculateR(C,F, R):
fori=0 to len-1:
    for j=0 to len-1:
        R[i][j] = C[i][j] - F[i][j]   // 容量减去流量
        if F[i][j] != 0:
R[j][i] += F[i][j]                    // 增加流量的反向边

4. 计算增广路径（广度优先搜索）

 由于在我的实验中，图的边的输入是完全随机的，导致整个容量图是比较稀疏或是比较稠密也是随机的，因此这里用深度优先搜索算法或是广度优先搜索算法没有太大差别，我选择了用广度优先搜索。

4.1 增广路径的定义

 首先遇到的问题是，如何将找到的增广路径返回，也就是路径该如何定义。这里我结合广度优先搜索的特点，定义图的节点结构，用来存储返回的路径。

 由于广度优先搜索后会生成广度优先搜索树，因此可以将该树记录下来，从终止节点开始寻找父亲节点，直到找到起始节点，就能找到路径。定义节点的结构，对每一个节点都有两个属性，一个是自己的标记index，另一个是它的父亲节点parent。

 假设经过广度优先搜索后，生成的节点列表如图2、3所示，由于节点0是起始节点，节点5是终止节点，因此我们从节点5开始，找到它的父亲节点3，再从节点3找到它的父亲节点4……这样下去，直到找到节点0，整个路径0-> 4 -> 3 -> 5就找到了。

图2 经过广度优先搜索后生成的节点列表

图3 根据广度优先搜索树找到增广路径

4.2 广度优先搜索

 根据广度优先搜索也上面的节点列表的结构，需要一个队列queue，和一个记录节点是否已经被访问过的数组。首先访问节点0也就是起始节点，将其加入队列，然后进入循环直到队列为空。在循环中，每次取出队首元素，遍历并访问其所有未访问过的邻接节点，将它们加入队列，同时将他们的父亲节点设为自己。

 由于这个步骤是为了找到增广路径，因此无需搜索所有的节点：当访问到终止节点时，说明已经找到从起始节点到终止节点的路径，即可停止搜索返回节点列表。若经过广度优先搜索没有访问到终止节点，说明起始节点到终止节点之间没有通路，因此返回空值。算法的流程如下面的伪代码所示。

伪代码
BFS():
visitnodeList[0]                                 // 访问起始节点
queue.push(nodeList[0])                          // 起始节点入队
whilequeue is not empty:                         // 循环至队列为空
    node = queue.pop()                           // 取出队首
    index = node.index
    for j=1 to len-1:                            // 遍历所有邻接节点
        if R[index][j] == 0: continue            // 没有邻接
        if nodeList[j] is visited: continue      // 已被访问
        nodeList[j].parent = node                // 将邻接节点的父亲设为自己
        queue.push(nodeList[j])                  // 邻接节点入队
        visitnodeList[j]                         // 访问邻接节点
        if j == N - 1:                           // 访问到终止节点，返回列表
             return nodeList
return NULL                                      // 没有访问到终止节点，返回空

5. 利用增广路径增加流

 经过上面的计算得到了记录了增广路径的节点列表，需要将这个路径提取出来，计算增广路径上能增加的流，即路径上所有边权值中的最小值。因此，需要一个结构来记录边，将路径转化为边的集合，计算其中权值的最小值cfp，并将流增加到流量网络中。由于权值在残余网络R中已经记录，因此这里的边只需记录流入会让流出的节点即可。增加流的步骤分为两步，首先需要计算cfp的值，再将cfp值增加到流量网络F中，因此需要一个栈，在计算cfp值的同时生成边压入栈，计算出cfp值后将栈中所有边在F中的权值加上cfp。

5.1 计算cfp值并生成边的栈

 首先将cfp值设为无穷大，从终止节点开始寻找父亲节点，生成边并将边压入栈，若该边的权值小于cfp，则更新cfp值，直到遍历到起始节点为止。

 

伪代码：
cfp= MAX                                 // cfp初始为无穷大
idx= len-1                               // 从终止节点开始
whilenodeList[idx].parent is not NULL:   // 循环至父亲节点为空，即初始节点
    node = nodeList[idx]                  
    edge.to = node.index                  // 边的流出为节点本身
    edge.from = node.parent               // 边的流入为父亲节点
    stack.push(edge)                      // 将边压入栈
    if R[edge.from][edge.to] < cfp:       // 若边的权值小于cfp，则更新cfp
        cfp = R[edge.from][edge.to]

5.2 将cfp添加到流量网络中

 这个过程非常简单，仅需将前面压入栈的边全部取出，依次在流量网络F中对应的边加上cfp即可。

伪代码：
while stack is not empty:
	edge = stack.pop()                // 取出边
	F[edge.from][edge.to] += cfp      // 在F中加上cfp 

6. 测试算法正确性

 为了检验整个算法的正确性，这里我做了一个简单的图，如图4所示，也就是初始的容量网络C，通过手动输入的方式来观察整个算法的执行过程和结果。

图4 初始的容量网络C

 执行程序，手动输入图4的所有边及权重，得到的结果如图5所示。可以看到，算法执行过程中，一共找到了4次增广路径，增加的流量cfp分别为4、2、2、2，直到最后找不到增广路径时，算法结束，得到的最大流为10。为了验证结果的正确性，我将自己手动计算增广路径，来验证该结果与手动计算的结果是否一致。

图5 算法执行结果

 如图6所示，是我按照广度优先搜索的顺序，找到的增广路径和增加的流量网络，直到在残余网络中找不到增广路径为止。将图中的结果与前面运行的结果比较，发现找到的增广路径、增加的流量以及最终得到的最大流都完全一致，说明算法正确。

图6 手动计算过程

7. 测试算法的性能

 对于前面实现的算法，我将测试节点规模为500到1000，随机生成边的条数接近节点数的平方时，算法运行的时间，得到的结果如表1所示。将表格转化成图表，如图7所示。可以看到，随着规模的增加，程序运行时间呈曲线增长。

 

表1 算法消耗的时间

规模	500	600	700	800	900	1000
时间(s)	2.541	4.28	6.541	9.941	13.718	18.595

图7 表1转化成图表

8. 算法的优化

 在自己手动计算的时候，可以发现新的残余网络其实可以通过增广路径在旧的残余网络上加减得到，而初始的残余网络就是容量网络C，因此可以考虑取消残余网络R，直接将残余网络的计算结果存到C中，这样能减少部分计算残余网络的时间，又能减少R所消耗的空间。

 具体的思路如图8所示，对于图中旧的残余网络R1，找到的增广路径为0 -> 1 -> 2 -> 5，该路径上的流量为4。以0 -> 1的边为例，在旧的残余网络中0 -> 1边的容量为6，没有反向边，也就是反向边1 -> 0的容量为0，由于增广路径经过0 -> 1边，因此只需将0 -> 1边的容量减去流量4得到2，其反向边1 -> 0的容量增加流量4得到4。其他经过的边1 -> 2、2 -> 5也进行相同的处理，得到新的残余网络R2。

图8 在旧的残余网络上修改得到新的残余网络

 因此，整个算法只需在原来增加流量的部分，增加对残余网络的操作，而将原来计算残余网络R的操作去掉。由于路径最多的边数为|V|-1，因此原来计算残余网络操作消耗的O(|V|²)降到了O(|V|)，总体上能减少许多时间。对算法主要的修改如下面伪代码所示，仅在原来的基础上增加两行代码。

伪代码：
whilestack is not empty:
    edge = stack.pop()                // 取出边
    F[edge.from][edge.to] += cfp      // 在F中加上cfp
 C[edge.from][edge.to]-= cfp      // 正向边减去cfp
 C[edge.to][edge.from]+= cfp      // 反向边加上cfp

 将修改后的程序运行，记录下运行时间后跟修改前的表格比较，比较结果如图9所示。可以看到，优化后的程序速度为优化前的5倍左右，说明对程序的优化是成功的。

图9 优化前后的时间比较 

9. 大规模测试

 最后我对程序进行大规模的测试，分别测试了5000和10000的节点，边的条数接近点数的平方，得到的结果如图10、11所示，可以看到，当节点规模为5000，边数接近25000000时，程序运行时间为7.5分钟；当节点规模为10000，边数接近1亿时，程序运行时间接近1个小时。

图10 规模为5000时的运行结果

图11 规模为10000时的运行结果

10. 代码

#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <queue>
#include <stack>

#define MAX_WEIGHT 1000

using namespace std;

int maxFlow = 0;

// 节点结构
struct Node {
    int idx;
    Node* parent;
};
// 边结构
struct Edge {
    int from;
    int to;
};

void createCF(int **C, int **F, int N);                // 新建容量和流网络
void deleteCF(int **C, int **F, int N);                // 删除容量和流网络
void printCF(int **C, int N);                          // 打印容量和流网络
void inputNet(int **C);                                // 手动输入容量网络，格式为【流出节点 流入节点 容量】，以-1 -1 -1 结束
void randomInput(int **C, int N);                      // 随机输入容量网络
Node* bfs(int **C, int N);                             // 用广度优先搜索，找出增广路径并返回路径，没有路径则返回NULL
void updateCF(int **C, int **F, int N, Node* result);  // 利用返回的增广路径，将路径流量增加到F，同时计算残余网络

int main() {
    srand((unsigned)time(NULL));
    clock_t start, finish;
    int N;
    int **C, **F;
    cout << "输入节点个数：";
    cin >> N;

    C = new int*[N];
    F = new int*[N];
    createCF(C, F, N);

    randomInput(C, N);
    // inputNet(C);

    start = clock();
    while (1) {
        Node* nodeList = bfs(C, N);
        if (nodeList == NULL) break;                   // 若找不到增广路径，说明已达到最大流
        updateCF(C, F, N, nodeList);
    }
    finish = clock();

    cout << "最大流：" << maxFlow << endl;
    double duration = double(finish - start) / CLOCKS_PER_SEC;
    cout << "需要的时间为：" << duration << "s" << endl;
    deleteCF(C, F, N);
	return 0;
}

Node* bfs(int **C, int N) {
    Node* nodeList = new Node[N];
    int* visited = new int[N];
    for (int i=0; i<N; ++i) {
        nodeList[i].idx = i;
        nodeList[i].parent = NULL;
        visited[i] = 0;
    }
    queue<Node> que;
    que.push(nodeList[0]);
    visited[0] = 1;
    while(!que.empty()) {
        Node node = que.front();
        int idx = node.idx;
        que.pop();
        for (int j=1; j<N; j++) {
            if (C[idx][j] == 0) continue;
            if (visited[j]) continue;
            visited[j] = 1;
            nodeList[j].parent = &nodeList[idx];
            que.push(nodeList[j]);
            if (j == N-1) {
                delete [] visited;
                return nodeList;
            }
        }
    }
    delete [] nodeList;
    delete [] visited;
    return NULL;
}

void updateCF(int **C, int **F, int N, Node* result) {
    stack<Edge> path;
    int cfp = MAX_WEIGHT;
    int idx = N-1;
    while(result[idx].parent) {
        Edge eg;
        eg.to = idx;
        eg.from = result[idx].parent->idx;
        if (C[eg.from][eg.to] == 0)
            cout << "should be zero" << endl;
        if (C[eg.from][eg.to] < cfp)
            cfp = C[eg.from][eg.to];
        path.push(eg);
        idx = result[idx].parent->idx;
    }
    while (!path.empty()) {
        Edge eg = path.top();
        F[eg.from][eg.to] += cfp;
        C[eg.from][eg.to] -= cfp;
        C[eg.to][eg.from] += cfp;
        path.pop();
    }
    maxFlow += cfp;
    delete [] result;

}

void randomInput(int **C, int N) {
    for (int i=0; i<N; ++i) {
        for (int j=0; j<N; ++j) {
            if (i == j) {
                C[i][j] = 0;
                continue;
            }
            C[i][j] = abs(rand()) % MAX_WEIGHT;
        }
    }
}

void createCF(int **C, int **F, int N) {               // 新建容量和流网络
    for (int i=0; i<N; ++i) {
        C[i] = new int[N];
        F[i] = new int[N];
    }
    for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j) {
        C[i][j] = 0;
        F[i][j] = 0;
    }
}

void deleteCF(int **C, int **F, int N) {               // 删除容量和流网络
    for (int i=0; i<N; ++i) {
        delete [] C[i];
        delete [] F[i];
    }
    delete [] C;
    delete [] F;
}

void inputNet(int **C) {                               // 手动输入容量网络
    int i, j, w;
    cout << "输入加权矩阵：" << endl;
    while(1) {
        cin >> i >> j >> w;
        if (i==-1)
            break;
        C[i][j] = w;
    }
}

void printCF(int **CF, int N) {                        // 打印容量和流网络
    for (int i=0; i<N; ++i) {
        for (int j=0; j<N; ++j) {
            cout << CF[i][j] << " ";
        }
        cout << endl;
    }
}