SciPy (pronunciat sai pie) és un paquet de matemàtiques basat en numpy que també inclou biblioteques C i Fortran. SciPy converteix la vostra sessió interactiva de Python en un entorn complet de ciència de dades com MATLAB, IDL, Octave, R o SciLab.
En aquest article, veurem les tècniques bàsiques de programació matemàtica: resoldre problemes d'optimització condicional per a una funció escalar de diverses variables mitjançant el paquet scipy.optimize. Els algorismes d'optimització sense restriccions ja s'han parlat a . Sempre es pot obtenir ajuda més detallada i actualitzada sobre les funcions scipy mitjançant l'ordre help(), Maj+Tab o a .
Introducció
La funció proporciona una interfície comuna per resoldre problemes d'optimització condicionals i no restringits al paquet scipy.optimize. minimize(). Tanmateix, se sap que no existeix un mètode universal per resoldre tots els problemes, per la qual cosa l'elecció d'un mètode adequat, com sempre, recau sobre les espatlles de l'investigador.
L'algorisme d'optimització adequat s'especifica mitjançant l'argument de la funció minimize(..., method="").
Per a l'optimització condicional d'una funció de diverses variables, hi ha disponibles implementacions dels mètodes següents:
trust-constr— cercar un mínim local a la regió de confiança. , ;SLSQP— programació quadràtica seqüencial amb restriccions, mètode newtonià per resoldre el sistema de Lagrange. .TNC- Newton truncat Restringit, nombre limitat d'iteracions, bo per a funcions no lineals amb un gran nombre de variables independents. .L-BFGS-B— un mètode de l'equip Broyden–Fletcher–Goldfarb–Shanno, implementat amb un consum de memòria reduït a causa de la càrrega parcial de vectors de la matriu de Hessian. , .COBYLA— Optimització restringida MARE per aproximació lineal, optimització restringida amb aproximació lineal (sense càlcul de gradients). .
Depenent del mètode escollit, les condicions i les restriccions per resoldre el problema s'estableixen de manera diferent:
- objecte de classe
Boundsper als mètodes L-BFGS-B, TNC, SLSQP, trust-constr; - la llista
(min, max)per als mateixos mètodes L-BFGS-B, TNC, SLSQP, trust-constr; - un objecte o una llista d'objectes
LinearConstraint,NonlinearConstraintper a COBYLA, SLSQP, mètodes de construcció de confiança; - diccionari o llista de diccionaris
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}per als mètodes COBYLA, SLSQP.
Esquema de l'article:
1) Considereu l'ús d'un algorisme d'optimització condicional a la regió de confiança (method="trust-constr") amb restriccions especificades com a objectes Bounds, LinearConstraint, NonlinearConstraint ;
2) Considereu la programació seqüencial utilitzant el mètode dels mínims quadrats (mètode = "SLSQP") amb restriccions especificades en forma de diccionari {'type', 'fun', 'jac', 'args'};
3) Analitzar un exemple d'optimització de productes fabricats utilitzant l'exemple d'un estudi web.
Mètode d'optimització condicional="trust-constr"
Implementació del mètode trust-constr basat en per a problemes amb limitacions de la forma d'igualtat i així per a problemes amb restriccions en forma de desigualtats. Ambdós mètodes estan implementats per algorismes per trobar un mínim local a la regió de confiança i són molt adequats per a problemes a gran escala.
Formulació matemàtica del problema de trobar un mínim en forma general:
Per a restriccions d'igualtat estrictes, el límit inferior s'estableix igual al límit superior 
.
Per a una restricció unidireccional, s'estableix el límit superior o inferior np.inf amb el signe corresponent.
Sigui necessari trobar el mínim d'una funció de Rosenbrock coneguda de dues variables:
En aquest cas, s'estableixen les restriccions següents al seu domini de definició:
En el nostre cas, hi ha una solució única en aquest punt 
, per al qual només són vàlides les restriccions primera i quarta.
Repassem les restriccions de baix a dalt i mirem com les podem escriure en scipy.
Restriccions 
и 
anem a definir-lo mitjançant l'objecte Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restriccions 
и 
Escrivim-ho en forma lineal:
Definim aquestes restriccions com a objecte LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])I finalment la restricció no lineal en forma de matriu:
Definim la matriu jacobiana per a aquesta restricció i una combinació lineal de la matriu hessiana amb un vector arbitrari 
:
Ara podem definir una restricció no lineal com a objecte 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)Si la mida és gran, les matrius també es poden especificar en forma escassa:
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 com a objecte 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)Quan es calcula la matriu hessiana 
requereix molt d'esforç, podeu utilitzar una classe . Les estratègies següents estan disponibles: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())El Hessià també es pot calcular mitjançant diferències finites:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')La matriu jacobiana per a les restriccions també es pot calcular utilitzant diferències finites. Tanmateix, en aquest cas, la matriu de Hess no es pot calcular mitjançant diferències finites. El Hessian s'ha de definir com una funció o utilitzant la classe HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())La solució al problema d'optimització és la següent:
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]Si cal, la funció per calcular l'hessià es pot definir mitjançant 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 el producte de l'hessià i un vector arbitrari a través del paràmetre 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)Alternativament, es poden aproximar les derivades primera i segona de la funció que s'està optimitzant. Per exemple, el Hessià es pot aproximar mitjançant la funció SR1 (aproximació quasi newtoniana). El gradient es pot aproximar mitjançant diferències finites.
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)Mètode d'optimització condicional="SLSQP"
El mètode SLSQP està dissenyat per resoldre problemes de minimització d'una funció en la forma:
On 
и 
— conjunts d'índexs d'expressions que descriuen restriccions en forma d'igualtats o desigualtats. 
— conjunts de límits inferior i superior per al domini de definició de la funció.
Les restriccions lineals i no lineals es descriuen en forma de diccionaris amb tecles 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 cerca del mínim es realitza de la següent manera:
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 ]Exemple d'optimització
En relació amb la transició a la cinquena estructura tecnològica, mirem l'optimització de la producció amb l'exemple d'un estudi web, que ens aporta uns ingressos petits però estables. Imaginem-nos com a directors d'una galera que produeix tres tipus de productes:
- x0 - venda de pàgines de destinació, a partir de 10 tr.
- x1 - llocs web corporatius, a partir de 20 tr.
- x2 - botigues en línia, a partir de 30 tr.
El nostre amable equip de treball està format per quatre júniors, dos mitjans i un sènior. El seu fons de temps de treball mensual:
- junys:
4 * 150 = 600 чел * час, - mitjans:
2 * 150 = 300 чел * час, - senyor:
150 чел * час.
Permet que el primer júnior disponible passi (0, 1, 2) hores en el desenvolupament i el desplegament d'un lloc de tipus (x10, x20, x30), mitjà - (7, 15, 20), sènior - (5, 10, 15). ) hores del millor moment de la teva vida.
Com qualsevol director normal, volem maximitzar els beneficis mensuals. El primer pas per a l'èxit és escriure la funció objectiu value com a quantitat d'ingressos dels productes produïts al mes:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Això no és un error en cercar el màxim, la funció objectiu es minimitza amb el signe contrari.
El següent pas és prohibir als nostres empleats l'excés de treball i introduir restriccions a l'horari de treball:
Què és equivalent:
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 restricció formal és que la producció del producte només ha de ser positiva:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])I, finalment, la hipòtesi més optimista és que a causa del baix preu i de l'alta qualitat, una cua de clients satisfets ens està constantment fent cua. Podem triar nosaltres mateixos els volums de producció mensuals, basant-nos en la resolució del problema d'optimització restringida 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]Arrodonim a nombres sencers i calculem la càrrega de treball mensual dels remers amb una distribució òptima dels productes x = (8, 6, 3) :
- junys:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - mitjans:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senyor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Conclusió: perquè el director rebi el seu màxim merescut, és òptim crear 8 landing pages, 6 llocs de mida mitjana i 3 botigues al mes. En aquest cas, el sènior ha de llaurar sense aixecar la vista de la màquina, la càrrega dels mitjans serà aproximadament de 2/3, els júniors menys de la meitat.
Conclusió
L'article descriu les tècniques bàsiques per treballar amb el paquet scipy.optimize, utilitzat per resoldre problemes de minimització condicional. Personalment faig servir scipy amb finalitats purament acadèmics, per això l'exemple donat és de caràcter tan còmic.
Es poden trobar una gran quantitat de teoria i exemples virtuals, per exemple, al llibre d'I.L Akulich "La programació matemàtica en exemples i problemes". Aplicació més hardcore scipy.optimize per construir una estructura 3D a partir d'un conjunt d'imatges () es pot veure a .
La principal font d'informació és aquells que vulguin contribuir a la traducció d'aquesta i d'altres seccions scipy Benvingut a .
Gràcies per participar en l'elaboració de la publicació.
Font: www.habr.com
