SciPy (izgovara se sai pie) je matematički paket baziran na numpyju koji također uključuje C i Fortran biblioteke. SciPy pretvara vašu interaktivnu Python sesiju u kompletno okruženje za nauku o podacima kao što je MATLAB, IDL, Octave, R ili SciLab.
U ovom članku ćemo se osvrnuti na osnovne tehnike matematičkog programiranja - rješavanje problema uvjetne optimizacije za skalarnu funkciju nekoliko varijabli pomoću paketa scipy.optimize. Algoritmi neograničene optimizacije su već razmatrani u . Detaljnija i ažurnija pomoć o scipy funkcijama uvijek se može dobiti pomoću naredbe help(), Shift+Tab ili u .
Uvod
Uobičajeno sučelje za rješavanje i uvjetnih i neograničenih problema optimizacije u paketu scipy.optimize obezbjeđuje funkcija minimize(). Međutim, poznato je da ne postoji univerzalna metoda za rješavanje svih problema, pa izbor adekvatne metode, kao i uvijek, pada na pleća istraživača.
Odgovarajući algoritam optimizacije je specificiran pomoću argumenta funkcije minimize(..., method="").
Za uslovnu optimizaciju funkcije nekoliko varijabli, dostupne su implementacije sljedećih metoda:
trust-constr— tražiti lokalni minimum u regiji povjerenja. , ;SLSQP— sekvencijalno kvadratno programiranje sa ograničenjima, Njutnova metoda za rešavanje Lagranžovog sistema. .TNC- Skraćeni Newton Ograničen, ograničen broj iteracija, dobar za nelinearne funkcije sa velikim brojem nezavisnih varijabli. .L-BFGS-B— metoda Broyden–Fletcher–Goldfarb–Shanno tima, implementirana uz smanjenu potrošnju memorije zbog djelomičnog učitavanja vektora iz Hessian matrice. , .COBYLA— MARE ograničena optimizacija linearnom aproksimacijom, ograničena optimizacija sa linearnom aproksimacijom (bez izračunavanja gradijenta). .
U zavisnosti od izabrane metode, različito se postavljaju uslovi i ograničenja za rešavanje problema:
- klasni objekat
Boundsza metode L-BFGS-B, TNC, SLSQP, trust-constr; - lista
(min, max)za iste metode L-BFGS-B, TNC, SLSQP, trust-constr; - objekat ili listu objekata
LinearConstraint,NonlinearConstraintza COBYLA, SLSQP, trust-constr metode; - rječnik ili lista rječnika
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}za COBYLA, SLSQP metode.
Pregled članka:
1) Razmotrite upotrebu algoritma uslovne optimizacije u regionu poverenja (metod=”trust-constr”) sa ograničenjima navedenim kao objekti Bounds, LinearConstraint, NonlinearConstraint ;
2) Razmotrite sekvencijalno programiranje koristeći metodu najmanjih kvadrata (metod = "SLSQP") sa 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"
Implementacija metode trust-constr na osnovu za probleme sa ograničenjima oblika jednakosti i dalje za probleme sa ograničenjima u obliku nejednakosti. Obje metode su implementirane pomoću algoritama za pronalaženje lokalnog minimuma u području povjerenja i dobro su prikladne za probleme velikih razmjera.
Matematička formulacija problema nalaženja minimuma u opštem obliku:
Za stroga ograničenja jednakosti, donja granica je postavljena jednakom gornjoj granici 
.
Za jednosmjerno ograničenje postavlja se gornja ili donja granica np.inf sa odgovarajućim znakom.
Neka je potrebno pronaći minimum poznate Rosenbrockove funkcije dvije varijable:
U ovom slučaju, sljedeća ograničenja su postavljena na njegovu domenu definicije:
U našem slučaju postoji jedinstveno rješenje u tom trenutku 
, 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 
и 
hajde da ga definišemo koristeći objekat 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 konačno, nelinearno ograničenje u matričnom obliku:
Definiramo Jacobian matricu za ovo ograničenje i linearnu kombinaciju Hessian matrice sa 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 objekat 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)Prilikom izračunavanja Hessove 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 korištenjem konačnih razlika:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Jacobian matrica za ograničenja se također može izračunati korištenjem konačnih razlika. Međutim, u ovom slučaju Hessian matrica se ne može izračunati korištenjem konačnih razlika. Hessian mora biti definiran kao funkcija ili korištenjem 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 Hessiana može se definirati korištenjem 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 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, prvi i drugi izvod funkcije koja se optimizira mogu se aproksimirati. Na primjer, Hessian se može aproksimirati pomoću funkcije SR1 (kvazi-njutnova 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 uslovne optimizacije="SLSQP"
SLSQP metoda je dizajnirana da riješi probleme minimiziranja funkcije u obliku:
Gde 
и 
— skupovi indeksa izraza koji opisuju ograničenja u obliku jednakosti ili nejednakosti. 
— skupovi donjih i gornjih granica za domenu definicije funkcije.
Linearna i nelinearna ograničenja su opisana 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])
}Potraga za minimumom se vrši 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 sa prelaskom na petu tehnološku strukturu, pogledajmo optimizaciju proizvodnje na primjeru web studija, koji nam donosi mali, ali stabilan prihod. Zamislimo sebe kao direktora kuhinje koja proizvodi tri vrste proizvoda:
- x0 - prodaja odredišnih stranica, od 10 tr.
- x1 - korporativne web stranice, od 20 tr.
- x2 - online prodavnice, od 30 tr.
Naš prijateljski radni tim čine četiri juniora, dva srednja i jedan senior. Njihov mjesečni fond radnog vremena:
- juni:
4 * 150 = 600 чел * час, - sredine:
2 * 150 = 300 чел * час, - senor:
150 чел * час.
Neka prvi raspoloživi junior potroši (0, 1, 2) sati na razvoj i implementaciju jednog sajta tipa (x10, x20, x30), srednji - (7, 15, 20), stariji - (5, 10, 15 ) sati najboljeg vremena u vašem životu.
Kao i svaki normalan direktor, želimo maksimizirati mjesečni profit. Prvi korak do uspjeha je zapisivanje ciljne funkcije value kao iznos prihoda od proizvedenih proizvoda mjesečno:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Ovo nije greška kada se traži maksimum, ciljna funkcija je minimizirana sa suprotnim predznakom.
Sljedeći korak je zabraniti našim zaposlenima preopterećenje i uvesti ograničenja radnog vremena:
Šta 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žočnija pretpostavka je da se zbog niske cijene i visokog kvaliteta za nas stalno kruži red zadovoljnih kupaca. Mjesečne količine proizvodnje možemo odabrati sami, na osnovu rješavanja problema ograničene optimizacije sa 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 uz optimalnu distribuciju proizvoda x = (8, 6, 3) :
- juni:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - sredine:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Zaključak: da bi direktor dobio svoj zasluženi maksimum, optimalno je mjesečno kreirati 8 landing stranica, 6 stranica srednje veličine i 3 trgovine. Istovremeno, seniori moraju orati bez gledanja sa mašine, opterećenje srednjih je otprilike 2/3, a juniori manje od polovine.
zaključak
U članku su navedene osnovne tehnike rada s paketom scipy.optimize, koji se koristi za rješavanje problema uvjetne minimizacije. Lično koristim scipy čisto u akademske svrhe, zbog čega je navedeni primjer tako komične prirode.
Mnogo teorije i virtuelnih primjera može se naći, na primjer, u knjizi I.L. Akulicha “Matematičko programiranje u primjerima i problemima”. Više hardcore aplikacija scipy.optimize da napravite 3D strukturu od skupa slika () možete pogledati u .
Glavni izvor informacija je oni koji žele da doprinesu prevođenju ovog i drugih delova scipy Dobrodošli .
Spasibo za učešće u pripremi publikacije.
izvor: www.habr.com
