SciPy, оптимизация с условиями

Ксюша Родионова17.04.2019

0 15 Время чтения: 6 мин.

С SciPy интерактивный сеанс Python превращается в такую же полноценную среду обработки данных, как MATLAB, IDL, Octave, R или SciLab. SciPy (произносится как сай пай) — это основанный на numpy математический пакет, включающий в себя также библиотеки на C и Fortran.

Алгоритмы безусловной оптимизации уже рассмотрены в прошлой статье. В этой статье рассмотрим основные приемы математического программирования — решения задач условной оптимизации для скалярной функции нескольких переменных с помощью пакета scipy.optimize. Более подробную и актуальную справку по функциям scipy всегда можно получить с помощью команды help(), Shift+Tab или в официальной документации.

Введение

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

Подходящий алгоритм оптимизации задается с помощью аргумента функции minimize(..., method="").

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

trust-constr — поиск локального минимума в доверительной области. Статья на wiki, статья на хабре;
SLSQP — последовательное квадратичное программирование с ограничениями, ньютоновский метод решения системы Лагранжа. Статья на вики.
TNC — Truncated Newton Constrained, ограниченное число итераций, хорош для нелинейных функций с большим числом независимых переменных. Статья на wiki.
L-BFGS-B — метод от четверки Broyden–Fletcher–Goldfarb–Shanno, реализованный с уменьшенным потреблением памяти за счет частичной загрузки векторов из матрицы Гессе. Статья на wiki, статья на хабре.
COBYLA — КОБЫЛА Constrained Optimization By Linear Approximation, ограниченная оптимизация с линейной аппроксимацией (без вычисления градиента). Статья на wiki.

В зависимости от выбранного метода, по-разному задаются условия и ограничения для решения задачи:

объектом класса Bounds для методов L-BFGS-B, TNC, SLSQP, trust-constr;
списком (min, max) для этих же методов L-BFGS-B, TNC, SLSQP, trust-constr;
объектом или списком объектов LinearConstraint, NonlinearConstraint для методов COBYLA, SLSQP, trust-constr;
словарем или списком словарей для методов COBYLA, SLSQP.

План статьи:
1) Рассмотреть применение алгоритма условной оптимизации в доверительной области (method="trust-constr") с ограничениями, заданными в виде объектов Bounds, LinearConstraint, NonlinearConstraint ;
2) Рассмотреть последовательное программирование методом наименьших квадратов (method="SLSQP") с ограничениями, заданными в виде словаря {'type', 'fun', 'jac', 'args'};
3) Разобрать пример оптимизации выпускаемой продукции на примере веб-студии.

Условная оптимизация method="trust-constr"

Оба метода реализованы алгоритмами поиска локального минимума в доверительной области и хорошо подходят для крупномасштабных задач. Реализация метода trust-constr основана на EQSQP для задач с ограничениями вида равенства и на TRIP для задач с ограничениями в виде неравенств.

Математическая постановка задачи поиска минимума в общем виде:

$\min_x f(x)$

$c^l \leq c(x) \leq c^u ,$

$x^l \leq x \leq x^u$

Для ограничений строгого равенства нижняя граница устанавливается равной верхней $c^l_j = c^u_j$ .
Для одностороннего ограничения верхняя или нижняя граница устанавливается np.inf с соответствующим знаком.
Пусть необходимо найти минимум известной функцию Розенброка от двух переменных:

$\min_{x_0, x_1} 100(x_0 - x_1^2)^2 + (1-x_0)^2$

При этом заданы следующие ограничения на ее область определения:

$x_0^2 + x_1 \leq 1$

$x_0^2 - x_1 \leq 1$

$ 2x_0 + x_1 = 1 $

$x_0 + 2x_1 \leq 1$

$0 \leq x_0 \leq 1$

5 \leq x_1 \leq 2. $ -0. 0 $

4149, 0. В нашем случае имеется единственное решение в точке $[x_0, x_1] = [0. 5 \leq x_1 \leq 2. 1701]$ , для которой справедливы только первое и четвертое ограничения.
Пройдемся по ограничениям снизу вверх и рассмотрим, как можно их записать в scipy.
Ограничения $0 \leq x_0 \leq 1$ и $-0. 0$ определим с помощью объекта Bounds.

from scipy.optimize import Bounds bounds = Bounds ([0, -0.5], [1.0, 2.0])

Ограничения $x_0 + 2 x_1 \leq 1$ и $2 x_0 + x_1 = 1$ запишем в линейной форме:

$\begin{bmatrix} - \infty \\ 1 \end{bmatrix} \leq \begin{bmatrix} 1 & 2 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \end{bmatrix} \leq \begin{bmatrix} 1 \\ 1 \end{bmatrix}$

Определим эти ограничения в виде объекта LinearConstraint:

import numpy as np from scipy.optimize import LinearConstraint linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])

И наконец нелинейное ограничение в матричной форме:

$c(x) = \begin{bmatrix} x_0^2 + x_1 \\ x_0^2 - x_1 \end{bmatrix} \leq \begin{bmatrix} 1 \\ 1 \end{bmatrix}$

Определим матрицу Якоби для этого ограничения и линейную комбинацию матрицы Гессе с произвольным вектором $v$ :

$J(x) = \begin{bmatrix} 2x_0 & 1 \\ 2x_0 & -1 \end{bmatrix}$

$H(x, v) = \sum\limits_{i=0}^1 v_i \nabla^2 c_i(x) = v_0 \begin{bmatrix} 2 & 0 \\ 2 & 0 \end{bmatrix} + v_1 \begin{bmatrix} 2 & 0 \\ 2 & 0 \end{bmatrix}$

Теперь нелинейное ограничение можем определить как объект NonlinearConstraint:

from scipy.optimize import NonlinearConstraint def cons_f(x): return [x[0]**2 + x[1], x[0]**2 - x[1]] def cons_J(x): return [[2*x[0], 1], [2*x[0], -1]] def cons_H(x, v): return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]]) nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)

Если размер велик, матрицы можно задавать и в разреженном виде:

from scipy.sparse import csc_matrix def cons_H_sparse(x, v): return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]]) nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H_sparse)

или как объект LinearOperator:

from scipy.sparse.linalg import LinearOperator def cons_H_linear_operator(x, v): def matvec(p): return np.array([p[0]*2*(v[0]+v[1]), 0]) return LinearOperator((2, 2), matvec=matvec) nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H_linear_operator)

Доступны следующие стратегии: BFGS и SR1. Когда вычисление матрицы Гессе $H (x, v)$ требует больших затрат, можно использовать класс HessianUpdateStrategy.

from scipy.optimize import BFGS nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())

Гессиан также может быть вычислен с помощью конечных разностей:

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')

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

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())

Решение оптимизационной задачи выглядит следующим образом:

from scipy.optimize import minimize from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prod x0 = np.array([0.5, 0]) res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess, constraints=[linear_constraint, nonlinear_constraint], options={'verbose': 1}, bounds=bounds) print(res.x)

`gtol` termination condition is satisfied. Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.033 s. [0.41494531 0.17010937]

При необходимости, функцию вычисления гессиана можно определить с помощью класса LinearOperator

def rosen_hess_linop(x): def matvec(p): return rosen_hess_prod(x, p) return LinearOperator((2, 2), matvec=matvec) res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop, constraints=[linear_constraint, nonlinear_constraint], options={'verbose': 1}, bounds=bounds) print(res.x)

или произведение Гессиана и произвольного вектора через параметр hessp:

res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_prod, constraints=[linear_constraint, nonlinear_constraint], options={'verbose': 1}, bounds=bounds) print(res.x)

Например, гессиан может быть аппроксимирован с помощью функции SR1 (квази-Ньютоновского приближения). Альтернативно, первая и вторая производные оптимизируемой функции могут быть вычислены приближенно. Градиент может быть аппроксимирован конечными разностями.

from scipy.optimize import SR1 res = minimize(rosen, x0, method='trust-constr', jac="2-point", hess=SR1(), constraints=[linear_constraint, nonlinear_constraint], options={'verbose': 1}, bounds=bounds) print(res.x)

Условная оптимизация method="SLSQP"

Метод SLSQP предназначен для решения задач минимизации функции в виде:

$\min_x f(x)$

$c_j(x) = 0, j \in \mathcal {E}$

$c_j(x) \geq 0, j \in \mathcal {I}$

$lb_i \leq x_i \leq ub_i, i=1,...,N$

$[lb_i, ub_i]$ — множества нижних и верхних границ для области определения функции. Где $\mathcal {E}$ и $\mathcal {I}$ — множества индексов выражений, описывающих ограничения в виде равенств или неравенств.

Линейные и нелинейные ограничения описываются в виде словарей с ключами type, fun и jac.

ineq_cons = {'type': 'ineq', 'fun': lambda x: np.array ([1 - x [0] - 2 * x [1], 1 - x [0] ** 2 - x [1], 1 - x [0] ** 2 + x [1]]), 'jac': lambda x: np.array ([[- 1.0, -2.0], [-2 * x [0], -1.0], [-2 * x [0], 1.0]]) } eq_cons = {'type': 'eq', 'fun': lambda x: np.array ([2 * x [0] + x [1] - 1]), 'jac': lambda x: np.array ([2.0, 1.0]) }

Поиск минимума осуществляется следующим образом:

x0 = np.array([0.5, 0]) res = minimize(rosen, x0, method='SLSQP', jac=rosen_der, constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True}, bounds=bounds) print(res.x)

Optimization terminated successfully. (Exit mode 0) Current function value: 0.34271757499419825 Iterations: 4 Function evaluations: 5 Gradient evaluations: 4 [0.41494475 0.1701105 ]

Пример оптимизации

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

x0 — продающие лэндинги, от 10 т.р.
x1 — корпоративные сайты, от 20 т.р.
x2 — интернет магазины, от 30 т.р.

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

джуны: 4 * 150 = 600 чел * час,
мидлы: 2 * 150 = 300 чел * час,
сеньор: 150 чел * час.

Пусть на разработку и деплой одного сайта типа (x0, x1, x2) первый попавшийся джуниор должен потратить (10, 20, 30) часов, мидл — (7, 15, 20), сеньор — (5, 10, 15) часов лучшего времени своей жизни.

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

def value(x): return - 10*x[0] - 20*x[1] - 30*x[2]

Это не ошибка, при поиске максимума целевая функция минимизируется с обратным знаком.

Следующий шаг — запрещаем перерабатывать своим сотрудникам и вводим ограничения на фонд рабочего времени:

$\begin{bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} \leq \begin{bmatrix} 600 \\ 300 \\ 150 \end{bmatrix}$

Что эквивалентно:

$\begin{bmatrix} 600 \\ 300 \\ 150 \end{bmatrix} - \begin{bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} \geq 0$

ineq_cons = {'type': 'ineq', 'fun': lambda x: np.array ([600 - 10 * x [0] - 20 * x [1] - 30 * x[2], 300 - 7 * x [0] - 15 * x [1] - 20 * x[2], 150 - 5 * x [0] - 10 * x [1] - 15 * x[2]]) }

Формальное ограничение — выпуск продукции должен быть только положительным:

bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])

Мы можем сами выбирать ежемесячные объемы производства продукции, исходя из решения задачи условной оптимизации с scipy.optimize: И наконец самое радужное допущение — из-за низкой цены и высокого качества к нам постоянно выстраивается очередь из довольных клиентов.

x0 = np.array([10, 10, 10]) res = minimize(value, x0, method='SLSQP', constraints=ineq_cons, bounds=bnds) print(res.x)

[7.85714286 5.71428571 3.57142857]

Нестрого округлим до целых и посчитаем месячную загрузку гребцов при оптимальном раскладе продукции x = (8, 6, 3) :

джуны: 8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час;
мидлы: 8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час;
сеньор: 8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.

Сеньор при этом должен пахать не отрываясь от станка, загрузка мидлов составит примерно 2/3, джунов меньше половины. Вывод: чтобы директор получал свой заслуженный максимум, оптимально делать в месяц по 8 лэндингов, 6 средних сайтов и 3 магазина.

Заключение

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

Более хардкорное применение scipy.optimize для построения 3D структуры по набору изображений (статья на хабре) можно посмотреть в scipy-cookbook. Много теории и винрарных примеров можно найти, например, в книге И.Л.Акулича "Математическое программирование в примерах и задачах".

Основной источник информации — docs.scipy.org, желающие поконтрибьютить в перевод этого и других разделов scipy добро пожаловать на GitHub.

Спасибо mephistopheies за участие в подготовке публикации.

Теги