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