SciPy (izgovorjeno sai pie) je matematični paket, ki temelji na numpyju in vključuje tudi knjižnici C in Fortran. SciPy vašo interaktivno sejo Python spremeni v popolno podatkovno znanstveno okolje, kot so MATLAB, IDL, Octave, R ali SciLab.
V tem članku si bomo ogledali osnovne tehnike matematičnega programiranja – reševanje problemov pogojne optimizacije za skalarno funkcijo več spremenljivk z uporabo paketa scipy.optimize. Algoritmi neomejene optimizacije so bili že obravnavani v . Podrobnejšo in posodobljeno pomoč pri funkcijah scipy lahko vedno dobite z uporabo ukaza help(), Shift+Tab ali v .
Predstavitev
Skupni vmesnik za reševanje problemov pogojne in neomejene optimizacije v paketu scipy.optimize zagotavlja funkcija minimize(). Znano pa je, da univerzalne metode za reševanje vseh problemov ni, zato je izbira ustrezne metode kot vedno na ramenih raziskovalca.
Ustrezen optimizacijski algoritem je določen z uporabo argumenta funkcije minimize(..., method="").
Za pogojno optimizacijo funkcije več spremenljivk so na voljo izvedbe naslednjih metod:
trust-constr— iskanje lokalnega minimuma v območju zaupanja. , ;SLSQP— zaporedno kvadratno programiranje z omejitvami, Newtonova metoda za reševanje Lagrangeovega sistema. .TNC- Prisekano Newtonovo omejeno, omejeno število ponovitev, dobro za nelinearne funkcije z velikim številom neodvisnih spremenljivk. .L-BFGS-B— metoda ekipe Broyden–Fletcher–Goldfarb–Shanno, implementirana z zmanjšano porabo pomnilnika zaradi delnega nalaganja vektorjev iz Hessove matrike. , .COBYLA— Omejena optimizacija MARE z linearno aproksimacijo, omejena optimizacija z linearno aproksimacijo (brez izračuna gradienta). .
Glede na izbrano metodo so različno postavljeni pogoji in omejitve za rešitev problema:
- predmet razreda
Boundsza metode L-BFGS-B, TNC, SLSQP, trust-constr; - seznam
(min, max)za iste metode L-BFGS-B, TNC, SLSQP, trust-constr; - predmet ali seznam predmetov
LinearConstraint,NonlinearConstraintza metode COBYLA, SLSQP, trust-constr; - slovar ali seznam slovarjev
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}za metode COBYLA, SLSQP.
Oris članka:
1) Razmislite o uporabi algoritma pogojne optimizacije v območju zaupanja (method=”trust-constr”) z omejitvami, določenimi kot objekti Bounds, LinearConstraint, NonlinearConstraint ;
2) Razmislite o zaporednem programiranju z uporabo metode najmanjših kvadratov (metoda = "SLSQP") z omejitvami, določenimi v obliki slovarja {'type', 'fun', 'jac', 'args'};
3) Analiziraj primer optimizacije izdelanih izdelkov na primeru spletnega studia.
Metoda pogojne optimizacije="trust-constr"
Izvedba metode trust-constr temelji na za težave z omejitvami oblike enakosti in na za probleme z omejitvami v obliki neenačb. Obe metodi izvajata algoritma za iskanje lokalnega minimuma v območju zaupanja in sta zelo primerni za probleme velikega obsega.
Matematična formulacija problema iskanja minimuma v splošni obliki:
Za stroge omejitve enakosti je spodnja meja enaka zgornji meji .
Za enosmerno omejitev je nastavljena zgornja ali spodnja meja np.inf z ustreznim znakom.
Naj bo treba najti minimum znane Rosenbrockove funkcije dveh spremenljivk:
V tem primeru so na domeni definicije določene naslednje omejitve:
V našem primeru gre za edinstveno rešitev , za katerega veljata samo prva in četrta omejitev.
Oglejmo si omejitve od spodaj navzgor in poglejmo, kako jih lahko zapišemo v scipy.
Omejitve и definirajmo ga z uporabo objekta Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Omejitve и Zapišimo ga v linearni obliki:
Definirajmo te omejitve kot objekt LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])In končno nelinearna omejitev v matrični obliki:
Definiramo Jacobijevo matriko za to omejitev in linearno kombinacijo Hessove matrike s poljubnim vektorjem :
Zdaj lahko definiramo nelinearno omejitev kot 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)Če je velikost velika, lahko matrike podate tudi v redki obliki:
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)ali kot predmet 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čunu Hessove matrike zahteva veliko truda, lahko uporabite razred . Na voljo so naslednje strategije: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian se lahko izračuna tudi z uporabo končnih razlik:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Jacobijevo matriko za omejitve je mogoče izračunati tudi z uporabo končnih razlik. Vendar v tem primeru Hessove matrike ni mogoče izračunati z uporabo končnih razlik. Hessian mora biti definiran kot funkcija ali z uporabo razreda HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Rešitev optimizacijskega problema izgleda takole:
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]Če je potrebno, lahko funkcijo za izračun Hessian definirate z uporabo razreda 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)ali zmnožek hessenovega in poljubnega vektorja skozi parameter 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)Druga možnost je, da sta prvi in drugi odvod funkcije, ki se optimizira, aproksimirana. Hessian se lahko na primer približa s funkcijo SR1 (kvazi-Newtonov približek). Gradient je mogoče aproksimirati s končnimi razlikami.
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 pogojne optimizacije="SLSQP"
Metoda SLSQP je zasnovana za reševanje problemov minimiziranja funkcije v obliki:
kjer je и — nizi indeksov izrazov, ki opisujejo omejitve v obliki enakosti ali neenakosti. — nizi spodnjih in zgornjih meja za področje definicije funkcije.
Linearne in nelinearne omejitve so opisane v obliki slovarjev s ključ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])
}Iskanje minimuma poteka na naslednji 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 ]Primer optimizacije
V povezavi s prehodom na peto tehnološko strukturo si poglejmo optimizacijo proizvodnje na primeru spletnega studia, ki nam prinaša majhen, a stabilen dohodek. Predstavljajmo si sebe kot direktorja kuhinje, ki proizvaja tri vrste izdelkov:
- x0 - prodaja ciljnih strani, od 10 tr.
- x1 - korporativne spletne strani, od 20 tr.
- x2 - spletne trgovine, od 30 tr.
Našo prijazno delovno ekipo sestavljajo štirje mlajši, dva srednja in en starejši. Njihov mesečni fond delovnega časa:
- junij:
4 * 150 = 600 чел * час, - sredina:
2 * 150 = 300 чел * час, - senor:
150 чел * час.
Naj prvi razpoložljivi mlajši porabi (0, 1, 2) ur za razvoj in uvajanje enega mesta tipa (x10, x20, x30), srednji - (7, 15, 20), višji - (5, 10, 15 ) ure najboljšega časa v vašem življenju.
Kot vsak običajen direktor, želimo maksimirati mesečni dobiček. Prvi korak do uspeha je zapis ciljne funkcije value kot znesek dohodka od proizvedenih izdelkov na mesec:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]To ni napaka, pri iskanju maksimuma se ciljna funkcija minimizira z nasprotnim predznakom.
Naslednji korak je prepoved prekomernega dela našim zaposlenim in uvedba omejitev delovnega časa:
Kaj je enakovredno:
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]])
}Formalna omejitev je, da mora biti rezultat izdelka le pozitiven:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])In končno, najbolj rožnata predpostavka je, da se zaradi nizke cene in visoke kakovosti za nas nenehno vije vrsta zadovoljnih strank. Mesečne količine proizvodnje lahko izbiramo sami, na podlagi reševanja problema omejene optimizacije z 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]Ohlapno zaokrožimo na cela števila in izračunamo mesečno obremenitev veslačev z optimalno porazdelitvijo zmnožkov x = (8, 6, 3) :
- junij:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - sredina:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Zaključek: da bi direktor prejel svoj zasluženi maksimum, je optimalno ustvariti 8 ciljnih strani, 6 srednje velikih mest in 3 trgovine na mesec. V tem primeru mora starejši orati, ne da bi dvignil pogled s stroja, obremenitev srednjih bo približno 2/3, mlajših manj kot polovica.
Zaključek
Članek opisuje osnovne tehnike dela s paketom scipy.optimize, ki se uporablja za reševanje problemov pogojne minimizacije. Osebno uporabljam scipy zgolj v akademske namene, zato je navedeni primer tako komičen.
Veliko teorije in virtualnih primerov lahko najdete na primer v knjigi I. L. Akulich "Matematično programiranje v primerih in problemih." Bolj zahtevna aplikacija scipy.optimize zgraditi 3D strukturo iz nabora slik () si lahko ogledate v .
Glavni vir informacij je tiste, ki želijo prispevati k prevodu tega in drugih razdelkov scipy Dobrodošli v .
Hvala za sodelovanje pri pripravi publikacije.
Vir: www.habr.com