SciPy (pronunciatu sai pie) hè un pacchettu matematicu basatu in numpy chì include ancu biblioteche C è Fortran. SciPy trasforma a vostra sessione Python interattiva in un ambiente cumpletu di scienza di dati cum'è MATLAB, IDL, Octave, R, o SciLab.
In questu articulu, fighjemu e tecniche basi di prugrammazione matematica - risolve i prublemi di ottimisazione cundizionale per una funzione scalare di parechje variàbili cù u pacchettu scipy.optimize. L'algoritmi di ottimisazione senza restrizioni sò digià discututi in . Aiutà più dettagliata è aghjurnata nantu à e funzioni scipy pò sempre esse ottenuta usendu u cumandimu help(), Shift+Tab o in .
Introduzione
Una interfaccia cumuna per risolve i prublemi di ottimisazione cundizionale è senza restrizioni in u pacchettu scipy.optimize hè furnita da a funzione minimize(). In ogni casu, hè cunnisciutu chì ùn ci hè micca un metudu universale per risolve tutti i prublemi, cusì l'scelta di un metudu adattatu, cum'è sempre, casca nantu à e spalle di l'investigatore.
L'algoritmu di ottimisazione adattatu hè specificatu cù l'argumentu di a funzione minimize(..., method="").
Per l'ottimisazione condicionale di una funzione di parechje variàbili, implementazioni di i metudi seguenti sò dispunibili:
trust-constr- cercate un minimu locale in a regione di cunfidenza. , ;SLSQP- prugrammazione quadratica sequenziale cù restrizioni, metudu newtonianu per risolve u sistema di Lagrange. .TNC- Truncated Newton Constrained, numeru limitatu di iterazioni, bonu per funzioni non lineari cù un gran numaru di variabili indipendenti. .L-BFGS-B- un metudu da a squadra Broyden–Fletcher–Goldfarb–Shanno, implementatu cù u cunsumu di memoria ridutta per via di a carica parziale di vettori da a matrice Hessian. , .COBYLA— MARE Constrained Optimization By Linear Approssimazione, ottimisazione restritta cù approssimazioni lineari (senza calculu di gradiente). .
Sicondu u metudu sceltu, e cundizioni è e restrizioni per risolve u prublema sò stabilite in modu diversu:
- oggettu di classa
Boundsper i metudi L-BFGS-B, TNC, SLSQP, trust-constr; - a lista
(min, max)per i stessi metudi L-BFGS-B, TNC, SLSQP, trust-constr; - un ughjettu o una lista d'oggetti
LinearConstraint,NonlinearConstraintper COBYLA, SLSQP, metudi trust-constr; - dizziunariu o lista di dizziunari
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}per i metudi COBYLA, SLSQP.
Schema di l'articulu:
1) Cunsiderate l'usu di un algoritmu di ottimisazione cundizionale in a regione di fiducia (metudu = "trust-constr") cù limitazioni specificate cum'è oggetti Bounds, LinearConstraint, NonlinearConstraint ;
2) Cunsiderate a prugrammazione sequenziale utilizendu u metudu di i minimi quadrati (metudu = "SLSQP") cù restrizioni specificate in forma di dizziunariu {'type', 'fun', 'jac', 'args'};
3) Analizà un esempiu di ottimisazione di i prudutti fabbricati cù l'esempiu di un studiu web.
Metudu di ottimisazione cundizionale = "trust-constr"
Implementazione di u metudu trust-constr bastu nantu à per prublemi cù limitazioni di a forma di ugualità è nantu per i prublemi cù limitazioni in forma di inuguaglianze. I dui metudi sò implementati da algoritmi per truvà un minimu locale in a regione di cunfidenza è sò adattati per prublemi à grande scala.
Formulazione matematica di u prublema di truvà un minimu in forma generale:
Per restrizioni strette di ugualità, u limite inferiore hè stabilitu uguale à u limite superiore .
Per una limitazione unidirezionale, u limitu superiore o inferiore hè stabilitu np.inf cù u signu currispundente.
Chì hè necessariu di truvà u minimu di una funzione Rosenbrock cunnisciuta di duie variàbili:
In questu casu, e seguenti restrizioni sò stabilite in u so duminiu di definizione:
In u nostru casu, ci hè una suluzione unica à u puntu , per quale solu a prima è a quarta restrizioni sò valide.
Andemu à traversu e restrizioni da u fondu à a cima è fighjemu cumu pudemu scrive in scipy.
Restrictions и definiscemu cù l'ughjettu Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restrictions и Scrivemu in forma lineare:
Definimu sti limitazioni cum'è un oggettu LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])È infine a limitazione non lineare in forma di matrice:
Definimu a matrice giacobiana per questa limitazione è una cumminazione lineale di a matrice Hessian cun un vettore arbitrariu :
Avà pudemu definisce una limitazione non lineale cum'è un oggettu 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)Se a dimensione hè grande, e matrici ponu ancu esse specificate in forma sparsa:
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)o cum'è un oggettu 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)Quandu calcula a matrice Hessian richiede assai sforzu, pudete aduprà una classa . E seguenti strategie sò dispunibili: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())U Hessian pò ancu esse calculatu cù differenzi finiti:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')A matrice giacobiana per e restrizioni pò ancu esse calculata utilizendu differenze finite. Tuttavia, in stu casu, a matrice Hessian ùn pò esse calculata cù diffirenzi finiti. U Hessian deve esse definitu cum'è una funzione o utilizendu a classa HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())A suluzione à u prublema di ottimisazione s'assumiglia cusì:
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]Se necessariu, a funzione per u calculu di Hessian pò esse definita cù a classa 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)o u pruduttu di u Hessian è un vettore arbitrariu attraversu u paràmetru 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)In alternativa, a prima è a seconda derivativa di a funzione chì hè ottimizzata pò esse apprussimata. Per esempiu, u Hessian pò esse apprussimata cù a funzione SR1 (approssimazione quasi newtoniana). U gradiente pò esse apprussimatu da differenze finite.
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)Metudu di ottimisazione cundizionale = "SLSQP"
U metudu SLSQP hè designatu per risolve i prublemi di minimizzà una funzione in a forma:
Induve и - setti di indici di espressioni chì descrizanu restrizioni in forma di uguaglianze o inuguaglianze. - insemi di limiti inferiori è superiori per u duminiu di definizione di a funzione.
E restrizioni lineari è non lineari sò descritti in forma di dizionari cù chjave 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])
}A ricerca di u minimu hè realizatu cusì:
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 ]Esempiu di ottimisazione
In cunnessione cù a transizione à a quinta struttura tecnologica, fighjemu l'ottimisazione di a produzzione cù l'esempiu di un studiu web, chì ci porta un ingressu chjucu, ma stabile. Imaginemu cum'è u direttore di una galera chì produce trè tippi di prudutti:
- x0 - vende pagine di destinazione, da 10 tr.
- x1 - siti web corporativi, da 20 tr.
- x2 - magazzini in linea, da 30 tr.
A nostra squadra di travagliu amichevule include quattru juniors, dui media è un senior. U so fondu di u tempu di travagliu mensile:
- ghjugnu:
4 * 150 = 600 чел * час, - mezu:
2 * 150 = 300 чел * час, - sentore:
150 чел * час.
Chì u primu junior dispunibile spende (0, 1, 2) ore nantu à u sviluppu è a implementazione di un situ di tipu (x10, x20, x30), mediu - (7, 15, 20), senior - (5, 10, 15). ) ore di u megliu tempu di a vostra vita.
Cum'è qualsiasi direttore normale, vulemu maximizà i profitti mensili. U primu passu per u successu hè di scrive a funzione objetiva value cum'è a quantità di ingressu da i prudutti pruduciuti per mese:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Questu ùn hè micca un errore; quandu cercate u massimu, a funzione objetiva hè minimizzata cù u signu oppostu.
U prossimu passu hè di pruibisce à i nostri impiegati di travaglià eccessivamente è di introduci restrizioni à l'ore di travagliu:
Cosa hè equivalente:
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]])
}Una restrizzione formale hè chì a pruduzzioni di u produttu deve esse solu pusitiva:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])È infine, l'assunzione più rosy hè chì per via di u prezzu bassu è di l'alta qualità, una fila di clienti soddisfatti hè sempre in linea per noi. Pudemu sceglie i volumi di produzzione mensili noi stessi, basatu nantu à risolve u prublema di ottimisazione restrittiva 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]Circundemu liberamente à numeri interi è calculemu a carica mensile di rematori cù una distribuzione ottima di i prudutti x = (8, 6, 3) :
- ghjugnu:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - mezu:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - sentore:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Conclusioni: per chì u direttore riceve u so massimu ben meritatu, hè ottimale per creà pagine di destinazione 8, siti di media 6 è magazzini 3 per mese. In questu casu, l'anzianu deve aratu senza guardà da a macchina, a carica di i medii serà apprussimatamente 2/3, i juniors menu di a mità.
cunchiusioni
L'articulu delinea e tecniche basiche per travaglià cù u pacchettu scipy.optimize, usatu per risolve i prublemi di minimizazione cundizionale. Personalmente aghju utilizatu scipy puramente per scopi accademici, chì hè per quessa chì l'esempiu datu hè di una natura cusì comica.
Una mansa di teoria è esempi virtuali ponu esse truvati, per esempiu, in u libru di IL Akulich "Programmazione matematica in esempi è prublemi". Applicazione più hardcore scipy.optimize per custruisce una struttura 3D da un inseme di imagine () pò esse vistu in .
A fonte principale di infurmazione hè quelli chì vulianu cuntribuisce à a traduzzione di sta è altre rùbbriche scipy Benvenuti à .
Спасибо per a participazione à a preparazione di a publicazione.
Source: www.habr.com