Поток минимальной стоимости (min-cost-flow). Алгоритм увеличивающих путей

Дана сеть G, состоящая из N вершин и M рёбер. У каждого ребра (вообще говоря, ориентированному, но по этому поводу см. ниже) указана пропускная способность (целое неотрицательное число) и стоимость единицы потока вдоль этого ребра (некоторое целое число). В графе указан исток S и сток T. Даётся некоторая величина K потока, требуется найти поток этой величины, причём среди всех потоков этой величины выбрать поток с наименьшей стоимостью ("задача min-cost-flow").

Иногда задачу ставят немного по-другому: требуется найти максимальный поток наименьшей стоимости ("задача min-cost-max-flow").

Обе эти задачи достаточно эффективно решаются описанным ниже алгоритмом увеличивающих путей.

Описание

Алгоритм очень похож на алгоритм Эдмондса-Карпа вычисления максимального потока.

Простейший случай

Рассмотрим для начала простейший случай, когда граф - ориентированный, и между любой парой вершин не более одного ребра (если есть ребро (i,j), то ребра (j,i) быть не должно).

Пусть U_ij - пропускная способность ребра (i,j), если это ребро существует. Пусть C_ij - стоимость единицы потока вдоль ребра (i,j). Пусть F_ij - величина потока вдоль ребра (i,j), изначально все величины потоков равны нулю.

Модифицируем сеть следующим образом: для каждого ребра (i,j) добавим в сеть так называемое обратное ребро (j,i) с пропускной способностью U_ji = 0 и стоимостью C_ji = - C_ij. Поскольку, по нашему предположению, ребра (j,i) до этого в сети не было, то модифицированная таким образом сеть по-прежнему не будет мультиграфом. Кроме того, на всём протяжении работы алгоритма будем поддерживать верным условие: F_ji = - F_ij.

Определим остаточную сеть для некоторого зафиксированного потока F следующим образом (собственно, так же, как и в алгоритме Форда-Фалкерсона): остаточной сети принадлежат только ненасыщенные рёбра (т.е. у которых F_ij < U_ij), а остаточную пропускную способность каждого такого ребра как UPI_ij = U_ij - F_ij.

Собственно алгоритм min-cost-flow заключается в следующем. На каждой итерации алгоритма находим кратчайший путь в остаточной сети из S в T (кратчайший относительно стоимостей C_ij). Если путь не был найден, то алгоритм завершается, поток F - искомый. Если же путь был найден, то мы увеличиваем поток вдоль него настолько, насколько это возможно (т.е. проходим вдоль этого пути, находим минимальную остаточную пропускную способность MIN_UPI среди рёбер этого пути, и затем увеличиваем поток вдоль каждого ребра пути на величину MIN_UPI, не забывая уменьшать на такую же величину поток вдоль обратных рёбер). Если в какой-то момент величина потока достигла величины K (данной нам по условию величины потока), то мы также останавливаем алгоритм (следует учесть, что тогда на последней итерации алгоритма при увеличении потока вдоль пути нужно увеличивать поток на такую величину, чтобы итоговый поток не превзошёл K, но это выполнить легко).

Нетрудно заметить, что если положить K равным бесконечности, то алгоритм найдёт максимальный поток минимальной стоимости, т.е. один и тот же алгоритм без изменений решает обе задачи min-cost-flow и min-cost-max-flow.

Случай неориентированных графов, мультиграфов

Случай неориентированных графов и мультиграфов в концептуальном плане ничем не отличается от вышеописанного, поэтому собственно алгоритм будет работать и на таких графах. Однако возникают некоторые сложности в реализации, на которые следует обратить внимание.

Неориентированное ребро (i,j) - это фактически два ориентированных ребра (i,j) и (j,i) с одинаковыми пропускными способностями и стоимостями. Поскольку вышеописанный алгоритм min-cost-flow требует для каждого неориентированного ребра создать обратное ему ребро, то в итоге получается, что неориентированное ребро расщепляется на 4 ориентированных ребра, и мы фактически получаем случай мультиграфа.

Какие проблемы вызывают кратные рёбра? Во-первых, поток по каждому из кратных рёбер должен сохраняться отдельно. Во-вторых, при поиске кратчайшего пути нужно учитывать, что важно то, какое именно из кратных рёбер выбрать при восстановлении пути по предкам. Т.е. вместо обычного массива предков для каждой вершины мы должны хранить вершину-предка и номер ребра, по котором мы из неё пришли. В-третьих, при увеличении потока вдоль некоторого ребра нужно, согласно алгоритму, уменьшить поток вдоль обратного ребра. Поскольку у нас могут быть кратные рёбра, то придётся для каждого ребра хранить номер ребра, обратного ему.

Других сложностей с неориентированными графами и мультиграфами нет.

Анализ времени работы

По аналогии с анализом алгоритма Эдмондса-Карпа, мы получаем такую оценку: O (N M) * T (N, M), где T (N, M) - время, необходимое для нахождения кратчайшего пути в графе с N вершинами и M рёбрами. Если это реализовать с помощью простейшего варианта алгоритма Дейкстры, то для всего алгоритма min-cost-flow получится оценка O (N³ M), правда, алгоритм Дейкстры придётся модифицировать, чтобы он работал на графах с отрицательными весами (это называется алгоритм Дейкстры с потенциалами).

Вместо этого можно использовать алгоритм Левита, который, хотя и асимптотически намного хуже, но на практике работает очень быстро (примерно за то же время, что и алгоритм Дейкстры).

Реализация

Здесь приведена реализация алгоритма min-cost-flow, базирующаяся на алгоритме Левита.

На вход алгоритма подаётся сеть (неориентированный мультиграф) с N вершинами и M рёбрами, и K - величина потока, который нужно найти. Алгоритм находит поток величины K минимальной стоимости, если такой существует. Иначе он находит поток максимальной величины минимальной стоимости.

В программе есть специальная функция для добавления ориентированного ребра. Если нужно добавить неориентированное ребро, то эту функцию нужно вызывать для каждого ребра (i,j) дважды: от (i,j) и от (j,i).

const int INF = 1000*1000*1000;

struct rib {
	int b, u, c, f;
	size_t back;
};

void add_rib (vector < vector<rib> > & g, int a, int b, int u, int c) {
	rib r1 = { b, u, c, 0, g[b].size() };
	rib r2 = { a, 0, -c, 0, g[a].size() };
	g[a].push_back (r1);
	g[b].push_back (r2);
}

int main()
{
	int n, m, k;
	vector < vector<rib> > g (n);
	int s, t;
	... чтение графа ...

	int flow = 0,  cost = 0;
	while (flow < k) {
		vector<int> id (n, 0);
		vector<int> d (n, INF);
		vector<int> q (n);
		vector<int> p (n);
		vector<size_t> p_rib (n);
		int qh=0, qt=0;
		q[qt++] = s;
		d[s] = 0;
		while (qh != qt) {
			int v = q[qh++];
			id[v] = 2;
			if (qh == n)  qh = 0;
			for (size_t i=0; i<g[v].size(); ++i) {
				rib & r = g[v][i];
				if (r.f < r.u && d[v] + r.c < d[r.b]) {
					d[r.b] = d[v] + r.c;
					if (id[r.b] == 0) {
						q[qt++] = r.b;
						if (qt == n)  qt = 0;
					}
					else if (id[r.b] == 2) {
						if (--qh == -1)  qh = n-1;
						q[qh] = r.b;
					}
					id[r.b] = 1;
					p[r.b] = v;
					p_rib[r.b] = i;
				}
			}
		}
		if (d[t] == INF)  break;
		int addflow = k - flow;
		for (int v=t; v!=s; v=p[v]) {
			int pv = p[v];  size_t pr = p_rib[v];
			addflow = min (addflow, g[pv][pr].u - g[pv][pr].f);
		}
		for (int v=t; v!=s; v=p[v]) {
			int pv = p[v];  size_t pr = p_rib[v],  r = g[pv][pr].back;
			g[pv][pr].f += addflow;
			g[v][r].f -= addflow;
			cost += g[pv][pr].c * addflow;
		}
		flow += addflow;
	}

	... вывод результата ...

}