SciPy (vyslovuje sa ako saj paj) je matematický balík založený na numpy, ktorý zahŕňa aj knižnice v C a Fortran. So SciPy sa interaktívna relácia Pythonu mení na rovnako plnohodnotné prostredie spracovania dát ako MATLAB, IDL, Octave, R alebo SciLab.
V tomto článku sa pozrieme na hlavné techniky matematického programovania — riešenie úloh podmieneného optimalizovania pre skalarni funkciu viacerých premenných pomocou balíka scipy.optimize. Algoritmy bezpodmienečného optimalizovania už boli preskúmané v . Podrobnejšie a aktuálne informácie o funkciách scipy sú vždy dostupné pomocou príkazu help(), Shift+Tab alebo v .
Úvod
Všeobecné rozhranie pre riešenie úloh ako podmieneného, tak aj bezpodmienečného optimalizovania v balíku scipy.optimize poskytuje funkcia minimize(). Je však známe, že neexistuje univerzálny spôsob na riešenie všetkých úloh, preto voľba adekvátnej metódy ako vždy spočíva na pleciach výskumníka.
Príslušný algoritmus optimalizácie je určený pomocou argumentu funkcie minimize(..., method="").
Pre podmienené optimalizovanie funkcie viacerých premenných sú dostupné implementácie nasledujúcich metód:
trust-constr— hľadanie lokálneho minima v dôveryhodnej oblasti. , ;SLSQP— sekvenčné kvadratické programovanie s obmedzeniami, Newtonova metóda riešenia systému Lagrangeových podmienok. .TNC— Truncated Newton Constrained, obmedzený počet iterácií, vhodný pre nelineárne funkcie s veľkým počtom nezávislých premenných. .L-BFGS-B— metóda od štvorice Broyden–Fletcher–Goldfarb–Shanno, realizovaná s nižšou spotrebou pamäte vďaka čiastočnému načítaniu vektorov z Hessovej matice. , .COBYLA— KOBYLA (Constrained Optimization By Linear Approximation), obmedzená optimalizácia s lineárnou aproximáciou (bez výpočtu gradientu). .
V závislosti od zvolenej metódy sa odlišne zadávajú podmienky a obmedzenia pre riešenie úlohy:
- objektom triedy
Boundspre metódy L-BFGS-B, TNC, SLSQP, trust-constr; - zoznamom
(min, max)pre tieto metódy L-BFGS-B, TNC, SLSQP, trust-constr; - objektom alebo zoznamom objektov
LinearConstraint,NonlinearConstraintpre metódy COBYLA, SLSQP, trust-constr; - súborom alebo zoznamom súborov
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}pre metódy COBYLA, SLSQP.
Plán článku:
1) Zvážiť aplikáciu algoritmu podmieneného optimalizovania v dôveryhodnej oblasti (method="trust-constr") s obmedzeniami zadanými v podobe objektov. Bounds, LinearConstraint, NonlinearConstraint ;
2) Zvážiť sekvenčné programovanie metódou najmenších štvorcov (method=»SLSQP») s obmedzeniami zadanými vo forme slovníka {'type', 'fun', 'jac', 'args'};
3) Rozobrať príklad optimalizácie vyrábanej produkcie na príklade webovej štúdie.
Podmienená optimalizácia method=»trust-constr»
Realizácia metódy trust-constr je založená na pre úlohy s obmedzeniami vo forme rovnosti a na pre úlohy s obmedzeniami vo forme nerovností. Obe metódy sú realizované algoritmami hľadania lokálneho minima v dôveryhodnej oblasti a sú dobre prispôsobené pre veľkokapacitné úlohy.
Matematická formulácia úlohy hľadania minima vo všeobecnom prevedení:
Pre obmedzenia prísnej rovnosti sa dolná hranica nastavuje rovná hornej 
.
Pre jednostranné obmedzenie sa horná alebo dolná hranica nastavuje np.inf s príslušným znamienkom.
Predpokladajme, že je potrebné nájsť minimum známej funkcie Rosenbrok od dvoch premenných:
Pri tom sú zadané nasledujúce obmedzenia na jej oblasť definovania:
V našom prípade existuje jediné riešenie v bode 
, pre ktorý platia iba prvé a štvrté obmedzenia.
Pôjdeme po obmedzeniach zdola nahor a zvážime, ako ich môžeme zapísať v scipy.
Obmedzenia 
a 
určíme pomocou objektu Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Obmedzenia 
a 
zapíšeme v lineárnej forme:
Určíme tieto obmedzenia vo forme objektu LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])A nakoniec nelineárne obmedzenie v maticovej forme:
Určíme maticu Jakobi pre toto obmedzenie a lineárnu kombináciu matice Hessiana s ľubovoľným vektorom 
:
Teraz môžeme nelineárne obmedzenie definovať ako objekt 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)Ak je veľkosť veľká, matice môžeme zadávať aj v riedkom formáte:
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)alebo ako objekt 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)Keď výpočet Hesseovej matice 
vyžaduje veľké náklady, môže sa použiť trieda . Dostupné sú nasledujúce stratégie: BFGS a SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hesseova matica môže byť tiež vypočítaná pomocou konečných diferencie:
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess='2-point')Matica Jacobian pre obmedzenia môže byť taktiež vypočítaná pomocou konečných diferencie. Avšak v tomto prípade Hesseovu maticu pomocou konečných diferencie už nie je možné vypočítať. Hesseova matica musí byť definovaná ako funkcia alebo pomocou triedy HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess=BFGS())Riešenie optimalizačného problému vyzerá nasledovne:
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` podmienka ukončenia je splnená.
Počet iterácií: 12, vyhodnotenia funkcie: 8, CG iterácie: 7, optimálnosť: 2.99e-09, porušenie obmedzenia: 1.11e-16, doba vykonávania: 0.033 s.
[0.41494531 0.17010937]Ak je to potrebné, funkciu na výpočet Hesseovej matice je možné definovať pomocou triedy 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)alebo násobenie Hesseovej matice a ľubovoľného vektora cez parameter 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)Alternatívne môžu byť prvé a druhé derivácie optimalizovanej funkcie vypočítané približne. Napríklad, Hesseova matica môže byť aproximovaná pomocou funkcie SR1 (kvazi-Newtonovskej aproximácie). Gradient môže byť aproximovaný konečnými diferenciami.
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)Podmienená optimalizácia method='SLSQP'
Metóda SLSQP je navrhnutá na riešenie úloh minimalizácie funkcie vo forme:
Kde 
a 
— množstva indexov výrazov, ktoré popisujú obmedzenia vo forme rovností alebo nerovností. 
— množstva dolných a horných hraníc pre oblasť definície funkcie.
Lineárne a nelineárne obmedzenia sú opísané vo formáte slovníkov s kľúčmi type, zábava a 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])
}Hľadanie minima sa uskutočňuje nasledujúcim spôsobom:
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)Optimalizácia bola úspešne ukončená. (Exit mode 0)
Aktuálna hodnota funkcie: 0.34271757499419825
Iterácie: 4
Vyhodnotenia funkcie: 5
Vyhodnotenia gradientu: 4
[0.41494475 0.1701105 ]Príklad optimalizácie
V súvislosti s prechodom na piaty technologický spôsob, zamerajme sa na optimalizáciu výroby na príklade webovej štúdie, ktorá nám prináša malý, ale stabilný príjem. Predstavme si, že sme riaditeľom galérie, v ktorej sa vyrábajú tri druhy produktov:
- x0 — predajné landing pages, od 10 t.r.
- x1 — firemné weby, od 20 t.r.
- x2 — internetové obchody, od 30 t.r.
Náš priateľský pracovný kolektív sa skladá zo štyroch juniorov, dvoch midlov a jedného seniora. Fond ich pracovného času za mesiac:
- juniori:
4 * 150 = 600 h * osoba, - midli:
2 * 150 = 300 h * osoba, - senior:
150 h * osoba.
Nech na vývoj a nasadenie jedného typu webu (x0, x1, x2) prvý najbližší junior musí stráviť (10, 20, 30) hodín, midl — (7, 15, 20), senior — (5, 10, 15) hodín najlepšieho času svojho života.
Ako každému normálnemu riaditeľovi, aj nám záleží na maximálnych mesačných ziskoch. Prvý krok k úspechu — zapíšeme cieľovú funkciu value ako súčet príjmov z výroby za mesiac:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]To nie je chyba, pri hľadaní maxima je cieľová funkcia minimalizovaná s opačným znamienkom.
Ďalší krok — zakážeme svojim zamestnancom nadčasy a zavádzame obmedzenia na fond pracovného času:
Čo je ekvivalentné:
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]])
}Formálne obmedzenie — výroba produktov musí byť iba kladná:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])A nakoniec najoptimistickejšia možnosť — kvôli nízkej cene a vysokej kvalite sa u nás neustále vytvára fronta spokojných klientov. Môžeme si sami vyberajú mesačné objemy výroby produktov, na základe riešenia úlohy podmieneného optimálneho. 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]Neprísne zaokrúhlime na celé čísla a spočítame mesačné zaťaženie veslovateľov pri optimálnom rozložení produktov. x = (8, 6, 3) :
- juniori:
8 * 10 + 6 * 20 + 3 * 30 = 290 os.*hod.; - midli:
8 * 7 + 6 * 15 + 3 * 20 = 206 os.*hod.; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 os.*hod..
Výsledok: aby riaditeľ dostal svoj zaslúžený maximum, optimálne je robiť mesačne 8 landing stránok, 6 stredných webov a 3 obchody. Senior by mal prácu vykonávať bez prestávky, zaťaženie midlov bude asi 2/3, juniorov menej než polovicu.
Záver
V článku sú uvedené základné techniky práce s balíkom scipy.optimize, ktoré sa používajú na riešenie úloh podmienenej minimalizácie. Osobne používam scipy čisto na akademické účely, preto uvedený príklad má vtipný charakter.
Mnoho teórie a skvelých príkladov nájdete napríklad v knihe I.L. Akulicha «Matematické programovanie v príkladoch a úlohách». Pre hardcore použitie scipy.optimize na vytvorenie 3D štruktúry podľa súboru obrázkov () si môžete pozrieť v .
Hlavný zdroj informácií — , tí, ktorí sa chcú zapojiť do prekladu tejto a iných sekcií, sú vítaní na scipy mephistopheies .
Ďakujem SciPy (vyslovuje sa ako saupai) — to je založené na.
Zdroj: habr.com
