news.shamcode.ru | Обзор пакетов SciPy, Pyomo и CVXPY для решения задач условной оптимизации

Привет, Habr! На связи Михаил Будылин и Антон Денисов, мы работаем в отделе аналитики данных X5 Tech.

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

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

Обзор пакетов Python

Для начала напомним основу — постановку задачи оптимизации:

Пусть - вектор размерности , $x \in X$ - допустимое множество значений этих переменных.
$f(x) \to \min(\max), f(\cdot)$ - целевая функция.

$g_i(x) \leqslant 0, \ i=1..m$ - ограничения вида неравенств.

$h_i(x) = 0, \ j=1..k$ - ограничения вида равенств.

А так же, классификацию оптимизационных задач, предложенную в предыдущей статье:

Безусловная оптимизация — отсутствуют, $X = \mathbb{R}^n$ ;
LP (linear programming) — линейное программирование. - линейные функции, $X = \mathbb{R}_+^n$ ;
MILP (mixed integer linear programming) — смешанное целочисленное линейное программирование, это задача LP в которой часть переменных являются целочисленными;
NLP (nonlinear programming) — нелинейное программирование, возникает когда хотя бы одна из функций $f(x),\ g_i(x),\ h_j(x)$ нелинейна;
MINLP (mixed integer nonlinear programming) — смешанное целочисленное нелинейное программирование, возникает как и в случае с MILP, когда часть переменных принимает целочисленные значения;

Теперь добавим в полученную картину солверы и отобразим в виде схемы:

Последовательность разбираемых в статье солверов

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

SciPy. Обзор.

SciPy — одна из первых библиотек в Python, знакомство с которой начинается у специалистов в области Data Science: она содержит большой набор функций для научных вычислений, в том числе имеет инструменты для решения оптимизационных задач, находящихся в модуле scipy.optimize. В этом модуле имеются методы для решения задач как нелинейного программирования (NLP), так и линейного программирования (LP), в том числе задач смешанного целочисленного линейного программирования (MILP).

Среди солверов, которые поддерживают решение задач условной оптимизации (NLP) можно выделить COBYLA, SLSQP, trust-constr. Все три метода реализованы в пакете
scipy.optimize.minimize и определяются параметром method ('cobyla','slsqp' и 'trust-constr', соответственно).

COBYLA — это метод, позволяющий производить оптимизацию функции, градиент которой неизвестен, т.е., по сути, заниматься оптимизацией «черного ящика».
Для практического применения важно, что этот метод:

не поддерживает ограничения типа равенства: . Связано это с особенностями реализации солвера, который внутри себя вызывает SciPy. Эту проблему можно обойти, заменив ограничение на два ограничения вида неравенства: $a(x) \leq b$ и $a(x) \geq b$ , которое равносильно равенству.
не поддерживает границы для переменных задаваемых через параметр bounds в функции minimize (из модуля scipy.optimize) — их необходимо задавать через линейные ограничения, например с помощью LinearConstraint.

SLSQP и trust-constr — это методы уже второго порядка, соответственно, для них требуется существование и непрерывность 1-ой и 2-ой производных.

Подробный обзор методов, литературы и примеры применения солверов можно найти тут.

SciPy. Решение задачи NLP.

Рассмотрим пример использование методов scipy для решения задачи ценообразования из статьи 1. Задача заключалась в максимизации выручки с сохранением текущего уровня маржи с ограниченным диапазоном изменения цены в пределах ±10% от текущей.

Постановка была следующей:

— количество товаров,

$P_{0, i}$ — текущая цена товара ,

$Q_{0, i}$ — текущие продажи товара ,

$x_{i} = P_{i} / P_{0, i}$ — отношение новой цены $P_{i}$ к текущей $P_{0, i}$ ,

$E_{i}$ — параметр для пересчёта спроса на товар по формуле $Q_{i}(x) = Q_{0, i} \cdot \exp(E_{i} \cdot (x_{i} - 1))$ ,

$R_{i}(x_{i}) = \sum_{i=1}^{n} P_{0, i} \cdot x_{i} \cdot Q_{i}(x_i)$ — выручка, которая для каждого товара определяется как произведение его цены на продажи,

$M_{i}(x_{i}) = \sum_{i=1}^{n} (P_{0, i} \cdot x_{i} - C_{i})\cdot Q_{0, i} \cdot Q_{i}(x_i)$ — маржа, которая считается как разность между ценой продажи и себестоимостью, умноженной на продажи,

$M_{0} = \sum_{i=1}^{n} M_{i}(x_{i}=1)$ — текущая маржа, т.е. её значение при текущей цене,

$\begin{cases} \tag{1} \sum_{i=1}^{n} R_{i}(x_{i}) \to \max, \\ \sum_{i=1}^{n} M_{i}(x_{i}) \geqslant M_{0}, \\ x_i \in [0.9, 1.1], \ i=1..n\\ \end{cases}$

Оформим все функции для решения задачи в виде кода и зададим начальные условия:

# мини пример из постановки первой статьи
import scipy.optimize as scopt
import numpy as np

# Количество товаров
N = 3
# задаем параметры E, используемые в формуле Q = Q0 * exp(E * (x - 1))
E = np.array([-3., -1., -0.5])
# текущие цены
P0 = np.array([100., 100., 100.])
# текущие продажи
Q0 = np.array([50000., 200000., 30000.0])
# себестоимость
C = np.array([90.0, 80.0, 70.0])
# текущая выручка
R0 = np.sum(P0 * Q0)
# текущая маржа
M0 = np.sum((P0 - C) * Q0)

# выручка - целевая функция, задаем возможность "управлять" масштабом через 'scale'
def f_obj(x, args):
    f = - args['scale'] * np.sum(Q0 * P0 * x * np.exp(E * (x - 1.)))
    return f
obj = f_obj

# функция для ограничения по марже, по умолчанию нормируем ограничения на текущую выручку
def f_cons(x):
    f = np.sum(Q0 * (P0 * x - C) * np.exp(E * (x - 1.0))) / R0
    return f

con_mrg = scopt.NonlinearConstraint(f_cons, lb=M0 / R0, ub=np.inf)
# поиск новой цены производим в диапазоне 90% - 110% от текущей цены
x_bounds = [(0.9, 1.1)] * 3
# ограничение для переменных в cobyla
con_bnd = scopt.LinearConstraint(np.eye(3), lb=[0.9] * 3, ub=[1.1] * 3)
# стартовая точка для поиска
x0 = [1.0] * 3

def print_results(model, solver_name, M0, R0):
    print(f"Решение {solver_name}: {list(np.round(model['x'], 4))}; model message: {model['message']}")
    print(f"Выручка: {round(-f_obj(model['x'], {'scale': 1}), 0)}")
    print(f"Ограничение на маржу M >= {M0}; значение маржи: {round(R0 * con_mrg.fun(model['x']), 0)}")

Применение методов SciPy:

# применение солвера cobyla
res_cobyla = scopt.minimize(obj, x0, constraints=[con_bnd, con_mrg], method='cobyla', args={'scale': 1})
print_results(res_cobyla, 'cobyla', M0, R0)

# output:
# Решение cobyla: [0.9, 1.0165, 1.1]; model message: Optimization terminated successfully.
# Выручка: 29210742.0
# Ограничение на маржу M >= 5400000.0; значение маржи: 5399998.0

# применение солвера trust-constr
res_trust_const = scopt.minimize(obj, x0, bounds=x_bounds, constraints=[con_mrg], method='trust-constr', args={'scale': 1})
print_results(res_trust_const, 'trust-constr', M0, R0)

# output:
# Решение trust-constr: [0.9, 1.0165, 1.1]; model message: `gtol` termination condition is satisfied.
# Выручка: 29210742.0
# Ограничение на маржу M >= 5400000.0; значение маржи: 5400005.0

# применение солвера slsqp
res_slsqp = scopt.minimize(obj, x0, bounds=x_bounds, constraints=[con_mrg], method='slsqp', args={'scale': 1})
print_results(res_slsqp, 'slsqp', M0, R0)

# output
# Решение slsqp: [1.0, 1.0, 1.0]; model message: Optimization terminated successfully
# Выручка: 28000000.0
# Ограничение на маржу M >= 5400000.0; значение маржи: 5400000.0

По результатам оптимизации можно выделить 2 основных момента:

методы COBYLA и trust-constr сходятся к одному и тому же решению. Также видно, что выполняются все ограничения (по границам и по марже - с точностью более 7-8 знаков),
SLSQP в данном случае «остался» в начальной точке — это связано с критерием на остановку при достижении определённого уровня изменения целевой функции и градиентов, которые регулируются через параметры ftol и eps в аргументе options.

Вообще говоря, в большинстве реализаций солверов производится масштабирование функций к значению порядка 1, что позволяет повысить устойчивость и сходимость методов. Попробуем применить масштабирование в задаче путём деления целевой функции (текущей выручки) на выручку при текущих ценах $R_{0}$ .

# применение солвера slsqp
res_slsqp = scopt.minimize(obj, x0, bounds=x_bounds, constraints=[con_mrg], method='slsqp', args={'scale': 1. / R0})
print_results(res_slsqp, 'slsqp', M0, R0)

# output
# Решение slsqp: [0.9, 1.0165, 1.1]; model message: Optimization terminated successfully
# Выручка: 29210742.0
# Ограничение на маржу M >= 5400000.0; значение маржи: 5399999.0

SciPy. Решение задачи MILP.

Для решения задач линейного программирования в подмодуле optimize имеется функция linprog, начиная с версии SciPy==1.9.0 появилась возможность решения задачи линейного целочисленного программирования с помощью функции milp и linprog.
В качестве солвера для LP(MILP) по умолчанию используется HiGHS — в нём реализованы симплекс метод (highs-ds) и метод внутренней точки (highs-ipm). При запуске решения по умолчанию выбирается один из методов, возможность задавать метод принудительно также присутствует.

Для примера рассмотрим небольшую задачку смешанного целочисленного рюкзака, где для наглядности применения метода linprog используем переменные разных типов:

$\begin{cases} \tag{2} x = (x_{1}, x_{2}, x_{3}, x_{4}) \\ \\ F(x) = 1 \cdot x_{1} + 2 \cdot x_{2} + 3 \cdot x_{3} + 1 \cdot x_{4} \to \max \\ \\ 3 \cdot x_1 + 1 \cdot x_2 + 2 \cdot x_3 + 2 \cdot x_4 \leqslant 7.0 \\ \\ x_{1} \in \{0, 1, 2\}, x_{2}, x_{3} \in \{0, 1\}, x_{4} \in [0.0; 0.5] \end{cases}$

# пример с "рюкзаком", у которого типы переменных различаются
# так как scipy решет задачу минимизации, не забываем ставить минус для коэффициентов в целевой функции для максимизации
c = -np.array([1., 2., 3., 1.])
A_ub = np.array([[3., 1., 2., 2.]])
b_ub = np.array([7.0])
# проставляем индикаторы для типа переменных, 0 - непрерывное, 1 - целое число
var_types = [1, 1, 1, 0]
# также указываем границы, в том числе и для целочисленных переменных
bounds = [(0, 2), (0, 1), (0, 1), (0, 0.5)]
res_milp = scopt.linprog(c=c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, integrality=var_types, method='highs')
print(f"Решение: x = {list(np.round(res_milp['x'], 2))}, f = {-res_milp['fun']}, {res_milp['message']}")

# output
# Решение: x = [1.0, 1.0, 1.0, 0.5], f = 6.5, Optimization terminated successfully. (HiGHS Status 7: Optimal)

Как мы видим, задача успешно решена. Попробуем убедиться в правильности, подобрав решение перебором, так как вариантов не так много: для первой переменной 3, второй и третьей — 2, итого 3 * 2 * 2 = 12 вариантов. И, поскольку, решается задача максимизации с положительными коэффициентами в целевой функции, то непрерывную переменную нужно брать «максимальной», пока не упрёмся либо в ограничение, либо в допустимые границы переменной. В результате перебора необходимо будет отсечь решения, которые нарушают ограничение на сумму, они будут отмечены меткой is_valid=False:

Код реализующий перебор

# реализуем полный перебор
import dataframe_image as dfi
import itertools
import pandas as pd

# генерируем всевозможные варианты для целочисленных переменных
sols = pd.DataFrame(list(itertools.product([0., 1., 2.], [0., 1.], [0., 1.])), columns=['x1', 'x2', 'x3'])
# часть суммы из ограничения, связанная с целочисленными переменными
sols['a'] = sols[['x1', 'x2', 'x3']].multiply(np.array([3., 1., 2.]), axis=1).sum(axis=1)
# непрерывная переменная формируется по остаточному принципу по ограничению для максимизации целевой функции
# здесь же учитываем, что максимальное значение равно 0.5
sols['x4'] = np.clip((7.0 - sols['a']) / 2.0, 0.0, 0.5)
# подсчет целевой функции и поиск оптимального решения с выводом всех вариантов
# где через is_valid обозначены допустимые решения
sols['f'] = sols[['x1', 'x2', 'x3', 'x4']].multiply(np.array([1., 2., 3., 1.]), axis=1).sum(axis=1)
sols['a'] = sols['a'] + 2 * sols['x4']
sols['is_valid'] = sols['a'] <= 7.0
i_max = sols[sols['is_valid']]['f'].idxmax()
sols = sols[['x1', 'x2', 'x3', 'x4', 'a', 'is_valid', 'f']].transpose()
res = sols.style.set_properties(**{'background-color': '#5eae76'}, subset=[i_max]).format(precision=2)
dfi.export(res, './images/milp_brute_force.png')

Здесь по вертикали значения переменных, a -значения в левой части ограничения, и последней строчкой значение функции. По горизонтали варианты перебора

Как мы видим оптимальным решением задачи является, то решение которые мы получили с помощью солвера.

Pyomo. Обзор.

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

Pyomo входит в проект COIN-OR, содержащий ряд солверов, среди которых выделим два:

Ipopt - находит локальные оптимумы в задаче NLP с помощью прямо-двойственного метода внутренней точки (подробнее в оригинальной статье).

Cbc - решает задачи MILP, на базе алгоритма, сочетающем в себе метод ветвей и границ и секущих плоскостей (подробнее на wiki).

Также для решения задач LP/MILP имеется поддержка пакета GLPK.
Отметим ещё один солвер bonmin, который построен на базе солверов Cbc и Ipopt, что позволяет браться за задачи MINLP.

Процесс построения оптимизационной модели в Pyomo состоит из следующих основных этапов:

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

Рассмотрим шаги на примере задач (1) и (2).

Pyomo. Решение задачи NLP.

# пример для задачи (1)
import pyomo.environ as pyo

# объявление объекта - модели 
model = pyo.ConcreteModel('model')

# задаем переменные, в данном случае они все непрерывные, инициализируем 1.0
model.x = pyo.Var(range(N), domain=pyo.Reals, bounds=x_bounds, initialize=1)

# объявление целевой функции и передача в модель
obj_expr = sum(P0[i] * model.x[i] * Q0[i] * pyo.exp(E[i] * (model.x[i] - 1)) for i in model.x)
model.obj = pyo.Objective(expr=obj_expr, sense=pyo.maximize)

# объявление ограничения и передача в модель
con_expr = sum((P0[i] * model.x[i] - C[i]) * Q0[i] * pyo.exp(E[i] * (model.x[i] - 1)) for i in model.x) >= M0
model.con = pyo.Constraint(expr=con_expr)

# запуск солвера ipopt для решения поставленной оптимизационной задачи
solver = pyo.SolverFactory('ipopt')
res = solver.solve(model)

# получение ответа - результата решения задачи
x_opt = [round(model.x[i].value, 3) for i in model.x]

print(f"Решение ipopt: {list(x_opt)} {res.solver.termination_condition}")
print(f"Выручка: {round(model.obj(), 0)}")
print(f"Ограничение на маржу M >= {M0}; значение маржи: {round(model.con(), 0)}")

# output
# Решение ipopt: [0.9, 1.016, 1.1] optimal
# Выручка: 29210742.0
# Ограничение на маржу M >= 5400000.0; значение маржи: 5400000.0

Pyomo. Решение задачи NLP.

'''
MILP пример с рюкзаком
пример с "рюкзаком", у которого типы переменных различаются
так как в pyomo не реализована возможность указывать в одном векторе значения разных типов
то необходимо их описывать отдельно и данные, соответственно, для удобства тоже
'''
c_i, c_c = np.array([1., 2., 3.]), np.array([1.])
A_i, A_c = np.array([3., 1., 2.]), np.array([2.])
b = 7.0
# объявление объекта - модели 
model = pyo.ConcreteModel('model')

# формирование переменных - отдельно целочисленные и непрерывные
bounds_i = [(0.0, 2.0), (0.0, 1.0), (0.0, 1.0)]
bounds_c = [(0.0, 0.5)]
model.x_i = pyo.Var(range(3), domain=pyo.Integers, bounds=bounds_i)
model.x_c = pyo.Var(range(1), domain=pyo.Reals, bounds=bounds_c)

# объявление целевой функции и передача в модель для максимизации
obj_expr = sum(c_i[i] * model.x_i[i] for i in model.x_i) +\
           sum(c_c[i] * model.x_c[i] for i in model.x_c)
model.obj = pyo.Objective(expr=obj_expr, sense=pyo.maximize)

# объявление
con_expr = sum(A_i[i] * model.x_i[i] for i in model.x_i) +\
           sum(A_c[i] * model.x_c[i] for i in model.x_c) <= b
model.con = pyo.Constraint(expr=con_expr)

solver = pyo.SolverFactory('glpk')
res = solver.solve(model)
x_opt = [model.x_i[i].value for i in model.x_i] + [model.x_c[i].value for i in model.x_c]

print(f"Решение: x = {x_opt}, f = {model.obj()}, {res.solver.termination_condition}")

# output
# Решение: x = [1.0, 1.0, 1.0, 0.5], f = 6.5, optimal

Видим, что солверы справились с обозначенными задачами и выдали результаты, аналогичные пакетам SciPy.

Pyomo. Дополнительные особенности.

Из дополнительного можно отметить, что Pyomo имеет подмодуль GDP (Generalized Disjunctive Programming), который позволяет моделировать логические правила и задавать ограничения, которые должны при этом выполняться. Например, в простейшем случае, данный подход может быть применён, когда необходимо выбрать одно из действий, каждое из которых описывается своей системой ограничений. Разберём, чем может быть полезен данный подход на примере следующей задачи:

$\begin{cases} \tag{3} f_1(x_1) = x_1^2 \\ f_2(x_2) = (x_2 - 1) ^ 2 \\ f_3(x_3) = (x_3 + 2) ^ 2 \\ F(x_1, x_2, x_3) = f_1(x_1) + f_2(x_2) + f_3(x_3) \\ x_i \in [-3; -1.5 ] \cup \{0\} \cup [ 1.5; 3] , \ i = 1..3 \end{cases}$

Как несложно заметить, здесь в определении области переменных есть «разрыв», который тривиальным образом описать невозможно, необходимо вводить условие на то, в какой из трёх областей $[-3; -1.5 ], \{0\}, [ 1.5; 3]$ необходимо производить поиск решения. Так как исходная функция является суммой независимых друг от друга функций, то минимум достигается в минимуме каждой из функций . Отобразим графики с учётом условий на ниже: на них можно видеть, что достигает своего минимума в , в , в .

Процесс формирования оптимизационной модели для задачи (3) концептуально похож на тот, что был для задач (1) и (2).
Опишем его кодом ниже, снабдив необходимыми комментариями:

# решение задачи о рюкзаке с помощью pyomo
import pyomo.environ as pyo
from pyomo.gdp import Disjunction

# формирование объекта модели pyomo
model = pyo.ConcreteModel('gdp_sample')
# общий диапазон для всех переменных [-3.0, 3.0]
model.x = pyo.Var(range(0, 3), domain=pyo.Reals, bounds=(-3., 3.))

# задание целевой функции
a = [0.0, 1.0, -2.0]
obj_expr = sum((model.x[i] - a[i]) ** 2 for i in model.x)
model.obj = pyo.Objective(expr=obj_expr, sense=pyo.minimize)

# формирование правил для диапазонов поиска переменных с помощью функции Disjunction.
# в объекте Disjunction необходимо перечислить условия в виде ограничений
# всего у нас три переменных, поэтому задаем три группы условий, в каждой из которых будет выполняться только одно
model.djn = Disjunction(range(3))

# непосредственно задание условий для каждой из переменных
for i in range(3):
    model.djn[i] = [model.x[i] <= -1.5, model.x[i] == 0, model.x[i] >= 1.5]

'''
Перед началом решения необходимо произвести преобразование модели для дальнейшей передачи в солвер.
В данном случае преобразование произведем т.н. методом big-M, который позволяет
перекодировать логические условия в алгебраические,
сводя задачу к MINLP. Более подробно, как устроен метод и пример формирования условий руками рассмотрим в следующей статье.
Здесь ниже оставим ссылку для ознакомления с подходом.
'''
pyo.TransformationFactory('gdp.bigm').apply_to(model)
solver = pyo.SolverFactory('bonmin')
res = solver.solve(model)
x_opt = [round(model.x[i].value, 3) for i in model.x]
print('Решение:', x_opt, ', f =', round(model.obj(), 2))

# output
# Решение: [-0.0, 1.5, -2.0] , f = 0.25

Как видим, решение для совпадает с ожидаемым. А вот ссылка на метод big-M.

CVXPY. Обзор.

CVXPY - данный пакет был реализован для решения задач выпуклой оптимизации.

Для решения задачи, по аналогии с Pyomo, необходимо выполнить несколько шагов: определить переменные, задать целевую функцию и ограничения для формирования объекта оптимизационной задачи. После того, как задача сформулирована, перед запуском солвера проверяется выпуклость целевой функции и ограничений с помощью правил DCP (Disciplined Convex Programming). По сути, это набор правил, которые однозначно гарантируют выпуклость функции, здесь стоит отметить, что эти правила - необходимые условия, но не достаточные, почему это так рассмотрим в примере ниже. После проверки задача преобразуется для передачи солверу в формат, поддерживающий решение задач выпуклой оптимизации, в том числе CVXPY поддерживает коммерческие солверы. Полный список солверов можно посмотреть на этой странице, там же указаны все типы задач, которые позволяет решать CVXPY. Обратим внимание, что для данного пакета не так важно наличие непрерывности первой и второй производной, важна только выпуклость.

Также по аналогии с предыдущими пакетами рассмотрим применение CVXPY на примерах.

CVXPY. Решение задачи NLP.

Задачу (1) с помощью CVXPY решить не удастся, так как целевая функция и ограничения в общем случае не являются выпуклыми и, соответственно, не подчиняются правилам DCP. Целевая функция состоит из суммы слагаемых вида $f(x) = x e^{E(x-1)}$ , вторая производная которой равна:

$f''(x) = E x e^{E(x-1)}(Ex + 2)$ и меняет свой знак в точке $x_{0} = \frac{-2}{E} = \frac{2}{|E|} \geqslant 0$

(с учётом того, что $E \leqslant 0$ ) — то есть в общем случае не является выпуклой.

Рассмотрим тогда другой пример: необходимо добраться из точки с координатами (0, 0) в точку (1, 2), область разделена на две части вдоль линии y = 1; до y = 1 максимальная скорость равна 1, а после в k (в примере k=1.5) раз меньше.
Очевидно, что самый быстрый путь — это движение по прямой линии с максимально скоростью на каждом из участков, при этом на границе не должно быть «разрывов» траектории. В таком случае оптимизационную задачу можно сформулировать следующим образом:

$\begin{cases} \tag{4} T = \sqrt{(x-x_{1})^2 + (y-y_{1})^2} + k \cdot \sqrt{(x-x_{2})^2 + (y-y_{2})^2} \to min, \\ x_{1} = 0,\ y_{1} = 0,\ x_{2} = 1,\ y_{2} = 2, \\ y = 1 \\ x \in [0; 2] \end{cases}$

Здесь в — первое и второе слагаемые — это путь по прямой от начальной/конечной точки до точки на границе, делённый на скорость, другими словами — время, затрачиваемое прохождение пути на 1-ом и 2-ом участках. Каждое из слагаемых является выпуклой функцией, а значит и их сумма тоже.

# решение задачи с разрывом в области определения с помощью cvxpy
import cvxpy as cp

X1, Y1 = 0.0, 0.0
X2, Y2 = 1.0, 2.0
K = 1.5

# задаем одну переменную x
X = cp.Variable(1)
Y_ = 1.0

# целевая функция - корни из суммы квадратов - являются выпуклыми
objective = cp.sqrt(cp.square(X - X1) + cp.square(Y_ - Y1) ** 2) +\
            cp.sqrt(cp.square(X2 - X) + cp.square(Y2 - Y_) ** 2)

# здесь единственное ограничение - это диапазон для x
constraints = []
constraints.extend([X >= 0.0, X <= X2])

# объявляем оптимизационную задачу, здесь берем минимизацию целевой функции
problem = cp.Problem(cp.Minimize(objective), constraints)

# совершаем проверку задачи на выпуклость согласно правилам DCP
print(f"is dcp: {problem.is_dcp()}")

# output
# is dcp: False

Как мы видим, задача не прошла проверку на выпуклость согласно правилам DCP, так как берётся вогнутая функция(квадратный корень) от выпуклой (квадратичной), но данная проблема решается просто вызовом функции norm для расчёта длины вектора (в нашем контексте расстояние), про которую CVXPY известно, что она является выпуклой.
В модуле есть также ещё дополнительный набор функций, которые не подчиняются DCP, но при этом выпуклые, например, log_sum_exp - логарифм от суммы экспонент. Такие моменты необходимо учитывать при формулировке задач в CVXPY.

# скорректируем целевую функцию через вызов norm
objective = cp.norm(cp.hstack([X - X1, Y_ - Y1]), 2) + K * cp.norm(cp.hstack([X2 - X, Y2 - Y_]), 2)

# формируем ограничения и формируем задачу
constraints = []
constraints.extend([X >= 0.0, X <= X2])
problem = cp.Problem(cp.Minimize(objective), constraints)

# проверка на выпуклость
print(f"is dcp: {problem.is_dcp()}")

# решаем задачу путем вызова солвера ECOS
sol = problem.solve('ECOS')

# извлекаем решение
x_opt = X.value[0]
print(f"Решение: x = {round(x_opt, 3)}, f = {round(sol, 2)}")

# output
# is dcp: True
# Решение: x = 0.623, f = 2.78

CVXPY. Решение задачи MILP

Также разберём, как решается задача (2) с помощью CVXPY, так как пакет позволяет решать задачи LP (MILP):

# объявление переменных, отдельно целочисленных и непрерывных
x_i = cp.Variable(3, integer=True)
x_c = cp.Variable(1, nonneg=True)

# коэффициенты для целевой функции и ограничений
c = np.array([1., 2., 3., 1.])
A = np.array([3., 1., 2., 2.])
b = 7.0

# максимизация функции - сумма по целочисленной и непрерывной части переменных
obj = cp.Maximize(cp.sum(c[:3] @ x_i) + cp.sum(c[3:4] @ x_c))

# ограничения на диапазон и на общую сумму
cons = [
    x_i[0] <= 2, x_i[1] <= 1, x_i[2] <= 1,
    x_c <= 1,
    ((A[:3] @ x_i) + (A[3:4] @ x_c)) <= b
]

# формирование задачи и ее решение
prb = cp.Problem(obj, cons)
sol = prb.solve(verbose=False, solver='GLPK_MI')
x_opt = np.concatenate([x_i.value, x_c.value])

print(f"Решение: x = {list(np.round(x_opt, 2))}, f = {sol}, {prb.status}")

# output
# Решение: x = [1.0, 1.0, 1.0, 0.5], f = 6.5, optimal

Здесь также наблюдаем полное совпадение решения с предыдущими пакетами.

Заключение

Мы познакомились с несколькими open-source библиотеками для решения оптимизационных задач.

Подводя итог, приведем некоторые рекомендации с чего можно начать, исходя из нашего опыта. Если в задаче NLP первые и вторые производные являются гладкими, то здесь хорошо себя зарекомендовал себя Ipopt в пакете Pyomo. Если в NLP функции не гладкие, но выпуклые, то для решения задачи подойдет CVXPY.
Scipy хорошо справляется с небольшими задачами, а COBYLA можно использовать для задач без градиента целевой функции. Для задачи MILP конкурирующими являются солверы GLPK, Cbc, HiGHS — здесь заранее сказать какой из них будет более производительным для нашей задачи непросто, надо пробовать, перебирать. Scipy с точки зрения интерфейса не очень удобен, если имеются сложные ограничения.

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

Пакеты в Python	Солвер(метод)	NLP	LP	MILP	MINLP
SciPy	COBYLA	+	-	-	-
SciPy	SLSQP	+	-	-	-
SciPy	trust-constr	+	-	-	-
SciPy	HiGHS	+	+	+	-
Pyomo	Ipopt	+	+	-	-
Pyomo, CVXPY	GLPK	-	+	+	-
Pyomo, CVXPY	Cbc	-	+	+	-
CVXPY	Ecos	+-	+	+	-
Pyomo	Bonmin	+	+	+	+

В статье мы акцентировали внимание на применение методов, но совсем не касались вопросов производительности. Об этом и не только мы планируем рассказать в следующей статье.