SciPy (вымаўляецца як сай пай) - гэта заснаваны на numpy матэматычны пакет, які ўключае ў сябе таксама бібліятэкі на C і Fortran. З SciPy інтэрактыўны сеанс Python ператвараецца ў такое ж паўнавартаснае асяроддзе апрацоўкі дадзеных, як MATLAB, IDL, Octave, R ці SciLab.
У гэтым артыкуле разгледзім асноўныя прыёмы матэматычнага праграмавання - рашэнні задач умоўнай аптымізацыі для скалярнай функцыі некалькіх зменных з дапамогай пакета scipy.optimize. Алгарытмы безумоўнай аптымізацыі ўжо разгледжаны ў . Больш падрабязную і актуальную даведку па функцый scipy заўсёды можна атрымаць з дапамогай каманды help(), Shift+Tab або ў .
Увядзенне
Агульны інтэрфейс для рашэння задач як умоўнай, так і безумоўнай аптымізацыі ў пакеце scipy.optimize падаецца функцыяй minimize(). Аднак вядома, што ўніверсальнага спосабу для рашэння ўсіх задач не існуе, таму выбар адэкватнага метаду як заўсёды кладзецца на плечы даследніка.
Падыходны алгарытм аптымізацыі задаецца з дапамогай аргументу функцыі minimize(..., method="").
Для ўмоўнай аптымізацыі функцыі некалькіх зменных даступныя рэалізацыі наступных метадаў:
trust-constr- Пошук лакальнага мінімуму ў давернай вобласці. , ;SLSQP- паслядоўнае квадратычнае праграмаванне з абмежаваннямі, ньютанаўскі метад рашэння сістэмы Лагранжа. .TNC- Truncated Newton Constrained, абмежаваны лік ітэрацый, добры для нелінейных функцый з вялікім лікам незалежных зменных. .L-BFGS-B- метад ад чацвёркі Broyden-Fletcher-Goldfarb-Shanno, рэалізаваны з паменшаным спажываннем памяці за кошт частковай загрузкі вектараў з матрыцы Гесэ. , .COBYLA- Кабыла Constrained Optimization By Linear Approximation, абмежаваная аптымізацыя з лінейнай апраксімацыяй (без вылічэння градыенту). .
У залежнасці ад абранага метаду, па-рознаму задаюцца ўмовы і абмежаванні для рашэння задачы:
- аб'ектам класа
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 заснавана на для задач з абмежаваннямі віду роўнасці і на для задач з абмежаваннямі ў выглядзе няроўнасцей. Абодва метаду рэалізаваны алгарытмамі пошуку лакальнага мінімуму ў давернай вобласці і добра падыходзяць для буйнамаштабных задач.
Матэматычная пастаноўка задачы пошуку мінімуму ў агульным выглядзе:
Для абмежаванняў строгай роўнасці ніжняя мяжа ўсталёўваецца роўнай верхняй .
Для аднабаковага абмежавання верхняя ці ніжняя мяжа усталёўваецца 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)Калі вылічэнне матрыцы Гесэ патрабуе вялікіх выдаткаў, можна выкарыстоўваць клас . Даступныя наступныя стратэгіі: 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')Матрыцу Якобі для абмежаванняў таксама можна вылічыць з дапамогай канчатковых рознасцяў. Аднак, у гэтым выпадку матрыцу Гесэ з дапамогай канчатковых рознасцяў ужо не вылічыць. Гесіян павінен быць вызначаны ў выглядзе функцыі або з дапамогай класа 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 прызначаны для рашэння задач мінімізацыі функцыі ў выглядзе:
Дзе и - Мноствы індэксаў выразаў, якія апісваюць абмежаванні ў выглядзе роўнасцяў або няроўнасцяў. - Мноствы ніжніх і верхніх межаў для вобласці вызначэння функцыі.
Лінейныя і нелінейныя абмежаванні апісваюцца ў выглядзе слоўнікаў з ключамі 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]Гэта не памылка, пры пошуку максімуму мэтавая функцыя мінімізуецца са зваротным знакам.
Наступны крок - забараняем перапрацоўваць сваім супрацоўнікам і ўводзім абмежаванні на фонд працоўнага часу:
Што эквівалентна:
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 сардэчна запрашаем на .
Дзякуй за ўдзел у падрыхтоўцы публікацыі.
Крыніца: habr.com