SciPy (vyslovuje sa sai koláč) je matematický balík založený na numpy, ktorý obsahuje aj knižnice C a Fortran. SciPy premení vašu interaktívnu reláciu Pythonu na kompletné prostredie pre vedu o údajoch, ako je MATLAB, IDL, Octave, R alebo SciLab.
V tomto článku sa pozrieme na základné techniky matematického programovania – riešenie problémov podmienenej optimalizácie pre skalárnu funkciu viacerých premenných pomocou balíka scipy.optimize. Neobmedzené optimalizačné algoritmy už boli diskutované v . Podrobnejšiu a aktuálnu pomoc k funkciám scipy môžete vždy získať pomocou príkazu help(), Shift+Tab alebo v .
Úvod
Spoločné rozhranie na riešenie podmienených aj neobmedzených optimalizačných problémov v balíku scipy.optimize poskytuje funkcia minimize(). Je však známe, že univerzálna metóda na riešenie všetkých problémov neexistuje, a tak výber adekvátnej metódy, ako vždy, leží na pleciach výskumníka.
Príslušný optimalizačný algoritmus je špecifikovaný pomocou argumentu funkcie minimize(..., method="").
Pre podmienenú optimalizáciu funkcie viacerých premenných sú k dispozícii implementácie nasledujúcich metód:
trust-constr— hľadať miestne minimum v oblasti spoľahlivosti. , ;SLSQP— sekvenčné kvadratické programovanie s obmedzeniami, newtonovská metóda na riešenie Lagrangeovho systému. .TNC- Skrátený Newton Obmedzený, 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 tímu Broyden–Fletcher–Goldfarb–Shanno, implementovaná so zníženou spotrebou pamäte v dôsledku čiastočného načítania vektorov z Hessovej matice. , .COBYLA— MARE Limited Optimization By Linear Approximation, obmedzená optimalizácia s lineárnou aproximáciou (bez výpočtu gradientu). .
V závislosti od zvolenej metódy sú podmienky a obmedzenia na riešenie problému nastavené inak:
- objekt triedy
Boundsдля методов L-BFGS-B, TNC, SLSQP, trust-constr; - zoznam
(min, max)pre rovnaké metódy L-BFGS-B, TNC, SLSQP, trust-constr; - objekt alebo zoznam objektov
LinearConstraint,NonlinearConstraintpre metódy COBYLA, SLSQP, trust-constr; - slovník alebo zoznam slovníkov
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}pre metódy COBYLA, SLSQP.
Náčrt článku:
1) Zvážte použitie podmieneného optimalizačného algoritmu v oblasti dôveryhodnosti (metóda = “trust-constr”) s obmedzeniami špecifikovanými ako objekty Bounds, LinearConstraint, NonlinearConstraint ;
2) Zvážte sekvenčné programovanie pomocou metódy najmenších štvorcov (metóda = "SLSQP") s obmedzeniami špecifikovanými vo forme slovníka {'type', 'fun', 'jac', 'args'};
3) Analyzujte príklad optimalizácie vyrábaných produktov na príklade webového štúdia.
Metóda podmienenej optimalizácie="trust-constr"
Implementácia metódy trust-constr založené na pre problémy s obmedzeniami formy rovnosti a ďalej pre problémy s obmedzeniami v podobe nerovností. Obidve metódy sú implementované pomocou algoritmov na nájdenie lokálneho minima v oblasti spoľahlivosti a sú vhodné pre rozsiahle problémy.
Matematická formulácia problému hľadania minima vo všeobecnej forme:
Pre prísne obmedzenia rovnosti sa spodná hranica rovná hornej hranici .
Pre jednosmerné obmedzenie sa nastavuje horná alebo dolná hranica np.inf s príslušným znakom.
Nech je potrebné nájsť minimum známej Rosenbrockovej funkcie dvoch premenných:
V tomto prípade sú na jeho definičnú doménu nastavené nasledujúce obmedzenia:
V našom prípade ide o unikátne riešenie , pre ktoré platí len prvé a štvrté obmedzenie.
Prejdime si obmedzenia zdola nahor a pozrime sa, ako ich môžeme napísať v scipy.
Obmedzenie и definujme ho pomocou objektu Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Obmedzenie и Napíšme to v lineárnej forme:
Definujme tieto obmedzenia ako objekt 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 maticovom tvare:
Definujeme Jacobovu maticu pre toto obmedzenie a lineárnu kombináciu Hessovej matice s ľubovoľným vektorom :
Teraz môžeme definovať nelineárne obmedzenie 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ôžu byť špecifikované aj v riedkej forme:
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 predmet 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)Pri výpočte Hessovej matice vyžaduje veľa úsilia, môžete použiť triedu . K dispozícii sú nasledujúce stratégie: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessián možno vypočítať aj pomocou konečných rozdielov:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Jacobovskú maticu pre obmedzenia možno vypočítať aj pomocou konečných rozdielov. V tomto prípade však Hessovu maticu nemožno vypočítať pomocou konečných rozdielov. Hessian 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á takto:
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]V prípade potreby je možné definovať funkciu pre výpočet Hessian 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 súčin Hessianu 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 možno prvú a druhú deriváciu funkcie, ktorá sa má optimalizovať, aproximovať. Napríklad Hessian môže byť aproximovaný pomocou funkcie SR1 (kvázi-newtonovská aproximácia). Gradient sa dá aproximovať konečnými rozdielmi.
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)Metóda podmienenej optimalizácie="SLSQP"
Metóda SLSQP je určená na riešenie problémov minimalizácie funkcie v tvare:
kde и — súbory indexov výrazov opisujúcich obmedzenia vo forme rovnosti alebo nerovností. — množiny dolných a horných hraníc pre definičný obor funkcie.
Lineárne a nelineárne obmedzenia sú opísané vo forme slovníkov s kľúčmi 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])
}Hľadanie minima sa vykonáva takto:
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 ]Príklad optimalizácie
V súvislosti s prechodom na piatu technologickú štruktúru sa pozrime na optimalizáciu výroby na príklade webového štúdia, ktoré nám prináša malý, ale stabilný príjem. Predstavme si seba ako riaditeľa kuchyne, ktorá vyrába tri druhy produktov:
- x0 - predaj vstupných stránok, od 10 tr.
- x1 - firemné stránky, od 20 tr.
- x2 - internetové obchody, od 30 tr.
Náš priateľský pracovný tím tvoria štyria juniori, dvaja strední a jeden senior. Ich mesačný fond pracovného času:
- júna:
4 * 150 = 600 чел * час, - stredy:
2 * 150 = 300 чел * час, - seňor:
150 чел * час.
Nechajte prvého dostupného juniora stráviť (0, 1, 2) hodín vývojom a nasadením jedného webu typu (x10, x20, x30), stredného - (7, 15, 20), staršieho - (5, 10, 15 ) hodiny najlepšieho obdobia vášho života.
Ako každý normálny riaditeľ, aj my chceme maximalizovať mesačné zisky. Prvým krokom k úspechu je zapísať si účelovú funkciu value ako výška príjmu z produktov vyrobených za mesiac:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Nejde o chybu, pri hľadaní maxima sa účelová funkcia minimalizuje s opačným znamienkom.
Ďalším krokom je zakázať našim zamestnancom prepracovanie a zaviesť obmedzenia 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álnym obmedzením je, že výstup produktu musí byť iba pozitívny:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])A nakoniec najružovejší predpoklad je, že kvôli nízkej cene a vysokej kvalite sa na nás neustále tlačí rad spokojných zákazníkov. Mesačné objemy výroby si môžeme zvoliť sami na základe riešenia problému s obmedzenou optimalizáciou 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]Voľne zaokrúhlime na celé čísla a vypočítame mesačnú záťaž veslárov s optimálnym rozložením produktov x = (8, 6, 3) :
- júna:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - stredy:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - seňor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Záver: na to, aby riaditeľ dostal svoje zaslúžené maximum, je optimálne vytvoriť 8 cieľových stránok, 6 stredne veľkých stránok a 3 obchody mesačne. V tomto prípade musí senior orať bez vzhliadnutia od stroja, zaťaženie stredov bude približne 2/3, juniorov menej ako polovicu.
Záver
V článku sú načrtnuté základné techniky práce s balíkom scipy.optimize, ktorý sa používa na riešenie problémov podmienenej minimalizácie. Osobne používam scipy čisto na akademické účely, a preto má uvedený príklad taký komický charakter.
Veľa teórie a virtuálnych príkladov možno nájsť napríklad v knihe I. L. Akulicha „Matematické programovanie v príkladoch a problémoch“. Viac hardcore aplikácie scipy.optimize vytvoriť 3D štruktúru zo sady obrázkov () si môžete pozrieť v .
Hlavným zdrojom informácií je tých, ktorí chcú prispieť k prekladu tejto a iných častí scipy Vitajte v .
Vďaka za spoluúčasť na príprave publikácie.
Zdroj: hab.com