SciPy (vyslovuje se jako saj paj) — to je matematický balíček založený na numpy, který zahrnuje také knihovny v C a Fortranu. S SciPy se interaktivní sezení Pythonu promění v plnohodnotné prostředí pro zpracování dat, jako jsou MATLAB, IDL, Octave, R nebo SciLab.
V tomto článku se zaměříme na základní techniky matematického programování — řešení úloh podmíněné optimalizace pro skalarovou funkci několika proměnných pomocí balíčku scipy.optimize. Algoritmy bezpodmínkové optimalizace byly již probrány v . Podrobnější a aktuální nápovědu k funkcím scipy je vždy možné získat pomocí příkazu help(), Shift+Tab nebo v .
Úvod
Obecné rozhraní pro řešení úloh jak podmíněné, tak bezpodmínkové optimalizace v balíčku scipy.optimize je poskytováno funkcí minimize(). Je však známo, že univerzální způsob pro řešení všech úloh neexistuje, a proto volba adekvátní metody jako vždy závisí na ramenou výzkumníka.
Vhodný optimalizační algoritmus se určuje pomocí argumentu funkce minimize(..., method="").
Pro podmíněnou optimalizaci funkce několika proměnných jsou k dispozici implementace následujících metod:
trust-constr— hledání lokálního minima v důvěryhodném oblasti. , ;SLSQP— sekvenční kvadratické programování s omezeními, Newtonova metoda pro řešení systému Lagrange. .TNC— Truncated Newton Constrained, omezený počet iterací, dobrý pro nelineární funkce s velkým množstvím nezávislých proměnných. .L-BFGS-B— metoda čtveřice Broyden–Fletcher–Goldfarb–Shanno, realizovaná s menší spotřebou paměti díky částečnému načítání vektorů z matice Hessovy. , .COBYLA— KOBYLA Constrained Optimization By Linear Approximation, omezená optimalizace s lineární aproximací (bez výpočtu gradientu). .
V závislosti na vybrané metodě se podmínky a omezení pro řešení úlohy stanovují různě:
- objektem třídy
Boundspro metody L-BFGS-B, TNC, SLSQP, trust-constr; - seznamem
(min, max)pro tyto metody L-BFGS-B, TNC, SLSQP, trust-constr; - objektem nebo seznamem objektů
LinearConstraint,NonlinearConstraintpro metody COBYLA, SLSQP, trust-constr; - slovníkem nebo seznamem slovníků
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}pro metody COBYLA, SLSQP.
Plán článku:
1) Zvážit použití algoritmu podmíněné optimalizace ve svěřené oblasti (method=»trust-constr») s omezeními stanovenými ve formě objektů Bounds, LinearConstraint, NonlinearConstraint ;
2) Zvážit sekvenční programování metodou nejmenších čtverců (method=»SLSQP») s omezeními stanovenými ve formě slovníku {'type', 'fun', 'jac', 'args'};
3) Prozkoumat příklad optimalizace vyráběného produktu na příkladu webového studia.
Podmíněná optimalizace method=»trust-constr»
Implementace metody trust-constr je založena na pro úkoly s omezeními ve formě rovnosti a na pro úkoly s omezeními ve formě nerovností. Obě metody jsou realizovány algoritmy hledání lokálního minima ve svěřené oblasti a dobře se hodí pro rozsáhlé úkoly.
Matematická formulace úlohy hledání minima v obecném vyjádření:
Pro omezení přísné rovnosti je dolní mez stanovena na hodnotu horní meze 
.
Pro jednostranné omezení je stanovena horní nebo dolní mez np.inf se odpovídajícím znaménkem.
Předpokládejme, že je třeba nalézt minimum známé funkce Rosenbrok od dvou proměnných:
Přitom jsou dána následující omezení na její definiční oblast:
V našem případě existuje jediné řešení v bodě 
, pro nějž platí pouze první a čtvrté omezení.
Projdeme si omezeními zdola vzhůru a zvážíme, jak je lze zapsat v scipy.
Omezení 
a 
určíme pomocí objektu Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Omezení 
a 
zapíšeme v lineární formě:
Určíme tato omezení ve formě objektu LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])A konečně nelineární omezení v maticové formě:
Určíme Jacobovu matici pro toto omezení a lineární kombinaci Hessovy matice s libovolným vektorem 
:
Nyní můžeme nelineární omezení definovat jako 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)Pokud je velikost velká, matice mohou být definovány také v řídkém formátu:
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)nebo jako objekt LinearOperator:
z 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)Když výpočet Hessovy matice 
vyžaduje velké náklady, lze použít třídu . K dispozici jsou následující strategie: BFGS a SR1.
z scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian lze také vypočítat pomocí konečných diferencí:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Jacobiho matice pro omezení může být také vypočítána pomocí konečných diferencí. V tomto případě však nelze Hessovu matici vypočítat pomocí konečných diferencí. Hessian musí být definován jako funkce nebo pomocí třídy HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Řešení optimalizačního úkolu vypadá následovně:
z scipy.optimize import minimize
z 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` terminační podmínka je splněna.
Počet iterací: 12, vyhodnocení funkcí: 8, CG iterací: 7, optimálnost: 2.99e-09, porušení omezení: 1.11e-16, doba provádění: 0.033 s.
[0.41494531 0.17010937]V případě potřeby lze funkci pro výpočet Hessovy matice definovat pomocí třídy 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)nebo produkt Hessovy matice a libovolného vektoru prostřednictvím parametru 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)Alternativně lze první a druhé derivace optimalizované funkce přibližně vypočítat. Například Hessian může být aproximován pomocí funkce SR1 (kvazi-Newtonovského přístupu). Gradient může být aproximován konečnými diferencemi.
z 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)Podmíněná optimalizace method='SLSQP'
Metoda SLSQP je určena k řešení problémů minimalizace funkce ve formě:
Kde 
a 
— množství indexů výrazů popisující omezení ve formě rovnic nebo nerovnic. 
— množství dolních a horních mezí pro definiční obor funkce.
Lineární a nelineární omezení jsou popsána ve formě slovníků s klíči 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])
}Hledání minima probíhá následujícím způsobem:
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)Optimalizace byla úspěšně dokončena. (Režim odchodu 0)
Aktuální hodnota funkce: 0.34271757499419825
Iterace: 4
Vyhodnocení funkce: 5
Vyhodnocení gradientu: 4
[0.41494475 0.1701105 ]Příklad optimalizace
V souvislosti s přechodem na pátou technologickou strukturu, se podíváme na optimalizaci výroby na příkladu webové studia, které přináší malý, ale stabilní příjem. Představme si, že jsme ředitelem galerie, která vyrábí tři druhy produktů:
- x0 — prodejní landing pages, od 10 tis. Kč.
- x1 — firemní webové stránky, od 20 tis. Kč.
- x2 — internetové obchody, od 30 tis. Kč.
Náš přátelský pracovní kolektiv se skládá ze čtyř juniorů, dvou midlů a jednoho seniora. Fond jejich pracovní doby za měsíc:
- junioři:
4 * 150 = 600 osobohodin, - mídla:
2 * 150 = 300 osobohodin, - senior:
150 osobohodin.
Předpokládejme, že na vývoj a nasazení jednoho webu typu (x0, x1, x2) musí první náhodný junior strávit (10, 20, 30) hodin, midl — (7, 15, 20), senior — (5, 10, 15) hodin svého nejlepšího času.
Jako každému normálnímu řediteli nám záleží na maximalizaci měsíčního zisku. Prvním krokem k úspěchu je zapsání cílové funkce value jako součet příjmů z výrobků vyrobených za měsíc:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]To není chyba, při hledání maxima je cílová funkce minimalizována s opačným znaménkem.
Dalším krokem je zakázat našim zaměstnancům přesčasy a zavést omezení na fond pracovní doby:
Co 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ální omezení — produkce by měla být pouze pozitivní:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])A nakonec to nejoptimističtější předpoklady — kvůli nízké ceně a vysoké kvalitě k nám neustále čeká fronta spokojených zákazníků. Můžeme si volit měsíční objemy produkce na základě řešení úkolu podmínkového optimalizace s 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]Nepřísně zaokroulíme na celé číslo a spočítáme měsíční zátěž veslařů při optimálním uspořádání produktů x = (8, 6, 3) :
- junioři:
8 * 10 + 6 * 20 + 3 * 30 = 290 osob * hodin; - mídla:
8 * 7 + 6 * 15 + 3 * 20 = 206 osob * hodin; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 osob * hodin.
Závěr: aby ředitel dostával svůj zasloužený maximum, je optimální dělat měsíčně 8 landingů, 6 středních stránek a 3 obchody. Senior by měl pracovat nepřetržitě, zátěž middle by měla být přibližně 2/3, juniorů méně než polovina.
Závěr
V článku jsou nastíněny základní techniky práce s balíčkem scipy.optimize, používané k řešení úkolů podmínkového minimalizování. Osobně používám scipy čistě pro akademické účely, proto má uvedený příklad takový humorný charakter.
Mnoho teorie a raritních příkladů lze nalézt například v knize I.L.Akuliče „Matematické programování v příkladech a úlohách“. Více hardcore aplikace scipy.optimize pro stavbu 3D struktury na základě souboru obrázků () lze vidět v .
Hlavní zdroj informací — , ti, kteří chtějí přispět k překladu této a dalších částí, scipy vítejte na .
Děkuji za účast na přípravě publikace.
Zdroj: habr.com
