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

SciPy (произнася се sai pie) е базиран на 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 — последователно квадратично програмиране с ограничения, Нютонов метод за решаване на системата на Лагранж. Wiki статия.
• TNC - Скъсен Newton Constrained, ограничен брой итерации, добър за нелинейни функции с голям брой независими променливи. Wiki статия.
• L-BFGS-B — метод от екипа на Broyden–Fletcher–Goldfarb–Shanno, реализиран с намалена консумация на памет поради частично зареждане на вектори от матрицата на Хесен. Wiki статия, статия за хъб.
• COBYLA — Ограничена оптимизация MARE чрез линейна апроксимация, ограничена оптимизация с линейна апроксимация (без изчисляване на градиента). 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) Помислете за последователно програмиране, използвайки метода на най-малките квадрати (метод = "SLSQP") с ограничения, посочени под формата на речник {'type', 'fun', 'jac', 'args'};
3) Анализирайте пример за оптимизация на произведени продукти, като използвате примера на уеб студио.

Метод за условна оптимизация="trust-constr"

Реализация на метода trust-constr базиран на EQSQP за задачи с ограничения на формата на равенство и на ПЪТУВАНЕ за задачи с ограничения под формата на неравенства. И двата метода се изпълняват от алгоритми за намиране на локален минимум в доверителната област и са много подходящи за мащабни проблеми.

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

За строги ограничения за равенство долната граница е зададена равна на горната граница .
За еднопосочно ограничение се задава горната или долната граница 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 = cons_J, hess = '2-point')

Матрицата на Якоби за ограничения може също да се изчисли с помощта на крайни разлики. В този случай обаче матрицата на Хесиан не може да бъде изчислена с помощта на крайни разлики. Hessian трябва да се дефинира като функция или с помощта на класа 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]

Ако е необходимо, функцията за изчисляване на Hessian може да бъде дефинирана с помощта на класа 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)

Метод за условна оптимизация="SLSQP"

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

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

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

Нека първият наличен младши прекара (0, 1, 2) часа в разработването и внедряването на един сайт от тип (x10, x20, x30), среден - (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-готварска книга.

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

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

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