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

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

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

Въведение

Общият интерфейс за решаване на задачи както за условна, така и за безусловна оптимизация в пакета scipy.optimize се предоставя чрез функцията minimize(). Въпреки това е известно, че не съществува универсален начин за решаване на всички задачи, затова изборът на адекватен метод отново е в ръцете на изследователя.
Подходящият алгоритъм за оптимизация се задава с помощта на аргумента на функцията minimize(..., method="").
За условна оптимизация на функция на множество променливи са налични реализации на следните методи:

• trust-constr — търсене на локален минимум в доверителна област. Статия на wiki, статия на хабре;
• SLSQP — последователно квадратно програмиране с ограничения, Нютонов метод за решаване на системата на Лагранж. Статия на вики.
• TNC — Truncated Newton Constrained, ограничено число итерации, добър е за нелинейни функции с много независими променливи. Статия на wiki.
• L-BFGS-B — метод от четворката на Бройдън–Флетчър–Голдфарб–Шанно, реализиран с намалено потребление на памет с помощта на частично зареждане на вектори от матрицата на Хес. Статия на 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;
• словарь или списък от речници {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} за методите COBYLA, SLSQP.

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

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

Реализация на метода trust-constr е основана на EQSQP за задачи с ограничения от вида равенство и на TRIP за задачи с ограничения под формата на неравенства. И двата метода са реализирани с алгоритми за търсене на локален минимум в доверителната област и са подходящи за мащабни задачи.

Математическа формулировка на задачата за намиране на минимум в общ вид:

За ограничения при стриктно равенство долната граница се установява равна на горната .
За едностранно ограничение горната или долната граница се определя np.inf с съответния знак.
Нека е необходимо да се намери минимум на известната функция на Розенброк от две променливи:

При това са зададени следните ограничения на нейното определение:

В нашия случай има единствено решение в точка , за която важат само първото и четвъртото ограничение.
Нека разгледаме ограниченията отдолу нагоре и видим как можем да ги запишем в scipy.
Ограничения и да определим с помощта на обекта Bounds.

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

Ограничения и да запишем в линейна форма:

Да определим тези ограничения под формата на обект LinearConstraint:

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

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

Да определим Якобиевата матрица за това ограничение и линейна комбинация от Хесовата матрица с произволен вектор :

Сега нелинейното ограничение можем да определим като обект 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:

от scipy.sparse.linalg импортировать 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)

непрерывное ограничение = Нелинейное ограничение(cons_f, -np.inf, 1,
                                jac=cons_J, hess=cons_H_linear_operator)

Когато изчисляването на Хесовата матрица изисква големи разходи, може да се използва класът HessianUpdateStrategy. Доступни са следните стратегии: BFGS и SR1.

от scipy.optimize импорт BFGS

непрерывное ограничение = Нелинейно ограничение(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())

Хесовата матрица може също да бъде изчислена с помощта на крайни разности:

непрерывное ограничение = Нелинейно ограничение(cons_f, -np.inf, 1, jac='2-point', hess='2-point')

Яковата матрица за ограниченията може също да бъде изчислена с помощта на крайни разности. Въпреки това, в този случай Хесовата матрица не може да бъде изчислена с крайни разности. Хесовата матрица трябва да бъде дефинирана като функция или с помощта на клас HessianUpdateStrategy.

непрерывное ограничение = Нелинейно ограничение(cons_f, -np.inf, 1, jac='2-point', hess=BFGS())

Решението на оптимизационната задача изглежда по следния начин:

от scipy.optimize импорт minimize
от scipy.optimize импорт 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` терминалното условие е удовлетворено.
Брой итерации: 12, оценки на функцията: 8, CG итерации: 7, оптималност: 2.99e-09, нарушение на условията: 1.11e-16, време на изпълнение: 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 (квази-Нютоново апроксимиране). Градиентът може да бъде апроксимирани с крайни разности.

от scipy.optimize импорт 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 е предназначен за решаване на задачи за минимизиране на функция под формата на:

Къде и — множества индекси на изрази, описващи ограничения под формата на равенства или неравенства. — множество долни и горни граници за областта на дефиниция на функцията.

Линейните и нелинейните ограничения се описват под формата на речници с ключове тип, забавление и 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)

Оптимизацията беше завършена успешно.    (Режим на изход 0)
            Текуща стойност на функцията: 0.34271757499419825
            Итерации: 4
            Оценки на функцията: 5
            Оценки на градиента: 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]

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

Следващата стъпка — забраняваме на служителите да работят над часовете и въвеждаме ограничения на фонда за работно време:

Което е еквивалентно на:

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 човека * час.

Извод: за да получи директорът своя заслужен максимум, е оптимално да правим по 8 лендинг страници, 6 средни сайта и 3 магазина на месец. Сеньорът при това трябва да работи упорито, натоварването на мидлите ще е около 2/3, а на джуновете — под половината.

Заключение

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

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

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

Благодаря mephistopheies за участие в подготовката на публикацията.

Източник: habr.com