SciPy (pronunciato sai pie) è un pacchetto matematico basato su Numpy che include anche le librerie C e Fortran. SciPy trasforma la tua sessione interattiva di Python in un ambiente completo di scienza dei dati come MATLAB, IDL, Octave, R o SciLab.
In questo articolo esamineremo le tecniche di base della programmazione matematica, risolvendo problemi di ottimizzazione condizionale per una funzione scalare di più variabili utilizzando il pacchetto scipy.optimize. Gli algoritmi di ottimizzazione non vincolata sono già stati discussi in . È sempre possibile ottenere un aiuto più dettagliato e aggiornato sulle funzioni scipy utilizzando il comando help(), Shift+Tab o in .
Introduzione
Un'interfaccia comune per risolvere problemi di ottimizzazione condizionale e 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, come sempre, ricade sulle spalle del ricercatore.
L'algoritmo di ottimizzazione appropriato viene specificato utilizzando l'argomento della funzione minimize(..., method="").
Per l'ottimizzazione condizionale di una funzione di più variabili sono disponibili implementazioni dei seguenti metodi:
trust-constr— ricerca di un minimo locale nella regione di fiducia. , ;SLSQP— programmazione quadratica sequenziale con vincoli, metodo newtoniano per la risoluzione del sistema di Lagrange. .TNC- Newton troncato Vincolato, numero limitato di iterazioni, ottimo per funzioni non lineari con un gran numero di variabili indipendenti. .L-BFGS-B— un metodo del team Broyden–Fletcher–Goldfarb–Shanno, implementato con un ridotto consumo di memoria dovuto al caricamento parziale di vettori dalla matrice dell'Assia. , .COBYLA— Ottimizzazione vincolata MARE per approssimazione lineare, ottimizzazione vincolata con approssimazione lineare (senza calcolo del gradiente). .
A seconda del metodo scelto, le condizioni e le restrizioni per la risoluzione del problema sono impostate in modo diverso:
- oggetto di classe
Boundsper i metodi L-BFGS-B, TNC, SLSQP, trust-constr; - la lista
(min, max)per gli stessi metodi L-BFGS-B, TNC, SLSQP, trust-constr; - un oggetto o un elenco di oggetti
LinearConstraint,NonlinearConstraintper COBYLA, SLSQP, metodi di costruzione della fiducia; - dizionario o elenco di dizionari
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}per COBYLA, metodi SLSQP.
Descrizione dell'articolo:
1) Considerare l'uso di un algoritmo di ottimizzazione condizionale nella regione di fiducia (method="trust-constr") con vincoli specificati come oggetti Bounds, LinearConstraint, NonlinearConstraint ;
2) Considerare la programmazione sequenziale utilizzando il metodo dei minimi quadrati (metodo = "SLSQP") con restrizioni specificate sotto forma di dizionario {'type', 'fun', 'jac', 'args'};
3) Analizzare un esempio di ottimizzazione dei prodotti fabbricati utilizzando l'esempio di uno studio web.
Metodo di ottimizzazione condizionale="trust-constr"
Implementazione del metodo trust-constr basato su per problemi con vincoli della forma di uguaglianza e via per problemi con vincoli sotto forma di disuguaglianze. Entrambi i metodi sono implementati da algoritmi per trovare un minimo locale nella regione di confidenza e sono adatti a problemi su larga scala.
Formulazione matematica del problema di trovare un minimo in forma generale:
Per vincoli di uguaglianza rigorosi, il limite inferiore è fissato uguale al limite superiore .
Per un vincolo unidirezionale, viene impostato il limite superiore o inferiore np.inf con il segno corrispondente.
Sia necessario trovare il minimo di una funzione di Rosenbrock nota di due variabili:
In questo caso, al suo ambito di definizione vengono poste le seguenti restrizioni:
Nel nostro caso esiste una soluzione unica , per il quale valgono solo la prima e la quarta restrizione.
Esaminiamo le restrizioni dal basso verso l'alto e vediamo come possiamo scriverle in Scipy.
Restrizioni и definiamolo utilizzando l'oggetto Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restrizioni и Scriviamolo in forma lineare:
Definiamo questi vincoli come un 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 Jacobiana per questo vincolo e una combinazione lineare della matrice Hessiana con un vettore arbitrario :
Ora possiamo definire un 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 specificate 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:
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 si calcola la matrice dell'Assia richiede molto impegno, puoi usare una classe . Sono disponibili le seguenti strategie: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())L'Assia può anche essere calcolato utilizzando le 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 utilizzando le differenze finite. Tuttavia, in questo caso la matrice dell'Assia non può essere calcolata utilizzando le differenze finite. L'Hessian 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 al problema di ottimizzazione è simile alla 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` 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 necessario, la funzione per il calcolo dell'Hessiana 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'Assia e un vettore arbitrario attraverso 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, è possibile approssimare la derivata prima e la seconda della funzione da ottimizzare. Ad esempio, l'Assia può essere approssimato utilizzando la funzione SR1 (approssimazione quasi newtoniana). Il gradiente può essere approssimato mediante 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)Metodo di ottimizzazione condizionale="SLSQP"
Il metodo SLSQP è progettato per risolvere problemi di minimizzazione di una funzione nella forma:
Где и — insiemi di indici di espressioni che descrivono restrizioni 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 и 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 si effettua come segue:
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 ]Esempio di ottimizzazione
In connessione con il passaggio alla quinta struttura tecnologica, consideriamo l'ottimizzazione della produzione utilizzando l'esempio di uno studio web, che ci porta un reddito piccolo ma stabile. Immaginiamoci come il direttore di una cucina che produce tre tipologie di prodotti:
- x0 - vendita di pagine di destinazione, da 10 tr.
- x1 - siti web aziendali, da 20 tr.
- x2 - negozi online, da 30 tr.
Il nostro amichevole team di lavoro comprende quattro junior, due intermedi e uno senior. Il loro fondo mensile per l’orario di lavoro:
- giugno:
4 * 150 = 600 чел * час, - medi:
2 * 150 = 300 чел * час, - signore:
150 чел * час.
Consentire al primo junior disponibile di dedicare (0, 1, 2) ore allo sviluppo e all'implementazione di un sito di tipo (x10, x20, x30), medio - (7, 15, 20), senior - (5, 10, 15 ) ore del momento più bello della tua vita.
Come ogni normale direttore, vogliamo massimizzare i profitti mensili. Il primo passo verso il successo è scrivere la funzione obiettivo value come importo del reddito derivante dai prodotti fabbricati al mese:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Questo non è un errore; nella ricerca del massimo la funzione obiettivo viene minimizzata con il segno opposto.
Il prossimo passo è vietare ai nostri dipendenti di lavorare troppo e introdurre restrizioni sull’orario di lavoro:
Cosa è 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 restrizione formale è che l'output del prodotto deve essere solo positivo:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])E infine, l'ipotesi più rosea è che a causa del prezzo basso e dell'alta qualità, per noi si forma costantemente una fila di clienti soddisfatti. Possiamo scegliere noi stessi i volumi di produzione mensile, in base alla risoluzione del problema di ottimizzazione vincolata 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 liberamente ai numeri interi e calcoliamo il carico mensile dei vogatori con una distribuzione ottimale dei prodotti x = (8, 6, 3) :
- giugno:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - medi:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - signore:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Conclusione: affinché il direttore riceva il suo meritato massimo, è ottimale creare 8 landing page, 6 siti di medie dimensioni e 3 negozi al mese. In questo caso, il senior deve arare senza alzare lo sguardo dalla macchina, il carico dei medi sarà di circa 2/3, dei junior meno della metà.
conclusione
L'articolo descrive le tecniche di base per lavorare con il pacchetto scipy.optimize, utilizzato per risolvere problemi di minimizzazione condizionale. Personalmente utilizzo scipy puramente per scopi accademici, motivo per cui l'esempio fornito è di natura così comica.
Molte teorie ed esempi virtuali si possono trovare, ad esempio, nel libro di I.L. Akulich “Programmazione matematica in esempi e problemi”. Applicazione più hardcore scipy.optimize per costruire una struttura 3D da una serie di immagini () può essere visualizzato in .
La principale fonte di informazione è coloro che desiderano contribuire alla traduzione di questa e di altre sezioni scipy Benvenuto a .
Grazie per la partecipazione alla preparazione della pubblicazione.
Fonte: habr.com