Целочисленная программа: записать задачу переменными 0/1

Если ты запишешь условия «если - то», «не больше одного», «хотя бы один» и «обслуживать только открытым» линейными неравенствами на переменные 0/1 и объяснишь, чем сильная запись лучше слабой, листай до «Теперь сам».

Как думаешь, что здесь нужно программировать

Шесть проектов, бюджет 13. У каждого польза и цена. Проект D имеет смысл только вместе с A; B и C используют одну команду, поэтому вместе нельзя; из E и F нужен хотя бы один - этого требует регулятор. Выбрать набор с наибольшей пользой.

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

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

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

Искусство - в заполнении. Одну и ту же заявку можно вписать так, что автомат решит её за миллисекунду, и так, что он будет думать час. В этом уроке разница выглядит как 49 против 61.

Механизм: шаблоны и релаксация

Переменная x_k ∈ {0, 1} - взят ли объект k. Ограничения и цель - линейные: сумма переменных с коэффициентами сравнивается с числом.

УСЛОВИЕ                              НЕРАВЕНСТВО
бюджет                               5·x_A + 4·x_B + ... <= 13
если D, то A                         x_D <= x_A
не больше одного из B, C             x_B + x_C <= 1
хотя бы один из E, F                 x_E + x_F >= 1
ровно один склад на клиента j        x_1j + x_2j + x_3j = 1
клиента j обслуживает только         x_ij <= y_i           (сильно)
  открытый склад i                   Σ_j x_ij <= n · y_i   (слабо)
и A, и B (новая переменная z)        z <= x_A,  z <= x_B,
                                     z >= x_A + x_B − 1

Последняя строка - как убрать произведение x_A · x_B, которое в линейную программу не пишется: новая переменная z зажата неравенствами так, что равна единице, только когда обе равны единице.

ЛП-релаксация. Заменим x ∈ {0, 1} на 0 ≤ x ≤ 1 - получится линейная программа, которая решается быстро (раздел о потоках и ЛП). Её допустимых точек больше, значит, её оптимум не хуже целого: для минимума - нижняя граница, для максимума - верхняя. Это то самое «ослабить задачу, чтобы получить границу» из главы о ветвях и границах.

Сильная и слабая запись. «Склад обслуживает клиента, только если открыт» можно записать одним неравенством на склад - Σ_j x_ij ≤ n · y_i - или неравенством на каждую пару - x_ij ≤ y_i. Для целых переменных это одно и то же. Для дробных нет: в слабой записи y_i = 0,2 разрешает складу обслуживать 0,2 · 5 = 1 клиента целиком, в сильной - каждого клиента не больше чем на 0,2. Сильная релаксация ближе к целому оптимуму, а чем ближе граница, тем больше ветвей отсекает решатель.

Правило

Целочисленная программа: переменные 0/1, линейные ограничения, линейная цель. Логические условия переводятся в неравенства по шаблонам, произведение переменных - через новую переменную с тремя неравенствами. Без требования целочисленности получается ЛП-релаксация - граница оптимума. Из равносильных записей выбирают ту, чья релаксация ближе к целому оптимуму: она сильнее.

Разбор: шесть проектов и четыре условия

Решатель - milp из библиотеки scipy, внутри у него HiGHS, открытый решатель целочисленных программ. Он минимизирует, поэтому пользу берём со знаком минус. Все ограничения записаны в виде «не больше»: «хотя бы один из E, F» стало −x_E − x_F ≤ −1. Эта программа запускается в обычном Python с установленным scipy, а не в движке курса.

import numpy as np
from scipy.optimize import milp, linprog, LinearConstraint, Bounds

imena = ["A", "B", "C", "D", "E", "F"]
polza = [12, 9, 10, 8, 5, 7]
cena = [5, 4, 5, 3, 3, 4]
A = np.array([cena,
              [-1, 0, 0, 1, 0, 0],
              [0, 1, 1, 0, 0, 0],
              [0, 0, 0, 0, -1, -1]])
verh = np.array([13, 0, 1, -1])
c = -np.array(polza)
celoe = milp(c, integrality=np.ones(6), bounds=Bounds(0, 1),
             constraints=[LinearConstraint(A, -np.inf, verh)])
vybrano = []
for k in range(6):
    if celoe.x[k] > 0.5:
        vybrano.append(imena[k])
print(vybrano, -round(celoe.fun, 2))
lp = linprog(c, A_ub=A, b_ub=verh, bounds=(0, 1))
print([abs(round(float(v), 2)) for v in lp.x], -round(lp.fun, 2))

['A', 'B', 'F'] 28.0
[1.0, 0.5, 0.0, 1.0, 1.0, 0.0] 29.5

Ответ - A, B, F: цена 5 + 4 + 4 = 13, польза 28. Условия соблюдены: D не взят, B без C, F закрывает требование регулятора.

Вторая строка - релаксация: A, D, E целиком и половина B, польза 29,5. Половину проекта не запустишь, но число полезно: лучше 29,5 целого ответа нет, значит, найденные 28 отстают от оптимума не больше чем на полтора. Именно так решатель и рассуждает: дробная точка даёт границу, а дробная переменная x_B = 0,5 подсказывает, по чему ветвиться.

Разбор: склады, слабая запись против сильной

Три возможных склада с ценой открытия 30, 40, 35 и пять клиентов с ценой доставки от каждого склада. y_i - открыт ли склад, x_ij - обслуживает ли склад i клиента j. Задача - размещение предприятий: открыть часть складов и прикрепить каждого клиента к открытому так, чтобы сумма открытий и доставок была наименьшей. Две записи связи «только открытый», для каждой - релаксация и целый оптимум. LinearConstraint(A, lb, ub) в milp - это lb ≤ A·x ≤ ub для каждой строки A; -np.inf означает «снизу не ограничено».

import numpy as np
from scipy.optimize import milp, linprog, LinearConstraint, Bounds

otkryt = [30, 40, 35]
dostavka = [[2, 9, 8, 6, 10],
            [7, 3, 4, 9, 5],
            [9, 8, 3, 2, 4]]
m = len(otkryt)
n = len(dostavka[0])
peremennyh = m + m * n

def x(i, j):
    return m + i * n + j

c = np.zeros(peremennyh)
for i in range(m):
    c[i] = otkryt[i]
    for j in range(n):
        c[x(i, j)] = dostavka[i][j]

kazhdyj = np.zeros((n, peremennyh))
for j in range(n):
    for i in range(m):
        kazhdyj[j, x(i, j)] = 1

slabaya = np.zeros((m, peremennyh))
for i in range(m):
    slabaya[i, i] = -n
    for j in range(n):
        slabaya[i, x(i, j)] = 1

silnaya = np.zeros((m * n, peremennyh))
for i in range(m):
    for j in range(n):
        silnaya[i * n + j, i] = -1
        silnaya[i * n + j, x(i, j)] = 1

for imya, svyaz in [("слабая", slabaya), ("сильная", silnaya)]:
    lp = linprog(c, A_ub=svyaz, b_ub=np.zeros(len(svyaz)),
                 A_eq=kazhdyj, b_eq=np.ones(n), bounds=(0, 1))
    celoe = milp(c, integrality=np.ones(peremennyh),
                 bounds=Bounds(0, 1),
                 constraints=[LinearConstraint(kazhdyj, 1, 1),
                              LinearConstraint(svyaz, -np.inf, 0)])
    sklady = [round(float(v), 2) for v in lp.x[:m]]
    print(imya, round(lp.fun, 2), round(celoe.fun, 2), sklady)

слабая 49.0 61.0 [0.2, 0.2, 0.6]
сильная 61.0 61.0 [0.0, 0.0, 1.0]

Целый оптимум у обеих записей один - 61: открыть третий склад (35) и возить всё с него (9 + 8 + 3 + 2 + 4 = 26). Модели равносильны.

Релаксации - нет. Слабая открыла все три склада «чуть-чуть»: по 0,2, 0,2 и 0,6. При y = 0,2 склад вправе обслужить целого клиента, и каждый клиент получил доставку с самого дешёвого для него склада, заплатив за открытие пятую часть. Граница 49 - на 20% ниже оптимума. Сильная запись не дала так схитрить: доставка от склада ограничена его долей открытия для каждого клиента отдельно, и граница совпала с оптимумом. Решателю, получившему сильную запись, ветвиться вообще не пришлось. На трёх складах это неважно, на трёхстах - разница между секундами и сутками.

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

«Большое M» взято с запасом. x ≤ 1000000 · y вместо x ≤ 10 · y верно для целых и бесполезно для дробных: y = 0,00001 разрешает всё. Релаксация слабеет, а решатель начинает путаться в округлениях. M берут наименьшим, при котором ограничение ещё верно.

Записали произведение переменных. x_A · x_B - уже не линейная программа. Нужна новая переменная и три неравенства из таблицы.

Забыли целочисленность. Без неё решатель честно вернёт x_B = 0,5. Это не ошибка решателя, это ответ на другую задачу.

Ждут от решателя чуда на слабой модели. Решатель тот же, а время отличается на порядки: от записи зависит, сколько ветвей он отсечёт.

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

Расписания смен, размещение складов, раскрой материала, выбор проектов в портфель, распределение рекламного бюджета - стандартные задачи, которые в промышленности решают целочисленными программами. Решатели: HiGHS (внутри scipy), SCIP, CBC, OR-Tools - открытые; Gurobi и CPLEX - коммерческие. Задачи с тысячами переменных они часто решают за секунды. Твоя работа в такой задаче - не алгоритм, а модель: правильные переменные, сильная запись связей, честные границы. И ещё одна польза: решатель отчитывается разрывом между найденным ответом и границей - «оптимум не больше чем на 0,3% лучше».

Теперь сам

1. Запиши неравенством: «если взят проект A, то не взят проект B».
2. Запиши: «из пяти кандидатов взять ровно двух».
3. Склад обслуживает до 4 клиентов, клиентов 10. Как выглядит сильная запись связи «только открытый» вместе с ограничением вместимости?