Тренды

Сравниваем целочисленное и линейное программирование в Python

26.06.2022

Почему линейное программирование называется так? Оба слова могут вводить в заблуждение.

Линейное подразумевает существование нелинейного.
В данном контексте программирование фактически означает “планирование”.

Если говорить вкратце, линейность не имеет никакого отношения к работе с кодом. Все дело в оптимизации переменных с различными ограничениями.

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

Запустить код из этого руководства можно с помощью этого блокнота Google Colab.

I. Типы задач по оптимизации

В статье “Введение в линейное программирование в Python” мы оптимизировали состав армии. Получился такой результат:

================= Solution =================
Solved in 87.00 milliseconds in 2 iterations

Optimal power = 1800.0 💪power
Army:
 - 🗡️Swordsmen = 6.0000000000000036
 - 🏹Bowmen = 0.0
 - 🐎Horsemen = 5.999999999999999

Почему получилось 5,999 всадников? С помощью VarInt мы определили, что переменные должны быть целочисленные. Что же пошло не так?

Проблема кроется не в модели, а в выборе решателя.

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

Это и есть разница между линейным программированием (linear programming, LP) и целочисленным линейным программированием (integer linear programming, ILP). Короче говоря, решатели LP в качестве переменных могут использовать только вещественные, а не целые числа. Так почему мы определили переменные целыми числами, если они не берутся в расчет?

GLOP не решает задачи ILP в отличие от других решателей. Многие из них — решатели частично-целочисленного линейного программирования (mixed-integer linear programming, MILP/MIP). Это значит, что они учитывают как непрерывные (действительные числа), так и дискретные (целочисленные) переменные. Частным случаем дискретных значений являются булевы переменные, представляющие решения со значениями 0–1.

Другие решатели, такие как SCIP и CBC, могут решать задачи как MIP, так и MINLP (mixed integer nonlinear programming, частично-целочисленное нелинейное программирование). Благодаря OR-Tools можно использовать ранее созданную модель и просто поменять решатель на SCIP или CBC.

# Создаем линейный решатель, используя бэкенд CBC
solver = pywraplp.Solver('Maximize army power', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

# 1. Создаем переменные для оптимизации
swordsmen = solver.IntVar(0, solver.infinity(), 'swordsmen')
bowmen = solver.IntVar(0, solver.infinity(), 'bowmen')
horsemen = solver.IntVar(0, solver.infinity(), 'horsemen')

# 2. Добавляем ограничения для каждого ресурса 
solver.Add(swordsmen*60 + bowmen*80 + horsemen*140 <= 1200)
solver.Add(swordsmen*20 + bowmen*10 <= 800)
solver.Add(bowmen*40 + horsemen*100 <= 600)

# 3. Максимизируем целевую функцию 
solver.Maximize(swordsmen*70 + bowmen*95 + horsemen*230)

# Решаем задачу
status = solver.Solve()

# Если оптимальное решение найдено, выводится результат 
if status == pywraplp.Solver.OPTIMAL:
  print('================= Solution =================')
  print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
  print()
  print(f'Optimal value = {solver.Objective().Value()} 💪power')
  print('Army:')
  print(f' - 🗡️Swordsmen = {swordsmen.solution_value()}')
  print(f' - 🏹Bowmen = {bowmen.solution_value()}')
  print(f' - 🐎Horsemen = {horsemen.solution_value()}')
else:
  print('The solver could not find an optimal solution.')

================= Solution =================
Solved in 3.00 milliseconds in 0 iterations

Optimal value = 1800.0 💪power
Army: 
 — 🗡️Swordsmen = 6.0
 — 🏹Bowmen = 0.0
 — 🐎Horsemen = 6.0

Переменные — все еще плавающие величины (type(swordsmen.solution_value()) = float). Однако у них больше нет странных десятичных дробей: CBC определил их как целые числа.

Таким образом, мы просто округлили эти значения, поскольку ошибка незначительна. Тем не менее нужно подходить с умом к выбору решателя:

LP — для непрерывных переменных;
MIP/MILP — когда есть и непрерывные, и дискретные переменные.

Существуют и другие типы задач: квадратичная (QP) или нелинейная (NLP/MINLP с экспоненциальной целевой функцией или ограничениями). Они применяются при различных обстоятельствах, но следуют таким же принципам работы, как решатели LP и MIP.

II. Построение общей модели

Что произойдет, если ресурсы или цена юнитов изменятся? Или если всадники буду модернизированы, а их сила увеличится?

Одним из лучших преимуществ OR-Tools является использование им языка программирования общего назначения Python. Вместо статических чисел хранить параметры можно в объектах типа словарей и списков.

Код станет не таким читаемым, но намного более гибким, что позволяет решать целый класс задач по оптимизации, не меняя модели (только параметры).

Преобразуем входные параметры в списки Python и через функцию передадим их решателю:

# Входные параметры
UNITS = ['🗡️Swordsmen', '🏹Bowmen', '🐎Horsemen']

DATA = [[60, 20, 0, 70],
        [80, 10, 40, 95],
        [140, 0, 100, 230]]

RESOURCES = [1200, 800, 600]


def solve_army(UNITS, DATA, RESOURCES):
  # Создаем линейный решатель, используя бэкенд CBC
  solver = pywraplp.Solver('Maximize army power', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

  # 1. Создаем переменные для оптимизации
  units = [solver.IntVar(0, solver.infinity(), unit) for unit in UNITS]

  # 2. Добавляем ограничения для каждого ресурса
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum(DATA[u][r] * units[u] for u, _ in enumerate(units)) <= RESOURCES[r])

  # 3. Максимизируем целевую функцию
  solver.Maximize(sum(DATA[u][-1] * units[u] for u, _ in enumerate(units)))

  # Решаем задачу
  status = solver.Solve()

  # Если оптимальное решение найдено, выводится результат
  if status == pywraplp.Solver.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
    print()
    print(f'Optimal value = {solver.Objective().Value()} 💪power')
    print('Army:')
    for u, _ in enumerate(units):
      print(f' - {units[u].name()} = {units[u].solution_value()}')
  else:
    print('The solver could not find an optimal solution.')

solve_army(UNITS, DATA, RESOURCES)

================= Solution =================
Solved in 2.00 milliseconds in 0 iterations

Optimal value = 1800.0 💪power 
Army:
 — 🗡️Swordsmen = 6.0
 — 🏹Bowmen = 0.0
 — 🐎Horsemen = 6.0

Мы получили такой же результат: код, похоже, работает. Теперь изменим параметры, чтобы немного усложнить задачу.

Представьте, что ресурсов намного больше: 183000 еды, 90512 древесины и 80150 золота, и теперь можно произвести множество дополнительных юнитов! Вот новая таблица:

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

Для примера возьмем 10: сила = 10 * атака + здоровье. Целевая функция выглядит так:

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

UNITS = [
    '🗡️Swordsmen',
    '🛡️Men-at-arms',
    '🏹Bowmen',
    '❌Crossbowmen',
    '🔫Handcannoneers',
    '🐎Horsemen',
    '♞Knights',
    '🐏Battering rams',
    '🎯Springalds',
    '🪨Mangonels',
]

DATA = [
    [60, 20, 0, 6, 70],
    [100, 0, 20, 12, 155],
    [30, 50, 0, 5, 70],
    [80, 0, 40, 12, 80],
    [120, 0, 120, 35, 150],
    [100, 20, 0, 9, 125],
    [140, 0, 100, 24, 230],
    [0, 300, 0, 200, 700],
    [0, 250, 250, 30, 200],
    [0, 400, 200, 12*3, 240]
]

RESOURCES = [183000, 90512, 80150]


def solve_army(UNITS, DATA, RESOURCES):
  # Создаем линейный решатель, используя бэкенд CBC
  solver = pywraplp.Solver('Maximize army power', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

  # 1. Создаем переменные для оптимизации
  units = [solver.IntVar(0, solver.infinity(), unit) for unit in UNITS]

  # 2. Добавляем ограничения для каждого ресурса
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum(DATA[u][r] * units[u] for u, _ in enumerate(units)) <= RESOURCES[r])

  # 3. Максимизируем новую целевую функцию
  solver.Maximize(sum((10*DATA[u][-2] + DATA[u][-1]) * units[u] for u, _ in enumerate(units)))

  # Решаем задачу
  status = solver.Solve()

  # Если оптимальное решение найдено, выводится результат
  if status == pywraplp.Solver.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
    print()
    print(f'Optimal value = {solver.Objective().Value()} 💪power')
    print('Army:')
    for u, _ in enumerate(units):
      print(f' - {units[u].name()} = {units[u].solution_value()}')
  else:
    print('The solver could not find an optimal solution.')

solve_army(UNITS, DATA, RESOURCES)

================= Solution =================
Solved in 74.00 milliseconds in 412 iterations

Optimal value = 1393145.0 💪power
Army:
 — 🗡️Swordsmen = 2.0
 — 🛡️Men-at-arms = 1283.0
 — 🏹Bowmen = 3.0
 — ❌Crossbowmen = 0.0
 — 🔫Handcannoneers = 454.0
 — 🐎Horsemen = 0.0
 — ♞Knights = 0.0
 — 🐏Battering rams = 301.0
 — 🎯Springalds = 0.0
 — 🪨Mangonels = 0.0

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

III. Комбинирование ограничений

Допустим, мы выяснили, что сила армии противника составляет 1000000. Можно создать более мощную армию, но ресурсы ценны, и это было бы не эффективно. Достаточно лишь собрать армию мощнее, чем 1000000 (даже 1000001 будет достаточно).

Другими словами, общая мощность теперь является ограничением (💪 > 1000000), а не объектом максимизации. Новая цель — минимизировать количество ресурсов, необходимых для создания армии. Однако можно использовать те же входные данные, поскольку они не поменялись.

Новое ограничение можно перевести как “сумма силы выбранных юнитов должна быть строго больше 1000000”:

Чтобы спроектировать это ограничение в коде, можно добавить цикл для юнитов и ресурсов.

Целевая функция тоже поменялась. Наша цель — минимизировать общее количество ресурсов, потраченных на создание армии.

И снова мы создадим цикл для ресурсов, чтобы внедрить их в OR-Tools:

def solve_army(UNITS, DATA, RESOURCES):
  # Создаем линейный решатель, используя бэкенд CBC
  solver = pywraplp.Solver('Minimize resource consumption', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

  # 1. Создаем переменные для оптимизации
  units = [solver.IntVar(0, solver.infinity(), unit) for unit in UNITS]

  # 2. Добавляем ограничения для каждого ресурса
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum((10 * DATA[u][-2] + DATA[u][-1]) * units[u] for u, _ in enumerate(units)) >= 1000001)

  # 3. Минимизируем целевую функцию
  solver.Minimize(sum((DATA[u][0] + DATA[u][1] + DATA[u][2]) * units[u] for u, _ in enumerate(units)))

  # Решаем задачу
  status = solver.Solve()

  # Если оптимальное решение найдено, выводится результат
  if status == pywraplp.Solver.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
    print()

    power = sum((10 * DATA[u][-2] + DATA[u][-1]) * units[u].solution_value() for u, _ in enumerate(units))
    print(f'Optimal value = {solver.Objective().Value()} 🌾🪵🪙resources')
    print(f'Power = 💪{power}')
    print('Army:')
    for u, _ in enumerate(units):
      print(f' - {units[u].name()} = {units[u].solution_value()}')
    print()

    food = sum((DATA[u][0]) * units[u].solution_value() for u, _ in enumerate(units))
    wood = sum((DATA[u][1]) * units[u].solution_value() for u, _ in enumerate(units))
    gold = sum((DATA[u][2]) * units[u].solution_value() for u, _ in enumerate(units))
    print('Resources:')
    print(f' - 🌾Food = {food}')
    print(f' - 🪵Wood = {wood}')
    print(f' - 🪙Gold = {gold}')
  else:
      print('The solver could not find an optimal solution.')

solve_army(UNITS, DATA, RESOURCES)

================= Solution =================
Solved in 4.00 milliseconds in 0 iterations

Optimal value = 111300.0 🌾🪵🪙resources
Power = 💪1001700.0 
Army:
 — 🗡️Swordsmen = 0.0
 — 🛡️Men-at-arms = 0.0
 — 🏹Bowmen = 0.0
 — ❌Crossbowmen = 0.0
 — 🔫Handcannoneers = 0.0
 — 🐎Horsemen = 0.0
 — ♞Knights = 0.0
 — 🐏Battering rams = 371.0
 — 🎯Springalds = 0.0
 — 🪨Mangonels = 0.0

Resources:
 — 🌾Food = 0.0
 — 🪵Wood = 111300.0
 — 🪙Gold = 0.0

Решатель нашел оптимальное решение: нужно построить 371 стенобитное орудие общей стоимостью 111300 древесины. Однако нам не хватает ресурсов! Нам доступно лишь 90512.

Возможно ли учитывать ограниченные ресурсы при создании мощнейшей армии? Это легко: нужно лишь скопировать/вставить ограничения из предыдущей секции.

В этой версии модели представлено два вида ограничений.

Общая сила должна быть больше 1000000.
Мы не можем тратить ресурсов больше, чем имеем.

def solve_army(UNITS, DATA, RESOURCES):
  # Создаем линейный решатель, используя бэкенд CBC
  solver = pywraplp.Solver('Minimize resource consumption', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

  # 1. Создаем переменные для оптимизации
  units = [solver.IntVar(0, solver.infinity(), unit) for unit in UNITS]

  # 2. Добавляем ограничения для каждого ресурса
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum((10 * DATA[u][-2] + DATA[u][-1]) * units[u] for u, _ in enumerate(units)) >= 1000001)

  # Старые ограничения для ограниченных ресурсов
  for r, _ in enumerate(RESOURCES):
    solver.Add(sum(DATA[u][r] * units[u] for u, _ in enumerate(units)) <= RESOURCES[r])

  # 3. Минимизируем целевую функцию
  solver.Minimize(sum((DATA[u][0] + DATA[u][1] + DATA[u][2]) * units[u] for u, _ in enumerate(units)))

  # Решаем задачу
  status = solver.Solve()

  # Если оптимальное решение найдено, выводится результат
  if status == pywraplp.Solver.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.wall_time():.2f} milliseconds in {solver.iterations()} iterations')
    print()

    power = sum((10 * DATA[u][-2] + DATA[u][-1]) * units[u].solution_value() for u, _ in enumerate(units))
    print(f'Optimal value = {solver.Objective().Value()} 🌾🪵🪙resources')
    print(f'Power = 💪{power}')
    print('Army:')
    for u, _ in enumerate(units):
      print(f' - {units[u].name()} = {units[u].solution_value()}')
    print()
    
    food = sum((DATA[u][0]) * units[u].solution_value() for u, _ in enumerate(units))
    wood = sum((DATA[u][1]) * units[u].solution_value() for u, _ in enumerate(units))
    gold = sum((DATA[u][2]) * units[u].solution_value() for u, _ in enumerate(units))
    print('Resources:')
    print(f' - 🌾Food = {food}')
    print(f' - 🪵Wood = {wood}')
    print(f' - 🪙Gold = {gold}')
  else:
      print('The solver could not find an optimal solution.')

solve_army(UNITS, DATA, RESOURCES)

================= Solution =================
Solved in 28.00 milliseconds in 1 iterations

Optimal value = 172100.0 🌾🪵🪙resources
Power = 💪1000105.0
Army:
 — 🗡️Swordsmen = 1.0
 — 🛡️Men-at-arms = 681.0
 — 🏹Bowmen = 0.0
 — ❌Crossbowmen = 0.0
 — 🔫Handcannoneers = 0.0
 — 🐎Horsemen = 0.0
 — ♞Knights = 0.0
 — 🐏Battering rams = 301.0
 — 🎯Springalds = 0.0
 — 🪨Mangonels = 0.0

Resources:
 — 🌾Food = 68160.0
 — 🪵Wood = 90320.0
 — 🪙Gold = 13620.0

Поскольку запасы древесины ограничены, то количество осадных орудий снижается с 371 до 301. Взамен мы получаем 681 солдата и одного мечника.

Общая стоимость армии 172100, что намного больше предыдущей стоимости 111300 (увеличение на 65%), но это оптимальное решение при данных ограничениях. Из этого результата следует, что необходимо производить больше древесины, потому что эти осадные орудия очень рентабельны!

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

IV. Линейное программирование в сравнении с машинным обучением

Можно ли использовать машинное обучение (в широком смысле) вместо линейного программирования? К примеру, эту задачу можно решить с помощью генетического алгоритма.

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

Линейное программирование может находить оптимальное решение за неопределенное количество времени (это может занять годы), а машинное обучение — аппроксимировать сложные функции в мгновение ока.
В LP нет обучения, но требуются знания эксперта, чтобы построить математическую модель. Машинному обучению необходимы данные, но для решения задачи модели могут использоваться как “черные коробки”.
Как правило, задачи, которые не имеют конкретного ограничения по времени и/или чрезвычайно сложные, преимущественно решаются с помощью линейного программирования.

Заключение

В этом руководстве мы разобрались в математической оптимизации.

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

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

Сравниваем целочисленное и линейное программирование в Python

I. Типы задач по оптимизации

II. Построение общей модели

III. Комбинирование ограничений

IV. Линейное программирование в сравнении с машинным обучением

Заключение

Читайте также

8 причин использовать Pydantic для улучшения парсинга и валидации данных

Как создать HTML-таблицу и PDF с Node и Google Puppeteer

Алгоритм YOLO простым языком

Тесты

Насколько хорошо вы разбираетесь в AngularJS?

Насколько хорошо вы разбираетесь в Windows?

Насколько хорошо вы умеете пользоваться Linux?