SciPy (pronunciato come sai pai) è un pacchetto matematico basato su numpy, che include anche librerie in C e Fortran. Con SciPy, una sessione interattiva di Python si trasforma in un ambiente di elaborazione dati completo, simile a MATLAB, IDL, Octave, R o SciLab.
In questo articolo esamineremo le principali tecniche di programmazione matematica: risoluzione di problemi di ottimizzazione vincolata per una funzione scalare di più variabili utilizzando il pacchetto scipy.optimize. Gli algoritmi di ottimizzazione senza vincoli sono già stati considerati in . È possibile ottenere una documentazione più dettagliata e aggiornata sulle funzioni di scipy in qualsiasi momento utilizzando il comando help(), Shift+Tab o in .
Introduzione
L'interfaccia generale per la risoluzione di problemi sia di ottimizzazione vincolata che non vincolata nel pacchetto scipy.optimize è fornita dalla funzione minimize(). Tuttavia, è noto che non esiste un metodo universale per risolvere tutti i problemi, quindi la scelta del metodo adeguato ricade sempre sulle spalle del ricercatore.
L'algoritmo di ottimizzazione appropriato è specificato mediante l'argomento della funzione minimize(..., method="").
Per l'ottimizzazione vincolata di funzioni di più variabili sono disponibili implementazioni dei seguenti metodi:
trust-constr— ricerca del minimo locale in un'area fidata. , ;SLSQP— programmazione quadratica sequenziale con vincoli, metodo di Newton per risolvere il sistema di Lagrange. .TNC— Truncated Newton Constrained, numero limitato di iterazioni, utile per funzioni non lineari con un elevato numero di variabili indipendenti. .L-BFGS-B— metodo Broyden–Fletcher–Goldfarb–Shanno, implementato con un consumo di memoria ridotto tramite il caricamento parziale dei vettori dalla matrice Hessiana. , .COBYLA— COBYLA Constrained Optimization By Linear Approximation, ottimizzazione vincolata con approssimazione lineare (senza calcolo del gradiente). .
A seconda del metodo scelto, le condizioni e i vincoli per risolvere il problema vengono definiti in modo diverso:
- oggetto della classe
Boundsper i metodi L-BFGS-B, TNC, SLSQP, trust-constr; - lista
(min, max)per gli stessi metodi L-BFGS-B, TNC, SLSQP, trust-constr; - oggetto o lista di oggetti
LinearConstraint,NonlinearConstraintper i metodi COBYLA, SLSQP, trust-constr; - dizionario o lista di dizionari
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}per i metodi COBYLA, SLSQP.
Piano dell'articolo:
1) Considerare l'applicazione dell'algoritmo di ottimizzazione condizionata nell'area di fiducia (method=»trust-constr») con vincoli definiti come oggetti Bounds, LinearConstraint, NonlinearConstraint ;
2) Esaminare la programmazione sequenziale mediante il metodo dei minimi quadrati (method=»SLSQP») con vincoli definiti come dizionario {'type', 'fun', 'jac', 'args'};
3) Analizzare un esempio di ottimizzazione della produzione nel contesto di uno studio web.
Ottimizzazione condizionata method=»trust-constr»
Implementazione del metodo trust-constr è basata su per problemi con vincoli di uguaglianza e su per problemi con vincoli di disuguaglianza. Entrambi i metodi sono implementati con algoritmi di ricerca del minimo locale nell'area di fiducia e si adattano bene a problemi su larga scala.
Formulazione matematica del problema di ricerca del minimo in forma generale:
Per vincoli di uguaglianza rigorosa, il limite inferiore è fissato uguale al limite superiore 
.
Per vincoli unilaterali, il limite superiore o inferiore è fissato np.inf con il segno corrispondente.
Si supponga di dover trovare il minimo della nota funzione di Rosenbrock di due variabili:
Siano dati i seguenti vincoli sul suo dominio di definizione:
Nel nostro caso c'è una sola soluzione nel punto 
, per il quale sono validi solo il primo e il quarto vincolo.
Esaminiamo i vincoli dal basso verso l'alto e vediamo come possiamo scriverli in scipy.
Limitazioni 
e 
definiamoli utilizzando l'oggetto Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Limitazioni 
e 
scriviamoli in forma lineare:
Definiamo questi vincoli come oggetto LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])E infine il vincolo non lineare in forma matriciale:
Definiamo la matrice di Jacobi per questo vincolo e una combinazione lineare della matrice Hessiana con un vettore arbitrario. 
:
Ora possiamo definire il vincolo non lineare come oggetto 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 le dimensioni sono grandi, le matrici possono essere definite anche 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)oppure come oggetto 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)Quando il calcolo della matrice Hessiana 
richiede risorse significative, è possibile utilizzare la classe . Sono disponibili le seguenti strategie: BFGS e SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())La matrice Hessiana può anche essere calcolata tramite differenze finite:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')La matrice Jacobiana per i vincoli può essere calcolata anche tramite differenze finite. Tuttavia, in questo caso, la matrice Hessiana non può più essere calcolata con le differenze finite. La Hessiana deve essere definita come funzione o tramite la classe HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())La soluzione del problema di ottimizzazione è la seguente:
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` condizione di terminazione soddisfatta.
Numero di iterazioni: 12, valutazioni della funzione: 8, iterazioni CG: 7, ottimalità: 2.99e-09, violazione dei vincoli: 1.11e-16, tempo di esecuzione: 0.033 s.
[0.41494531 0.17010937]Se necessario, la funzione per calcolare l' Hessiano può essere definita utilizzando la classe 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 il prodotto dell'Hessiano e di un vettore arbitrario tramite il parametro 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, le prime e seconde derivate della funzione ottimizzata possono essere calcolate approssimativamente. Ad esempio, l'Hessiano può essere approssimato usando la funzione SR1 (approssimazione quasi-Newton). Il gradiente può essere approssimato tramite 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)Ottimizzazione condizionata method=»SLSQP»
Il metodo SLSQP è progettato per risolvere problemi di minimizzazione di una funzione sotto forma di:
Dove 
e 
— un insieme di indici di espressioni che descrivono le vincoli in forma di uguaglianze o disuguaglianze. 
— un insieme di limiti inferiori e superiori per il dominio della funzione.
I vincoli lineari e non lineari sono descritti sotto forma di dizionari con le chiavi type, divertimento e 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])
}La ricerca del minimo viene effettuata nel seguente modo:
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)Ottimizzazione terminata con successo. (Exit mode 0)
Valore attuale della funzione: 0.34271757499419825
Iterazioni: 4
Valutazioni della funzione: 5
Valutazioni del gradiente: 4
[0.41494475 0.1701105 ]Esempio di ottimizzazione
Con il passaggio al quinto paradigma tecnologico, consideriamo l'ottimizzazione della produzione prendendo come esempio uno studio web che ci offre un guadagno modesto ma costante. Immaginiamo di essere il direttore di una galleria che produce tre tipi di prodotti:
- x0 — landing page di vendita, a partire da 10.000 €.
- x1 — siti aziendali, a partire da 20.000 €.
- x2 — negozi online, a partire da 30.000 €.
Il nostro affiatato team di lavoro è composto da quattro junior, due mid-level e un senior. Il monte ore lavorativo del mese:
- junior:
4 * 150 = 600 ore * persona, - mid-level:
2 * 150 = 300 ore * persona, - senior:
150 ore * persona.
Supponiamo che per lo sviluppo e il deploy di un sito del tipo (x0, x1, x2) un junior qualsiasi debba impiegare (10, 20, 30) ore, un mid-level — (7, 15, 20), e un senior — (5, 10, 15) ore del miglior tempo della sua vita.
Come a qualsiasi normale direttore, vogliamo massimizzare il profitto mensile. Il primo passo verso il successo è scrivere la funzione obiettivo value come somma dei ricavi derivanti dai prodotti realizzati nel mese:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Non è un errore, nella ricerca del massimo la funzione obiettivo viene minimizzata con il segno opposto.
Il passo successivo è vietare la produzione e introdurre limiti sul fondo orario lavorativo:
Cosa equivale a:
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]])
}Il vincolo formale è che la produzione deve essere solo positiva:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])E infine, l'assunzione più ottimista è che a causa del basso costo e dell'alta qualità, abbiamo sempre una fila di clienti soddisfatti. Possiamo scegliere noi stessi i volumi di produzione mensili, basandoci sulla soluzione del problema di ottimizzazione vincolata con 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]Arrotondiamo all'intero più vicino e calcoliamo il carico mensile dei rematori con la migliore disposizione della produzione x = (8, 6, 3) :
- junior:
8 * 10 + 6 * 20 + 3 * 30 = 290 persone * ora; - mid-level:
8 * 7 + 6 * 15 + 3 * 20 = 206 persone * ora; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 persone * ora.
Risultato: affinché il direttore riceva il massimo meritato, è ottimale realizzare 8 landing page, 6 siti medi e 3 negozi al mese. Il senior, nel frattempo, deve lavorare incessantemente, e il carico dei mid deve essere circa 2/3, mentre quello dei junior meno della metà.
Conclusione
L'articolo espone le principali tecniche di lavoro con il pacchetto scipy.optimize, utilizzate per risolvere problemi di minimizzazione condizionale. Personalmente utilizzo scipy puramente a scopi accademici, quindi l'esempio fornito ha un carattere umoristico.
Molta teoria e esempi rari possono essere trovati, ad esempio, nel libro di I.L. Akulich «Programmazione matematica in esempi e problemi». Applicazioni più hardcore scipy.optimize per la costruzione di strutture 3D basate su una serie di immagini () possono essere consultate in .
La fonte principale di informazioni è , chi desidera contribuire alla traduzione di questa e altre sezioni scipy è il benvenuto su .
Grazie per la partecipazione alla preparazione della pubblicazione.
Fonte: habr.com
