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 — метод от четворката 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;
• речник или списък от речници {'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:

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)

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

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='2-point', 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` условието за прекратяване е удовлетворено.
Брой итерации: 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 (квази-Нютоново апроксиране). Градиентът може да бъде апроксимирани с краен диференциал.

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

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

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

Оптимизацията приключи успешно. (Exit mode 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) часа от най-доброто време на живота си.

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

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 чисто за академични цели, затова посоченият пример има хумористичен характер.

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

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

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

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