Сетевой симплекс: симплекс-метод, который ходит по деревьям

Если ты знаешь, почему базис задачи о потоке - остовное дерево, как считают потенциалы и как устроен оборот по циклу, листай до «Теперь сам».

Сначала разберём транспортную задачу без верхних ограничений на отдельные маршруты: в ней небазисные рёбра действительно имеют нулевой поток. Затем добавим конечные ёмкости. Программа ниже реализует первый вариант; конечные ёмкости она не принимает.

Как думаешь, как выглядит угол многогранника для потока

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

Прежде чем читать дальше, ответь себе: в базисном оптимальном плане перевозок без верхних ёмкостей шесть складов и восемь магазинов. Сколько рёбер из сорока восьми будут нести ненулевой поток? Подсказка: если ненулевые рёбра образуют цикл, поток вокруг цикла можно чуть увеличить в одну сторону или уменьшить в другую, и одно из двух не дороже.

На что это похоже

Перестройка маршрутов по одному. Логист держит действующую схему - по каким дорогам что возят - и для каждого склада и магазина знает «внутреннюю цену» товара в этом месте. Он смотрит на дорогу, которой сейчас не пользуются: если по ней товар приедет дешевле внутренней цены, её стоит открыть. Открыв, он перекидывает часть перевозок по образовавшемуся кольцу, пока какая-то из старых дорог не опустеет, и закрывает её. Схема снова древовидная, цены пересчитаны, ищем следующую выгодную дорогу.

Механизм: дерево, потенциалы, оборот

базис: остовное дерево, поток идёт только по его рёбрам
       (n вершин → n − 1 рёбер; рёбра вне дерева - ноль)

потенциалы: π(корня) = 0; на ребре дерева u → v: π(v) = π(u) + c(u,v)
приведённая цена ребра:  r(u,v) = c(u,v) + π(u) − π(v)
       на рёбрах дерева r = 0; вне дерева - выгода открыть ребро

шаг (оборот):
    вход: ребро вне дерева с r < 0  (нет таких - план оптимален)
    цикл: это ребро + путь в дереве между его концами
    по направлению u → v рёбра цикла «попутные» (+θ) или «встречные» (−θ)
    θ = наименьший поток на встречных рёбрах
    выход: встречное ребро, ставшее нулём

Почему базис - дерево в этом варианте. Если ненулевые рёбра содержат цикл, поток по нему можно сдвинуть на ±ε; цена меняется на ±ε · (цена цикла), и сдвиг в одну из сторон не хуже. Сдвигаем, пока ребро цикла не обнулится: циклов становится меньше, цена не растёт. В итоге ненулевые рёбра - лес, и до дерева его дополняют рёбрами с нулевым потоком.

Как получить стартовое дерево. Добавляют вершину-корень и от каждого склада к корню (от корня к каждому магазину) искусственное ребро очень большой цены - все запасы идут через корень. Симплекс сам вытеснит дорогие искусственные рёбра; если какое-то осталось с потоком, задача невыполнима.

«Большую цену» нельзя выбирать наугад. Общий надёжный способ - две фазы: сначала реальным дугам назначить цену 0, искусственным 1 и минимизировать суммарный искусственный поток. Положительный минимум означает невыполнимость. При нуле фиксируем потоки искусственных дуг нулём и возвращаем реальные цены. Нельзя оставить их доступными с ценой 1: они могли бы снова стать выгоднее настоящих маршрутов. Если нулевая искусственная дуга ещё нужна для описания дерева, её можно временно держать в базисе, сохраняя нулевую ёмкость. В целочисленной программе ниже применяется одна фаза с M=10⁶: весь реальный поток не больше 240 единиц, цены не больше 30, стоимость допустимого реального плана не больше 7200. Любой положительный целый искусственный поток дороже этой границы.

Цены r - те же потенциалы Джонсона из прошлого урока, только посчитанные по дереву, а не кратчайшими путями. Отрицательная приведённая цена у ребра вне дерева значит, что цикл через него имеет отрицательную стоимость: сетевой симплекс - это отмена отрицательных циклов, где циклы ищутся не Беллманом-Фордом, а по дереву.

Правило

В варианте без верхних ёмкостей сетевой симплекс держит остовное дерево с потоком и потенциалы, при которых приведённая цена c(u,v) + π(u) − π(v) на рёбрах дерева нулевая. Ребро вне дерева с отрицательной приведённой ценой входит в дерево, замыкая цикл; по циклу пускают θ - наименьший поток на встречных рёбрах, - и обнулившееся встречное ребро выходит. Нет рёбер с отрицательной ценой - план оптимален. Стартовое дерево дают дорогие искусственные рёбра через корень.

Разбор: два склада, два магазина

склады A 5, B 5;  магазины X 4, Y 6
цены: A→X 2  A→Y 4  B→X 3  B→Y 1

дерево: A→Y 5, B→X 4, B→Y 1       стоимость 20 + 12 + 1 = 33
потенциалы: π(A) = 0, π(Y) = 4, π(B) = 4 − 1 = 3, π(X) = 3 + 3 = 6
вне дерева: r(A→X) = 2 + 0 − 6 = −4  < 0 → входит

цикл: A→X (попутное), X…B по B→X (встречное), B→Y (попутное),
      Y…A по A→Y (встречное)
θ = min(B→X 4, A→Y 5) = 4
новый поток: A→X 4, B→X 0 (выходит), B→Y 5, A→Y 1
стоимость 8 + 5 + 4 = 17 - на 4 · 4 = 16 меньше

проверка: π(A) 0, π(X) 2, π(Y) 4, π(B) 3
          r(B→X) = 3 + 3 − 2 = 4 ≥ 0 → оптимум

Один оборот снизил стоимость ровно на |r| · θ = 4 · 4. Так и должно быть: цикл с приведённой ценой −4 прошли четыре единицы.

Разбор: конечные ёмкости меняют правило

Когда у рёбер есть верхние ёмкости

Теперь зададим 0≤f(e)≤u(e). В базисе свободные рёбра, чей поток может изменяться в обе стороны, не содержат цикла; их дополняют до дерева T. Рёбра вне дерева делятся на L с потоком 0 и U с потоком u(e). Ненулевые рёбра теперь могут содержать цикл: насыщенная небазисная дуга не обязана принадлежать дереву.

На дереве по-прежнему r(e)=c(e)+π(u)−π(v)=0. Для e из L можно только увеличивать поток, поэтому улучшение возможно при r(e)<0. Для e из U можно только уменьшать его: остаточное направление обратное, с ценой −r(e). Улучшение возможно при r(e)>0. Если нарушений нет, все остаточные направления имеют неотрицательную приведённую цену, отрицательных циклов нет и план оптимален.

После выбора направления его добавляют к дереву и меняют поток по циклу. Допустимая величина шага - минимум из u(e)−f(e) на попутных рёбрах и f(e) на встречных. Дуга, достигшая границы, может выйти из дерева. Если ограничивающей оказалась сама входящая дуга, дерево остаётся прежним, а дуга переходит между L и U. Это перестановка границы, а не обязательный обмен двух рёбер дерева.

Например, есть две параллельные дуги A→B и спрос на 5 единиц. Дорогая дуга стоит 5, имеет ёмкость 10 и входит в дерево с потоком 3. Дешёвая стоит 1, имеет ёмкость 2 и насыщена потоком 2 вне дерева. При π(A)=0,π(B)=5 её цена r=−4. Это уже оптимум: дешёвую дугу увеличить нельзя, её обратная остаточная цена +4. Правило «каждое отрицательное r вне дерева требует входа» здесь было бы неверно.

Другой старт: все 5 единиц на дорогой дуге, дешёвая пуста. Тогда по циклу можно перенести две единицы: запас дешёвой дуги 2 ограничивает шаг раньше, чем опустеет дорогая. Дешёвая насыщается и остаётся небазисной, стоимость падает с 25 до 17. Это проверка, что границы нужны и в выборе направления, и в ограничении шага.

Разбор: шесть складов, восемь магазинов

Программа решает транспортную задачу сетевым симплексом со стартом через искусственный корень и сверяет стоимость с алгоритмом дешёвых путей из прошлого урока.

import random

rng = random.Random(22)
sklady, magaziny = 6, 8
n = sklady + magaziny                 # вершины 0..n−1, корень R = n
zapas = [rng.randint(10, 40) for _ in range(sklady)]
spros = [0] * magaziny
for _ in range(sum(zapas)):           # спрос в сумме равен запасу
    spros[rng.randrange(magaziny)] += 1
b = zapas + [-x for x in spros] + [0]
dugi = [(i, sklady + j, rng.randint(1, 30))       # (откуда, куда, цена)
        for i in range(sklady) for j in range(magaziny)]


def setevoj_simpleks():
    R, M = n, 10 ** 6                 # искусственные дуги - дорогие
    vse = list(dugi)
    potok = {}
    for v in range(n):                # старт: всё через корень
        d = (v, R, M) if b[v] >= 0 else (R, v, M)
        vse.append(d)
        potok[len(vse) - 1] = abs(b[v])
    derevo = set(potok)
    oborotov = 0
    while True:
        sos = {x: [] for x in range(n + 1)}
        for k in derevo:
            u, v, _ = vse[k]
            sos[u].append((v, k))
            sos[v].append((u, k))
        pi, otec, glub = {R: 0}, {R: None}, {R: 0}
        stek = [R]
        while stek:                   # потенциалы: на дереве r = 0
            x = stek.pop()
            for y, k in sos[x]:
                if y not in pi:
                    u, v, c = vse[k]
                    pi[y] = pi[x] + c if u == x else pi[x] - c
                    otec[y], glub[y] = (x, k), glub[x] + 1
                    stek.append(y)
        r = lambda k: vse[k][2] + pi[vse[k][0]] - pi[vse[k][1]]
        vhod = min((k for k in range(len(vse)) if k not in derevo),
                   key=r)
        u, v, c = vse[vhod]
        if r(vhod) >= 0:              # отрицательных цен нет - оптимум
            break
        vpered, nazad = [vhod], []    # цикл: дуга u→v и путь v … u
        a, z = u, v
        while a != z:                 # поднимаемся к общему предку
            if glub[a] >= glub[z]:
                x, k = otec[a]
                (vpered if vse[k][1] == a else nazad).append(k)
                a = x
            else:
                x, k = otec[z]
                (vpered if vse[k][0] == z else nazad).append(k)
                z = x
        teta = min(potok[k] for k in nazad)
        for k in vpered:
            potok[k] = potok.get(k, 0) + teta
        for k in nazad:
            potok[k] -= teta
        uhod = next(k for k in nazad if potok[k] == 0)
        derevo.add(vhod)
        derevo.remove(uhod)
        oborotov += 1
    assert all(potok.get(k, 0) == 0 for k in range(len(dugi), len(vse)))
    return sum(potok.get(k, 0) * vse[k][2] for k in range(len(dugi))), \
        oborotov


def ssp():                            # проверка: дешёвые пути по одному
    INF = float("inf")
    S, T = n, n + 1
    ost, cena = {}, {}
    def dobavit(u, v, cap, c):
        ost[u, v] = ost.get((u, v), 0) + cap
        ost.setdefault((v, u), 0)
        cena[u, v], cena[v, u] = c, -c
    for u, v, c in dugi:
        dobavit(u, v, INF, c)
    for i in range(sklady):
        dobavit(S, i, zapas[i], 0)
    for j in range(magaziny):
        dobavit(sklady + j, T, spros[j], 0)
    itog = shagov = 0
    while True:
        d, ot = {S: 0}, {}
        for _ in range(n + 2):        # Беллман-Форд по остаточной сети
            for (u, v), cap in ost.items():
                if cap > 0 and u in d:
                    if d[u] + cena[u, v] < d.get(v, INF):
                        d[v], ot[v] = d[u] + cena[u, v], u
        if T not in d:
            return itog, shagov
        put, v = [], T
        while v != S:
            put.append((ot[v], v))
            v = ot[v]
        f = min(ost[e] for e in put)
        for u, v in put:
            ost[u, v] -= f
            ost[v, u] += f
        itog, shagov = itog + f * d[T], shagov + 1


a, ob = setevoj_simpleks()
c2, sh = ssp()
print("сетевой симплекс: стоимость", a, " оборотов", ob)
print("дешёвые пути:     стоимость", c2, " путей", sh)

сетевой симплекс: стоимость 997  оборотов 14
дешёвые пути:     стоимость 997  путей 14

Стоимости совпали. Четырнадцать оборотов: каждый - один обход дерева на пятнадцать вершин и просмотр рёбер вне дерева, без поиска путей по всей сети. В промышленных решателях на сетях в миллионы рёбер результат зависит от структуры сети, правила выбора рёбер и реализации. Гарантии полиномиального времени у простейших правил выбора входящего ребра нет.

Что ломается без этого

Первая ошибка - считать потенциалы не по дереву, а кратчайшими путями: на рёбрах дерева приведённая цена должна быть ровно нулём, иначе оборот выбирает не тот цикл. Вторая - перепутать попутные и встречные рёбра цикла: уменьшать надо встречные, и θ - минимум именно по ним; перепутав, получишь отрицательный поток. Третья - вырожденные обороты с θ = 0: дерево меняется, а план - нет, и простое правило выбора может ходить по кругу бесконечно; промышленные реализации держат «сильно допустимые» деревья, где такой круг невозможен. Четвёртая - забыть, что искусственные рёбра с потоком в конце - это ответ «план невыполним», а не часть плана.

Зачем это в работе

Транспортные задачи и планирование поставок, распределение мощностей между заводами и рынками, задачи о назначениях большого размера, расписания экипажей - всё, что сводится к потоку минимальной стоимости, промышленные решатели решают сетевым симплексом. Его же ценят за «тёплый старт»: если изменились цены нескольких дорог, оборачивать можно от прежнего оптимального дерева, и новый ответ находится за несколько оборотов. А потенциалы в конце - это двойственные цены: сколько стоит лишняя единица товара в каждом пункте.

Теперь сам

В разборе после оборота дерево - A→X, A→Y, B→Y. Склад B получил ещё одну единицу запаса, а магазин X - ещё одну единицу спроса. На сколько вырастет оптимальная стоимость и как это видно по потенциалам?