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, come MATLAB, IDL, Octave, R o SciLab.
In questo articolo esploreremo le principali tecniche di programmazione matematica — soluzioni per problemi di ottimizzazione vincolata per una funzione scalare di più variabili utilizzando il pacchetto scipy.optimize. Gli algoritmi di ottimizzazione non vincolata sono già stati trattati in . Una documentazione più dettagliata e aggiornata sulle funzioni di scipy è sempre disponibile tramite il comando help(), Shift+Tab o in .
Introduzione
L'interfaccia generale per risolvere sia problemi 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 di un metodo adeguato ricade sempre sulle spalle del ricercatore.
L'algoritmo di ottimizzazione appropriato è definito tramite 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 di fiducia. , ;SLSQP— programmazione quadratica sequenziale con vincoli, metodo di Newton per risolvere il sistema di Lagrange. .TNC— Truncated Newton Constrained, numero limitato di iterazioni, buono per funzioni non lineari con un alto numero di variabili indipendenti. .L-BFGS-B— metodo della quadrupla Broyden–Fletcher–Goldfarb–Shanno, implementato con un ridotto consumo di memoria grazie al caricamento parziale dei vettori dalla matrice Hessiana. , .COBYLA— KOBYLA 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 sono definiti in modi diversi:
- oggetto della classe
Boundsper i metodi L-BFGS-B, TNC, SLSQP, trust-constr; - lista
(min, max)per questi 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 vincolata nella regione di fiducia (method=»trust-constr») con vincoli definiti sotto forma di oggetti Bounds, LinearConstraint, NonlinearConstraint ;
2) Considerare la programmazione sequenziale mediante metodo dei minimi quadrati (method=»SLSQP») con vincoli definiti in forma di dizionario {'type', 'fun', 'jac', 'args'};
3) Analizzare un esempio di ottimizzazione della produzione in un'agenzia web.
Ottimizzazione vincolata method=»trust-constr»
L'implementazione del metodo trust-constr si basa su per problemi con vincoli di uguaglianza e su per problemi con vincoli di disuguaglianza. Entrambi i metodi sono implementati mediante algoritmi di ricerca del minimo locale nella regione di fiducia e sono adatti per problemi su larga scala.
La formulazione matematica del problema di ricerca del minimo in forma generale:
Per vincoli di uguaglianza rigorosa, il limite inferiore è impostato uguale al limite superiore 
.
Per un vincolo unilaterale, il limite superiore o inferiore è impostato np.inf con il segno corrispondente.
Supponiamo di dover trovare il minimo della nota funzione di Rosenbrock di due variabili:
In questo caso, sono imposti i seguenti vincoli sulla sua area di definizione:
Nel nostro caso, c'è una sola soluzione nel punto 
, che soddisfa solo il primo e il quarto vincolo.
Esaminiamo i vincoli dall'alto verso il basso e vediamo come possono essere registrati 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 
scriviamo 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 un vincolo non lineare in forma matriciale:
Definiamo la matrice Jacobiana 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 la dimensione è grande, 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)o come oggetto LinearOperator:
da 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 elevati costi, è possibile utilizzare la classe . Sono disponibili le seguenti strategie: BFGS e SR1.
da scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())L'Hessiano può anche essere calcolato usando le differenze finite:
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess='2-point')La matrice Jacobiana per i vincoli può essere calcolata anche usando le differenze finite. Tuttavia, in questo caso l'Hessiano non può più essere calcolato con le differenze finite. L'Hessiano deve essere definito come una funzione o utilizzando la classe HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess=BFGS())La soluzione del problema di ottimizzazione appare come segue:
da scipy.optimize import minimize
da 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, optimalità: 2.99e-09, violazione del vincolo: 1.11e-16, tempo di esecuzione: 0.033 s.
[0.41494531 0.17010937]Se necessario, la funzione per calcolare l'Hessiano può essere definita usando 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-newtoniana). Il gradiente può essere approssimato con le differenze finite.
da 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 della funzione nella forma:
Dove 
e 
— insiemi di indici di espressioni che descrivono vincoli sotto forma di uguaglianze o disuguaglianze. 
— insiemi di limiti inferiori e superiori per il dominio di definizione della funzione.
I vincoli lineari e non lineari sono descritti sotto forma di dizionari con chiavi type, fun 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 avviene 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. (Modalità di uscita 0)
Valore attuale della funzione: 0.34271757499419825
Iterazioni: 4
Valutazioni della funzione: 5
Valutazioni del gradiente: 4
[0.41494475 0.1701105 ]Esempio di ottimizzazione
In relazione al passaggio al quinto paradigma tecnologico, consideriamo l'ottimizzazione della produzione prendendo come esempio uno studio web che ci fornisce un reddito modesto ma stabile. Immaginiamo di essere il direttore di una galleria, dove vengono prodotti tre tipi di prodotti:
- x0 — landing page di vendita, a partire da 10.000.
- x1 — siti web 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 e un senior. Il loro tempo di lavoro disponibile per mese è:
- junior:
4 * 150 = 600 ore * uomo, - mid:
2 * 150 = 300 ore * uomo, - senior:
150 ore * uomo.
Supponiamo che per sviluppare e rilasciare un sito di tipo (x0, x1, x2) il primo junior a disposizione debba impiegare (10, 20, 30) ore, il mid — (7, 15, 20), il senior — (5, 10, 15) ore del miglior tempo della sua vita.
Come a qualsiasi normale direttore, desideriamo massimizzare il profitto mensile. Il primo passo verso il successo è scrivere la funzione obiettivo value come la somma dei ricavi dei prodotti realizzati durante il 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 segno opposto.
Il passo successivo è vietare il lavoro straordinario ai dipendenti e introdurre vincoli sul tempo di lavoro disponibile:
Che 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]])
}La limitazione 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, grazie al prezzo basso e all'alta qualità, si crea continuamente una fila di clienti soddisfatti. Possiamo scegliere autonomamente i volumi di produzione mensili, basandoci sulla soluzione del problema di ottimizzazione condizionata 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 in modo impreciso ai numeri interi e calcoliamo il carico mensile degli rematori nella situazione ottimale di produzione x = (8, 6, 3) :
- junior:
8 * 10 + 6 * 20 + 3 * 30 = 290 persone * ora; - mid:
8 * 7 + 6 * 15 + 3 * 20 = 206 persone * ora; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 persone * ora.
Conclusione: per far sì che il direttore riceva il suo meritato massimo, è ottimale realizzare 8 landing page, 6 siti medi e 3 negozi al mese. Il senior, nel frattempo, deve lavorare incessantemente, il carico dei mid-level sarà di circa 2/3, mentre per i junior sarà meno della metà.
Conclusione
L'articolo espone le principali tecniche di lavoro con il pacchetto scipy.optimize, utilizzate per risolvere problemi di minimizzazione condizionata. Personalmente utilizzo scipy esclusivamente a fini accademici, quindi l'esempio fornito ha un carattere piuttosto giocoso.
Molta teoria e esempi in WinRAR possono essere trovati, ad esempio, nel libro di I.L. Akulich "Programmazione matematica in esempi e compiti". Un'applicazione più hardcore scipy.optimize per costruire strutture 3D da un insieme di immagini () può essere vista in .
La principale fonte di informazioni è , chi desidera contribuire alla traduzione di questa e di altre sezioni scipy è benvenuto su .
Grazie per la partecipazione nella preparazione della pubblicazione.
Fonte: habr.com
