SciPy (izgovara se sai pie) je matematički paket temeljen na numpyju koji također uključuje C i Fortran biblioteke. SciPy pretvara vašu interaktivnu Python sesiju u potpuno okruženje znanosti o podacima kao što su MATLAB, IDL, Octave, R ili SciLab.
U ovom ćemo članku pogledati osnovne tehnike matematičkog programiranja – rješavanje problema uvjetne optimizacije za skalarnu funkciju nekoliko varijabli pomoću paketa scipy.optimize. O algoritmima neograničene optimizacije već je bilo riječi u . Detaljnija i ažurirana pomoć o scipy funkcijama uvijek se može dobiti pomoću naredbe help(), Shift+Tab ili u .
Uvod
Zajedničko sučelje za rješavanje problema uvjetne i neograničene optimizacije u paketu scipy.optimize pruža funkcija minimize(). No, poznato je da ne postoji univerzalna metoda za rješavanje svih problema, pa je izbor adekvatne metode, kao i uvijek, na plećima istraživača.
Odgovarajući algoritam optimizacije specificiran je pomoću argumenta funkcije minimize(..., method="").
Za uvjetnu optimizaciju funkcije nekoliko varijabli dostupne su implementacije sljedećih metoda:
trust-constr— traženje lokalnog minimuma u području povjerenja. , ;SLSQP— sekvencijalno kvadratno programiranje s ograničenjima, Newtonova metoda za rješavanje Lagrangeova sustava. .TNC- Krnji Newton ograničen, ograničen broj ponavljanja, dobar za nelinearne funkcije s velikim brojem nezavisnih varijabli. .L-BFGS-B— metoda tima Broyden–Fletcher–Goldfarb–Shanno, implementirana sa smanjenom potrošnjom memorije zbog djelomičnog učitavanja vektora iz Hessian matrice. , .COBYLA— MARE ograničena optimizacija linearnom aproksimacijom, ograničena optimizacija s linearnom aproksimacijom (bez proračuna gradijenta). .
Ovisno o odabranoj metodi različito se postavljaju uvjeti i ograničenja za rješavanje problema:
- objekt klase
Boundsza metode L-BFGS-B, TNC, SLSQP, trust-constr; - popis
(min, max)za iste metode L-BFGS-B, TNC, SLSQP, trust-constr; - objekt ili popis objekata
LinearConstraint,NonlinearConstraintza COBYLA, SLSQP, trust-constr metode; - rječnik ili popis rječnika
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}za metode COBYLA, SLSQP.
Pregled članka:
1) Razmotrite upotrebu algoritma uvjetne optimizacije u regiji povjerenja (method=”trust-constr”) s ograničenjima navedenim kao objekti Bounds, LinearConstraint, NonlinearConstraint ;
2) Razmotrite sekvencijalno programiranje korištenjem metode najmanjih kvadrata (metoda = "SLSQP") s ograničenjima navedenim u obliku rječnika {'type', 'fun', 'jac', 'args'};
3) Analizirati primjer optimizacije proizvedenih proizvoda na primjeru web studija.
Metoda uvjetne optimizacije="trust-constr"
Provedba metode trust-constr na temelju za probleme s ograničenjima oblika jednakosti i na za probleme s ograničenjima u obliku nejednakosti. Obje su metode implementirane algoritmima za pronalaženje lokalnog minimuma u području pouzdanosti i dobro su prikladne za probleme velikih razmjera.
Matematička formulacija problema nalaženja minimuma u općem obliku:
Za stroga ograničenja jednakosti, donja granica je jednaka gornjoj granici .
Za jednosmjerno ograničenje postavlja se gornja ili donja granica np.inf s pripadajućim znakom.
Neka je potrebno pronaći minimum poznate Rosenbrockove funkcije dviju varijabli:
U ovom slučaju, sljedeća su ograničenja postavljena na njegovu domenu definicije:
U našem slučaju, u točki postoji jedinstveno rješenje , za koje vrijede samo prvo i četvrto ograničenje.
Prođimo kroz ograničenja odozdo prema gore i pogledajmo kako ih možemo napisati u scipyju.
Ograničenja и definirajmo ga pomoću objekta Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Ograničenja и Zapišimo to u linearnom obliku:
Definirajmo ova ograničenja kao objekt LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])I na kraju nelinearno ograničenje u matričnom obliku:
Definiramo Jacobijevu matricu za ovo ograničenje i linearnu kombinaciju Hessianove matrice s proizvoljnim vektorom :
Sada možemo definirati nelinearno ograničenje 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 također mogu specificirati 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:
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 izračunavanju Hessianove matrice zahtijeva puno truda, možete koristiti klasu . Dostupne su sljedeće strategije: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian 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')Jacobianova matrica za ograničenja također se može izračunati korištenjem konačnih razlika. Međutim, u ovom slučaju Hessova matrica ne može se izračunati pomoću konačnih razlika. Hessian se mora definirati kao funkcija ili pomoću klase HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Rješenje problema optimizacije izgleda ovako:
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]Ako je potrebno, funkcija za izračunavanje Hessana 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 umnožak Hessiana 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 funkcije koja se optimizira može se aproksimirati. Na primjer, Hessian se može aproksimirati pomoću funkcije SR1 (kvazi-Newtonova aproksimacija). Gradijent se može aproksimirati konačnim razlikama.
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 uvjetne optimizacije="SLSQP"
Metoda SLSQP dizajnirana je za rješavanje problema minimiziranja funkcije u obliku:
gdje и — skupovi indeksa izraza koji opisuju ograničenja u obliku jednakosti ili nejednakosti. — skupovi donjih i gornjih granica za domenu definiranja funkcije.
Linearna i nelinearna ograničenja opisana su u obliku rječnika s ključevima 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])
}Traženje minimuma provodi se na sljedeći način:
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 ]Primjer optimizacije
U vezi s prelaskom na petu tehnološku strukturu, pogledajmo optimizaciju proizvodnje na primjeru web studija koji nam donosi mali, ali stabilan prihod. Zamislimo sebe kao direktora brodske kuhinje koja proizvodi tri vrste proizvoda:
- x0 - prodaja odredišnih stranica, od 10 tr.
- x1 - korporativne web stranice, od 20 tr.
- x2 - online trgovine, od 30 tr.
Naš prijateljski radni tim sastoji se od četiri juniora, dva srednje i jednog seniora. Njihov mjesečni fond radnog vremena:
- lipnja:
4 * 150 = 600 чел * час, - sredina:
2 * 150 = 300 чел * час, - senjor:
150 чел * час.
Neka prvi raspoloživi mlađi provede (0, 1, 2) sati na razvoju i postavljanju jednog mjesta tipa (x10, x20, x30), srednji - (7, 15, 20), viši - (5, 10, 15 ) sati najboljeg provoda u vašem životu.
Kao svaki normalan direktor, želimo maksimizirati mjesečnu dobit. Prvi korak do uspjeha je zapisati funkciju cilja value kao iznos prihoda od proizvoda proizvedenih mjesečno:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Ovo nije pogreška, pri traženju maksimuma funkcija cilja se minimizira s suprotnim predznakom.
Sljedeći korak je zabrana prekomjernog rada našim zaposlenicima i uvođenje ograničenja 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 je da izlaz proizvoda mora biti samo pozitivan:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])I na kraju, najružičastija pretpostavka je da se zbog niske cijene i visoke kvalitete za nama stalno ređa red zadovoljnih kupaca. Možemo sami odabrati mjesečne količine proizvodnje na temelju rješavanja problema ograničene 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]Zaokružimo slobodno na cijele brojeve i izračunajmo mjesečno opterećenje veslača s optimalnom raspodjelom proizvoda x = (8, 6, 3) :
- lipnja:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - sredina:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senjor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Zaključak: kako bi direktor dobio svoj zasluženi maksimum, optimalno je izraditi 8 odredišnih stranica, 6 stranica srednje veličine i 3 trgovine mjesečno. U ovom slučaju, senior mora orati bez podizanja pogleda sa stroja, opterećenje srednjih bit će otprilike 2/3, juniori manje od polovice.
Zaključak
U članku se opisuju osnovne tehnike rada s paketom scipy.optimize, koristi se za rješavanje problema uvjetne minimizacije. Osobno koristim scipy čisto u akademske svrhe, zbog čega je navedeni primjer tako komičan.
Puno teorije i virtualnih primjera može se pronaći, na primjer, u knjizi I.L. Akulich "Matematičko programiranje u primjerima i problemima." Više hardcore aplikacija scipy.optimize izgraditi 3D strukturu iz skupa slika () možete pogledati u .
Glavni izvor informacija je onima koji žele pridonijeti prijevodu ovog i drugih odjeljaka scipy Dobrodošli u .
Hvala za sudjelovanje u pripremi publikacije.
Izvor: www.habr.com