c语言最小费用流_用网络单纯形法(network simplex)解最小费用最大流(mincost)问题

该算法的原理来自线性规划，笔者还有没搞懂的地方，这里只谈实现。后面会给出完整的C++代码。

前提假设与问题点总结

事先有所了解的看官可以略过此小节。

如果一个节点是某些链路（或称有向边）的起点且不是任何链路的终点，则该节点为源节点(source)；如果一个节点是某些链路的终点且不是任何链路的起点，则该节点为宿节点(sink).

每条链路都可以有一定的流量(flow)，其方向必须与链路方向相同；一条链路上的流量允许的最大值称为该链路的容量(capacity)。

根据当前的流量情况，链路可分为三类:

满的（full）, 流量已达到最大值，即等于容量

空的（empty）, 流量为0

部分满的（partial）, 流量大于0但小于容量

最大流问题：对于一个存在源、宿节点的网络，选取其中一个源节点s和一个宿节点t，要求除了这两点外，其他节点上出、入的流量相等。安排各链路的流量，使s的流出、t的流入(显然，两者必相等)达到最大值。

最小费用最大流问题：假设每一单位的流量经过一条链路时都会产生一定的费用(cost), 已知各条链路的费用，给出解决方案，在达到上述最大流的同时，总费用最小。

这里约定流量、费用都是整数值。

解决方法概述

主要分以下步骤:

1. 构造一条虚拟的从s到t的链路，它的容量足够大，大于原本从s流出的所有链路的容量之和，也大于流入t的所有链路的容量之和；它的费用特别高，大于其他所有链路的费用的总和。使其他链路的流量都为0，虚拟链路流量等于其容量。

此时若在原有的链路上发现任意一条从s到t的路径、并将虚拟链路上的一部分流量导入到该条路径上，总费用必然减小。如此最大流和最小费用的问题就统一起来了：为了让费用达到最小，必定要尽量避开那条假设的链路，使更多的流量从原先的链路上经过。

以下图作为例子，cap和cost分别代表容量和费用，s = 0, t = 5：

s和t之间构造一个假设的链路，其容量和费用这里设为7和16，如下图所示：

初始阶段，只有虚拟链路是满的，流量为7，等于容量；其他链路流量均为0：

2. 构造一个以t节点为根的生成树。且该树必须包含虚拟链路。

比如对于刚才的例子构造生成树如下（具体的构造步骤等具体到代码时再谈，这里可以不用关心）

更直观的等效画法如下图所示。

注意加入生成树时不计链路方向，上两张图中的箭头只是为了说明树中节点间的父子关系。

这里再引入两个概念：

设对于一个节点u，都有一个值与其对应，称为该节点的potential（不知道标准中文翻译的词笔者就直接用原文了。后面也会有类似的情况), 记为φ(u).

设一从节点u到节点v的链路e，费用为C(e), 令 R(e) = C(e) - φ(u) + φ(v)， 称为链路e 的reduced cost.

做如下规定:

如果e在生成树上，则必有 R(e) = 0，由此推得 φ(v) = φ(u) - C(e).

可见只要给定生成树的根节点一个初值，按树的下沉方向即可确定所有节点的φ值.

3. 给定生成树的根节点一个任意的取值，由此推出所有节点的potential.

遍历所有不在生成树上的链路e，按 R(e) = C(e) - φ(u) + φ(v) 计算reduced cost, 若R(e) < 0，则做如下操作：

1). 找到u和v的最小共同祖先节点a, u到v、沿着树的方向从v到a的路径、沿着树的边从a到u的路径，三者组成了一个环路。

2). 沿着环路方向将尽量多的流推入该链路

3). 在环中选择一个空的或满的且离a最近的链路，将它移出生成树，而把e加入生成树。如果只有e是空的或满的，则保持树的结构不变；否则优先把e加入到生成树中；当树中有多个空的或满的链路时，优先把离a最近的链路移出。

“共同祖先节点”即既是u的祖先、又是v的祖先的节点。“最小共同祖先节点”即这些节点中离叶子最近的，简称LCA(Least Common Ancester)。

比如示例中的网络，遍历到链路1->4时，发现其reduced cost小于0, 分别从节点1和4向上追溯，发现其LCA为节点5，因而组成环路1->4->5->0->1。按1->4的方向推流，遍历1->4->5->0->1上的各个链路，如果链路方向与推流方向相同，则增加该链路的流量；如果相反，则减少流量。发现1->4、4->5 容量只有2，即最多推2个单位的流量。流量图现在变为:

接着将1->4移出生成树，将4->5加入到生成树。

注意即使环中无流量可推，重构生成树的步骤也不能省。

5. 重复步骤4，直到找不到任何reduce cost小于0且不在生成树中的链路。

代码:基本数据结构

完整代码可从这里查看

如下用一个vector代表网络中的所有节点，一个链表代表网络中的所有链接。

struct network {
  std::vector<vertex> vers; // 网络中的所有节点
  link_info *links; // 网络中的所有链表

  network(int V);
  virtual ~network();
  void copy_from(const network &other);
  network(const network &other);
  network &operator=(const network &other);
  int get_vercnt() const { return (int)vers.size(); }

  link_info *insert_link(int cap, int cost, int from, int to);
  void del_link(link_info *info);
protected:
  void insert_link_info(link_info *new_info);
  void del_link_info(link_info *to_del);

public:

  bool is_source(int i) const;
  bool is_sink(int i) const;
};

节点结构中只有一个链表, 表示它邻接的链路(即进入该节点的和从该节点出的所有链路)

struct vertex {
  // list of adjacent edges(links)
  link *adj_list;
  vertex() : adj_list(0) {}

  void insert_adj(link *new_link);
  void del_adj(link *to_del);
};

network::links为linkinfo类型的链表，每个linkinfo保存着链路的基本信息；vertex::adjlist为抽象类link的链表，它有两种派生类：link_forward和link_reverse. 每添加一个链路会在它起点的链表中增加一个实际类型为link_forward的节点，它终点的链表中增加一个实际类型为link_reverse的节点。link::info指向network::link中的一个节点，记录该节点的详细信息。

struct link;
// information of link as in a network
struct link_info {
  int cap; // capacity
  int flow; // flow 
  int cost; // cost per flow
  int from; // orderal number of source vertex in the network
  int to; // orderal number of destination vertex in the network

  link *forward;
  link *reverse;

  link_info *nxt, *prv;

  link_info(int c, int f, int cst, int src, int dst) :
    cap(c), flow(f), cost(cst), from(src), to(dst), 
    forward(0), reverse(0), nxt(0), prv(0)
  {}
};

struct link {
  link_info *info;
  link *prv, // the previous one in the adjacent list of a vertex
    *nxt; // the next one in the adjacent list of a vertex
  link(link_info *p) : info(p), prv(0), nxt(0)
  {
  }

  virtual int get_cap() const = 0;
  virtual int get_flow() const = 0;
  virtual int get_cost() const = 0;
  virtual int self_no() const = 0;
  virtual int peer_no() const = 0;
  virtual link *peer_link() const = 0;

  virtual int usable_flow() const = 0;
  virtual int use_flow(int d) = 0;

  virtual bool is_forward() const = 0;
  virtual bool is_reverse() const = 0;

  virtual ~link() {}
};

struct link_forward : public link {
  link_forward(link_info *p) : link(p)
  {
  }

  bool is_forward() const { return true; }
  bool is_reverse() const { return false; }

  int get_cap() const
  {
    return info->cap;
  }

  int get_flow() const
  {
    return info->flow;
  }

  int get_cost() const
  {
    return info->cost;
  }

  int self_no() const
  {
    return info->from;
  }

  int peer_no() const
  {
    return info->to;
  }

  link *peer_link() const
  {
    return info->reverse;
  }

  int usable_flow() const
  {
    return info->cap - info->flow;
  }

  int use_flow(int d)
  {
    info->flow += d;
    return usable_flow();
  }
};

struct link_reverse : public link {
  link_reverse(link_info *p) : link(p)
  {
  }

  bool is_forward() const { return false; }
  bool is_reverse() const { return true; }

  int get_cap() const
  {
    return info->cap;
  }

  int get_flow() const
  {
    return -info->flow;
  }

  int get_cost() const
  {
    return -info->cost;
  }

  int self_no() const
  {
    return info->to;
  }

  int peer_no() const
  {
    return info->from;
  }

  link *peer_link() const
  {
    return info->forward;
  }

  int usable_flow() const
  {
    return info->flow;
  }

  int use_flow(int d)
  {
    info->flow -= d;
    return usable_flow();
  }
};

代码：solution类

mincost_solution 表示最小流问题解决方案的基类，network_simplex为一个具体实现。

struct mincost_solution {
  virtual int solve(int s, int t) = 0;
  mincost_solution(network &ntw);

  virtual ~mincost_solution();

protected:
  std::vector<vertex> &vers;
  link_info *links;
  network &puzzle;
};

struct network_simplex : public mincost_solution {
  network_simplex(network &gra);

  int solve(int s, int t);

private:
  int phi_lv, lca_lv;
  int dummy_link_cap, dummy_link_cost;
  std::vector<mincost_verinfo> verinfos;


  link_info *prepare(int s, int t);
  void build_tree
  (link *cur_link, std::list<link *> &stk);
  int phi(int no);
  int lca(int y, int z);
  link *augment(link *to_add);
  bool on_path(int key, int from, int to) const;
  void reverse(int a, int dst);
  void update(link *to_del, link *to_add);
  int reduced_cost(link *p);
  void print_tree();
  void print_flow();
  void calc_flow_cost(int &cost);
};

代码：solve函数

把辅助函数写完了，主体的solve函数就可按前面说的流程展开了

int network_simplex::solve(int s, int t)
{
  int mincost = 0;
  if (!puzzle.is_source(s))
    throw std::runtime_error
      ("mincost : start vertex "
       "isn't a source");
  if (!puzzle.is_sink(t))
    throw std::runtime_error
      ("mincost : destination vertex "
       "isn't a sink");

  link_info *dummy_link_info = prepare(s, t);
  if (!dummy_link_info) {
    mincost = 0;
    return mincost;
  }
  print_flow();
  print_tree();
  int i, V = puzzle.get_vercnt();
  link *p;
  bool is_dec = true;
  for (phi_lv = 1; is_dec; phi_lv++) {
    is_dec = false;
    for (i = 0; i < V; i++) {
      if (!verinfos[i].related)
	continue;
      for (p = vers[i].adj_list; 
	   p != 0;
	   p = p->nxt) {
	if (!verinfos[p->peer_no()].related)
	  continue;
	int r = reduced_cost(p),
	  u = p->usable_flow();
	printf("%d->%d : r = %d, u = %dn", 
	       p->self_no(), p->peer_no(), r, u);
	if (r >= 0 || u <= 0)
	  continue;
	link *to_del = augment(p);
	print_flow();
	if (to_del != p)
	  update(to_del, p);
	print_tree();
	is_dec = true;
	break;
      }
      if (is_dec)
	break;
    }
  }

  puzzle.del_link(dummy_link_info);
  calc_flow_cost(mincost);
  return mincost;
}