SciPy (izgovara se kao saj paj) je matematički paket koji se temelji na numpy-u i uključuje biblioteke napisane na C-u i Fortranu. Uz SciPy, interaktivna Python sesija postaje potpuna okolina za obradu podataka, poput MATLAB-a, IDL-a, Octave-a, R-a ili SciLab-a.
U ovom članku razmotrit ćemo osnovne tehnike matematičkog programiranja — rješavanje problema uvjetne optimizacije za skalarne funkcije više varijabli pomoću paketa scipy.optimize. Algoritmi bezuvjetne optimizacije već su obrađeni u . Detaljna i aktualna pomoć o funkcijama scipy uvijek se može dobiti pomoću naredbe help(), Shift+Tab ili u .
Uvod
Opći интерфейс za rješavanje problema i uvjetne i bezuvjetne optimizacije u paketu scipy.optimize pruža funkcija minimize(). Međutim, poznato je da ne postoji univerzalni način za rješavanje svih problema, stoga izbor adekvatne metode kao i uvijek pada na teret istraživača.
Pogodan algoritam optimizacije određuje se pomoću argumenta funkcije minimize(..., method="").
Za uvjetnu optimizaciju funkcija više varijabli dostupne su implementacije sljedećih metoda:
trust-constr— traženje lokalnog minimuma u povjerenstvenom području. , ;SLSQP— sekvencijsko kvadratno programiranje s ograničenjima, Newtonova metoda rješavanja sustava Lagrangea. .TNC— Truncated Newton Constrained, ograničen broj iteracija, dobar za nelinearne funkcije s velikim brojem neovisnih varijabla. .L-BFGS-B— metoda četvorke Broyden-Fletcher-Goldfarb-Shanno, implementirana s smanjenom potrošnjom memorije korištenjem parcijalnog učitavanja vektora iz Hesseove matrice. , .COBYLA— COBYLA Constrained Optimization By Linear Approximation, ograničena optimizacija s linearnom aproksimacijom (bez izračunavanja gradijenta). .
Ovisno o odabranoj metodi, različito se postavljaju uvjeti i ograničenja za rješavanje problema:
- objektom klase
Boundsza metode L-BFGS-B, TNC, SLSQP, trust-constr; - popisom
(min, max)za iste metode L-BFGS-B, TNC, SLSQP, trust-constr; - objektom ili popisom objekata
LinearConstraint,NonlinearConstraintza metode COBYLA, SLSQP, trust-constr; - rječnikom ili popisom rječnika
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}za metode COBYLA, SLSQP.
Plan članka:
1) Razmotrite primjenu algoritma uvjetne optimizacije u povjerenoj oblasti (method=»trust-constr») s ograničenjima koja su zadana u obliku objekata Bounds, LinearConstraint, NonlinearConstraint ;
2) Razmotrite sekvencijsko programiranje metodom najmanjih kvadrata (method=»SLSQP») s ograničenjima koja su zadana u obliku rječnika {'type', 'fun', 'jac', 'args'};
3) Razmotrite primjer optimizacije proizvodnje na primjeru web-studija.
Uvjetna optimizacija method=»trust-constr»
Implementacija metode trust-constr temelji se na za probleme s ograničenjima u obliku jednadžbe i na za probleme s ograničenjima u obliku nejednakosti. Ove metode su implementirane algoritmima pretraživanja lokalnog minimuma u povjerenoj oblasti i dobro su prilagođene za velike probleme.
Matematička formulacija problema pronalaženja minimuma u općem obliku:
Za ograničenja stroge jednakosti, donja granica se postavlja jednaka gornjoj 
.
Za jednostrano ograničenje, gornja ili donja granica se postavlja np.inf s odgovarajućim znakom.
Neka je potrebno pronaći minimum poznate funkcije Rosenbrock od dva varijable:
Pri tome su zadana sljedeća ograničenja na njezinu oblast definicije:
U našem slučaju postoji jedinstveno rješenje u točki 
, za koju su valjana samo prva i četvrta ograničenja.
Proći ćemo kroz ograničenja od dna prema vrhu i razmotriti kako ih možemo zapisati u scipy.
Ograničenja 
i 
određujemo pomoću objekta Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Ograničenja 
i 
zapisujemo u linearnoj formi:
Odredit ćemo ta ograničenja u obliku objekta LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])I konačno, nelinearno ograničenje u matričnom obliku:
Odredimo Jacobi matricu za ovo ograničenje i linearne kombinacije Hessove matrice s proizvoljnim vektorom 
:
Sada nelinearno ograničenje možemo odrediti kao 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)Ako je veličina velika, matrice se mogu zadati i u rijetkom obliku:
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)ili kao objekt LinearOperator:
iz 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)Kada izračunavanje Hesijanove matrice 
zahtijeva velike troškove, može se upotrijebiti klasa . Dostupne su sljedeće strategije: BFGS i SR1.
iz scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hesijan se također može izračunati pomoću konačnih razlika:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Jacobi matrica za ograničenja također se može izračunati pomoću konačnih razlika. Međutim, u ovom slučaju Hesijan se ne može izračunati pomoću konačnih razlika. Hesijan mora biti definiran kao funkcija ili putem klase HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Rješenje optimizacijskog problema izgleda ovako:
iz scipy.optimize import minimize
iz 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` uvjet za prekid je zadovoljen.
Broj različitih iteracija: 12, evaluacije funkcije: 8, CG iteracije: 7, optimalnost: 2.99e-09, kršenje ograničenja: 1.11e-16, vrijeme izvođenja: 0.033 s.
[0.41494531 0.17010937]Ako je potrebno, funkciju za izračunavanje Hesijana može se definirati pomoću klase 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)ili proizvod Hesijana i proizvoljnog vektora kroz parametar 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)Alternativno, prva i druga derivacija optimizirane funkcije mogu se približno izračunati. Na primjer, Hesijan može biti aproksimiran pomoću funkcije SR1 (kvazi-Njutonske aproksimacije). Gradient se može aproksimirati konačnim razlikama.
iz 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)Kondicionalna optimizacija method=»SLSQP»
Metoda SLSQP namijenjena je rješavanju problema minimizacije funkcije u obliku:
Gdje 
i 
— skupovi indeksa izraza koji opisuju ograničenja u obliku jednakosti ili nejednakosti. 
— skupovi donjih i gornjih granica za domen funkcije.
Linearna i nelinearna ograničenja opisuju se u obliku rečnika sa ključevima type, fun i 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])
}Postupak traženja minimuma je sledeći:
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)Optimizacija je uspešno završena. (Izlazni mod 0)
Trenutna vrednost funkcije: 0.34271757499419825
Iteracije: 4
Evaluacije funkcije: 5
Evaluacije gradijenta: 4
[0.41494475 0.1701105 ]Primer optimizacije
U svetlu prelaska na peti tehnološki ciklus, razmotrićemo optimizaciju proizvodnje na primeru veb-studija koja nam donosi mali, ali stabilan prihod. Zamislimo se kao direktor galerije koja proizvodi tri vrste proizvoda:
- x0 — prodajni landing strani, od 10.000 kn.
- x1 — korporativne veb stranice, od 20.000 kn.
- x2 — internet prodavnice, od 30.000 kn.
Naš tim se sastoji od četiri juniora, dva midla i jednog seniora. Fond njihovog radnog vremena za mesec:
- juniori:
4 * 150 = 600 čovek-sati, - midlovi:
2 * 150 = 300 čovek-sati, - senior:
150 čovek-sati.
Neka za razvoj i implementaciju jednog tipa veb sajta (x0, x1, x2) prvi deo juniora potroši (10, 20, 30) sati, midl — (7, 15, 20), senior — (5, 10, 15) sati svog najboljeg vremena.
Kao što svakom normalnom direktoru priliči, želimo maksimizovati mesečni profit. Prvi korak ka uspehu — zapisujemo ciljnu funkciju value kao zbir prihoda od proizvedene robe tokom meseca:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]To nije greška, prilikom traženja maksimuma, ciljana funkcija se minimizira sa obrnutim znakom.
Sledeći korak — zabranjujemo prekovremeni rad svojim zaposlenima i uvodimo ograničenja na fond radnog vremena:
Što je ekvivalentno:
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]])
}Formalno ograničenje - proizvodnja mora biti isključivo pozitivna:
bnds = Boundaries ([0, 0, 0], [np.inf, np.inf, np.inf])I konačno, najoptimističnija pretpostavka - zbog niske cijene i visoke kvalitete, konstantno se formira red zadovoljnih klijenata. Možemo sami odabrati mjesečne obujme proizvodnje, temeljeći se na rješavanju zadatka uvjetne optimizacije 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]Neslužbeno ćemo zaokružiti na cijele brojeve i izračunati mjesečnu opterećenost veslača pri optimalnoj raspodjeli proizvoda x = (8, 6, 3) :
- juniori:
8 * 10 + 6 * 20 + 3 * 30 = 290 ljudi * sati; - midlovi:
8 * 7 + 6 * 15 + 3 * 20 = 206 ljudi * sati; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 ljudi * sati.
Zaključak: kako bi direktor ostvario svoj zasluženi maksimum, optimalno je proizvoditi 8 odredišnih stranica, 6 srednje velikih web stranica i 3 trgovine mjesečno. Senior bi trebao raditi neprekidno, opterećenje mid-level bi iznosilo otprilike 2/3, dok bi juniori imali manje od polovice.
Zaključak
U članku su izložene osnovne tehnike rada s paketom scipy.optimize, koje se koriste za rješavanje zadataka uvjetne minimalizacije. Osobno koristim scipy isključivo u akademske svrhe, pa primjer koji slijedi ima za cilj duhovitost.
Mnogo teorije i raritetnih primjera može se pronaći, na primjer, u knjizi I.L. Akuliča "Matematičko programiranje u primjerima i zadacima". Teži korištenje scipy.optimize za izradu 3D strukture prema skupu slika () može se vidjeti u .
Glavni izvor informacija - , oni koji žele doprinijeti prevodu ovog i drugih odjeljaka scipy dobrodošli su na .
Hvala za sudjelovanje u pripremi publikacije.
Izvor: habr.com
