SciPy (vyslovováno sai pie) je matematický balíček založený na numpy, který také obsahuje knihovny C a Fortran. SciPy promění vaši interaktivní relaci Pythonu na kompletní prostředí pro vědu o datech, jako je MATLAB, IDL, Octave, R nebo SciLab.
V tomto článku se podíváme na základní techniky matematického programování - řešení podmíněných optimalizačních problémů pro skalární funkci několika proměnných pomocí balíčku scipy.optimize. O neomezených optimalizačních algoritmech již byla řeč . Podrobnější a aktuální nápovědu k funkcím scipy lze vždy získat pomocí příkazu help(), Shift+Tab nebo v .
úvod
Společné rozhraní pro řešení podmíněných i neomezených optimalizačních problémů v balíčku scipy.optimize poskytuje funkce minimize(). Je však známo, že neexistuje univerzální metoda pro řešení všech problémů, takže výběr adekvátní metody jako vždy leží na bedrech výzkumníka.
Vhodný optimalizační algoritmus je specifikován 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— hledat místní minimum v oblasti spolehlivosti. , ;SLSQP— sekvenční kvadratické programování s omezeními, Newtonova metoda pro řešení Lagrangeova systému. .TNC- Zkrácený Newton Constrained, omezený počet iterací, vhodný pro nelineární funkce s velkým počtem nezávislých proměnných. .L-BFGS-B— metoda od týmu Broyden–Fletcher–Goldfarb–Shanno, implementovaná se sníženou spotřebou paměti díky částečnému načítání vektorů z hessovské matice. , .COBYLA— MARE Constrained Optimization By Linear Approximation, omezená optimalizace s lineární aproximací (bez výpočtu gradientu). .
V závislosti na zvolené metodě jsou podmínky a omezení pro řešení problému nastaveny odlišně:
- objekt třídy
Boundspro metody L-BFGS-B, TNC, SLSQP, trust-constr; - seznam
(min, max)pro stejné metody L-BFGS-B, TNC, SLSQP, trust-constr; - objekt nebo seznam objektů
LinearConstraint,NonlinearConstraintpro metody COBYLA, SLSQP, trust-constr; - slovník nebo seznam slovníků
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}pro metody COBYLA, SLSQP.
Přehled článku:
1) Zvažte použití podmíněného optimalizačního algoritmu v oblasti důvěryhodnosti (method=”trust-constr”) s omezeními určenými jako objekty Bounds, LinearConstraint, NonlinearConstraint ;
2) Zvažte sekvenční programování pomocí metody nejmenších čtverců (metoda = "SLSQP") s omezeními specifikovanými ve formě slovníku {'type', 'fun', 'jac', 'args'};
3) Analyzujte příklad optimalizace vyráběných produktů na příkladu webového studia.
Metoda podmíněné optimalizace="trust-constr"
Implementace metody trust-constr na základě pro problémy s omezeními formy rovnosti a dále pro problémy s omezeními v podobě nerovností. Obě metody jsou implementovány pomocí algoritmů pro nalezení lokálního minima v oblasti spolehlivosti a jsou vhodné pro rozsáhlé problémy.
Matematická formulace problému hledání minima v obecné podobě:
Pro přísná omezení rovnosti je dolní mez nastavena stejně jako horní mez 
.
Pro jednosměrné omezení je nastavena horní nebo dolní mez np.inf s odpovídajícím znakem.
Nechť je třeba najít minimum známé Rosenbrockovy funkce dvou proměnných:
V tomto případě jsou na její definiční doménu nastavena následující omezení:
V našem případě jde o unikátní řešení 
, pro které platí pouze první a čtvrté omezení.
Projdeme si omezení zdola nahoru a podíváme se, jak je můžeme napsat ve scipy.
Omezení 
и 
pojďme jej definovat pomocí objektu Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Omezení 
и 
Zapišme to v lineárním tvaru:
Pojďme definovat tato omezení jako objekt LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])A nakonec nelineární omezení ve formě matice:
Definujeme Jacobovu matici pro toto omezení a lineární kombinaci Hessovy matice s libovolným vektorem 
:
Nyní můžeme definovat nelineární omezení 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á, lze matice zadat také v řídkém tvaru:
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 předmět 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)Při výpočtu hessovské matice 
vyžaduje hodně úsilí, můžete použít třídu . K dispozici jsou následující strategie: BFGS и SR1.
from 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 rozdílů:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Jakobiánskou matici pro omezení lze také vypočítat pomocí konečných rozdílů. V tomto případě však nelze Hessovu matici vypočítat pomocí konečných rozdílů. 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 problému vypadá 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 případě potřeby lze funkci pro výpočet Hessian 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 součin Hessianu 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 druhou derivaci funkce, která se optimalizuje, aproximovat. Pomocí funkce lze například aproximovat Hessian SR1 (kvazi-newtonská aproximace). Gradient lze aproximovat konečnými rozdíly.
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)Metoda podmíněné optimalizace="SLSQP"
Metoda SLSQP je navržena pro řešení problémů minimalizace funkce ve tvaru:
Kde 
и 
— soubory indexů výrazů popisujících omezení ve formě rovnosti nebo nerovností. 
— množiny dolní a horní meze pro definiční obor funkce.
Lineární a nelineární omezení jsou popsána ve formě slovníků s klíči 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])
}Hledání minima se provádí následovně:
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 ]Příklad optimalizace
V souvislosti s přechodem na pátou technologickou strukturu se podívejme na optimalizaci výroby na příkladu webového studia, které nám přináší malý, ale stabilní příjem. Představme si sami sebe jako ředitele kuchyně, která vyrábí tři druhy produktů:
- x0 - prodej vstupních stránek, od 10 tr.
- x1 - firemní weby, od 20 tr.
- x2 - internetové obchody, od 30 tr.
Náš přátelský pracovní tým tvoří čtyři junioři, dva střední a jeden senior. Jejich měsíční fond pracovní doby:
- června:
4 * 150 = 600 чел * час, - středy:
2 * 150 = 300 чел * час, - senor:
150 чел * час.
Nechte prvního dostupného juniora strávit (0, 1, 2) hodin vývojem a nasazením jednoho webu typu (x10, x20, x30), střední - (7, 15, 20), senior - (5, 10, 15 ) hodiny nejlepšího období vašeho života.
Jako každý normální ředitel chceme maximalizovat měsíční zisky. Prvním krokem k úspěchu je zapsat si účelovou funkci value jako výše příjmu z produktů vyrobených za měsíc:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Nejde o chybu, při hledání maxima se cílová funkce minimalizuje s opačným znaménkem.
Dalším krokem je zakázat našim zaměstnancům přepracování a zavést omezení 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ím omezením je, že výstup produktu musí být pouze kladný:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])A nakonec nejrůžovějším předpokladem je, že kvůli nízké ceně a vysoké kvalitě se na nás neustále staví fronta spokojených zákazníků. Můžeme si sami zvolit měsíční objemy výroby na základě vyřešení problému s omezenou optimalizací 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]Volně zaokrouhleme na celá čísla a spočítejme si měsíční zátěž veslařů s optimálním rozložením produktů x = (8, 6, 3) :
- června:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - středy:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Závěr: aby ředitel dostal své zasloužené maximum, je optimální vytvořit 8 vstupních stránek, 6 středně velkých webů a 3 obchody měsíčně. V tomto případě musí senior orat bez zvednutí od stroje, zatížení středů bude přibližně 2/3, juniorů méně než polovina.
Závěr
Článek nastiňuje základní techniky práce s balíčkem scipy.optimize, který se používá k řešení problémů podmíněné minimalizace. Osobně používám scipy čistě pro akademické účely, a proto je uvedený příklad tak komického charakteru.
Spoustu teorie a virtuálních příkladů lze nalézt například v knize I. L. Akulicha „Matematické programování v příkladech a problémech“. Více hardcore aplikace scipy.optimize vytvořit 3D strukturu ze sady obrázků () si můžete prohlédnout v .
Hlavním zdrojem informací je kteří chtějí přispět k překladu této a dalších částí scipy Vítejte v .
Díky za podíl na přípravě publikace.
Zdroj: www.habr.com
