SciPy (se pronuncia como sai pai) es un paquete matemático basado en numpy, que también incluye bibliotecas en C y Fortran. Con SciPy, una sesión interactiva de Python se convierte en un entorno completo para el procesamiento de datos, similar a MATLAB, IDL, Octave, R o SciLab.
En este artículo, revisaremos las técnicas básicas de programación matemática: soluciones a problemas de optimización con restricciones para funciones escalares de múltiples variables utilizando el paquete scipy.optimize. Los algoritmos de optimización sin restricciones ya se han revisado en . Se puede obtener siempre una guía más detallada y actualizada sobre las funciones de scipy utilizando el comando help(), Shift+Tab o en .
Introducción
La interfaz general para resolver problemas tanto de optimización con restricciones como sin restricciones en el paquete scipy.optimize se proporciona mediante la función minimize(). Sin embargo, se sabe que no existe un método universal para resolver todos los problemas, por lo que la elección del método adecuado recae, como siempre, en el investigador.
El algoritmo de optimización adecuado se determina mediante el argumento de la función minimize(..., method="").
Para la optimización con restricciones de funciones de múltiples variables, están disponibles implementaciones de los siguientes métodos:
trust-constr— búsqueda de un mínimo local en un área de confianza. , ;SLSQP— programación cuadrática secuencial con restricciones, método de Newton para resolver el sistema de Lagrange. .TNC— Truncated Newton Constrained, número limitado de iteraciones, bueno para funciones no lineales con un gran número de variables independientes. .L-BFGS-B— método de la cuádruple Broyden–Fletcher–Goldfarb–Shanno, implementado con un menor consumo de memoria mediante la carga parcial de vectores de la matriz Hessiana. , .COBYLA— COBYLA (Constrained Optimization By Linear Approximation), optimización restringida con aproximación lineal (sin cálculo del gradiente). .
Dependiendo del método elegido, las condiciones y restricciones para resolver el problema se establecen de diferentes maneras:
- con un objeto de la clase
Boundspara los métodos L-BFGS-B, TNC, SLSQP, trust-constr; - con una lista
(min, max)para los mismos métodos L-BFGS-B, TNC, SLSQP, trust-constr; - con un objeto o lista de objetos
LinearConstraint,NonlinearConstraintpara los métodos COBYLA, SLSQP, trust-constr; - con un diccionario o lista de diccionarios
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}para los métodos COBYLA, SLSQP.
Plan del artículo:
1) Considerar la aplicación del algoritmo de optimización con restricciones en un área de confianza (method="trust-constr") con restricciones establecidas como objetos. Bounds, LinearConstraint, NonlinearConstraint ;
2) Considerar la programación secuencial mediante el método de los mínimos cuadrados (method=»SLSQP») con restricciones definidas como un diccionario {'type', 'fun', 'jac', 'args'};
3) Analizar un ejemplo de optimización de la producción a través de un estudio web.
Optimización condicional method=»trust-constr»
Implementación del método trust-constr se basa en para problemas con restricciones en forma de igualdad y en para problemas con restricciones en forma de desigualdad. Ambos métodos están implementados por algoritmos que buscan el mínimo local en un área de confianza y son adecuados para problemas a gran escala.
Planteamiento matemático del problema de búsqueda de mínimos en forma general:
Para las restricciones de igualdad estricta, el límite inferior se establece igual al límite superior 
.
Para la restricción unilateral, se establece el límite superior o inferior np.inf con el signo correspondiente.
Supongamos que es necesario encontrar el mínimo de la conocida función de Rosenbrock de dos variables:
Dichas restricciones son las siguientes sobre su dominio:
En nuestro caso hay una única solución en el punto 
, para el cual solo se cumplen la primera y la cuarta restricciones.
Pasemos por las restricciones de abajo hacia arriba y veamos cómo se pueden expresar en scipy.
Limitaciones 
y 
definir con el objeto Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Limitaciones 
y 
escribamos en forma lineal:
Definamos estas restricciones como un objeto LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])Y finalmente, la restricción no lineal en forma matricial:
Definimos la matriz Jacobiana para esta restricción y la combinación lineal de la matriz Hessiana con un vector arbitrario 
:
Ahora podemos definir la restricción no lineal como un objeto 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 el tamaño es grande, las matrices se pueden definir también en forma dispersa:
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 como un objeto 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)Cuando el cálculo de la matriz Hessiana 
requiere altos costos, se puede utilizar la clase Las siguientes estrategias están disponibles: BFGS y SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())La Hessiana también se puede calcular usando diferencias finitas:
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess='2-point')La matriz Jacobiana para las restricciones también se puede calcular usando diferencias finitas. Sin embargo, en este caso, la matriz Hessiana ya no se puede calcular con diferencias finitas. La Hessiana debe ser definida como una función o mediante la clase HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess=BFGS())La solución del problema de optimización es la siguiente:
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` condición de terminación satisfecha.
Número de iteraciones: 12, evaluaciones de función: 8, iteraciones CG: 7, optimalidad: 2.99e-09, violación de restricciones: 1.11e-16, tiempo de ejecución: 0.033 s.
[0.41494531 0.17010937]Si es necesario, la función que calcula la Hessiana se puede definir utilizando la clase 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 producto de la Hessiana y un vector arbitrario a través del parámetro 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)Alternativamente, las primeras y segundas derivadas de la función optimizada pueden ser calculadas aproximadamente. Por ejemplo, la Hessiana puede ser aproximada mediante una función SR1 (aproximación cuasi-Newton). El gradiente puede ser aproximado utilizando diferencias finitas.
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)Optimización condicional método='SLSQP'
El método SLSQP está diseñado para resolver problemas de minimización de funciones en forma de:
Donde 
y 
— conjuntos de índices de expresiones que describen restricciones en forma de igualdades o desigualdades. 
— conjuntos de límites inferiores y superiores para el dominio de definición de la función.
Las restricciones lineales y no lineales se describen en forma de diccionarios con claves tipo, fun y 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 búsqueda del mínimo se realiza de la siguiente 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)La optimización se ha terminado con éxito. (Modo de salida 0)
Valor actual de la función: 0.34271757499419825
Iteraciones: 4
Evaluaciones de función: 5
Evaluaciones de gradiente: 4
[0.41494475 0.1701105 ]Ejemplo de optimización
Con la transición hacia la quinta revolución tecnológica, consideremos la optimización de la producción a través del ejemplo de una agencia web que nos genera ingresos pequeños pero estables. Pongámonos en el lugar del director de una galería, donde se producen tres tipos de productos:
- x0 — landing pages que venden, desde 10 mil R.
- x1 — sitios web corporativos, desde 20 mil R.
- x2 — tiendas en línea, desde 30 mil R.
Nuestro amistoso equipo de trabajo consiste en cuatro juniors, dos mediadores y un senior. El fondo de su tiempo laboral al mes:
- juniors:
4 * 150 = 600 horas-persona, - mediadores:
2 * 150 = 300 horas-persona, - senior:
150 horas-persona.
Supongamos que para desarrollar y desplegar un sitio de tipo (x0, x1, x2), el primer junior que se encuentra necesita dedicar (10, 20, 30) horas, el mediador — (7, 15, 20), y el senior — (5, 10, 15) horas del mejor tiempo de su vida.
Como a cualquier director normal, nos gustaría maximizar las ganancias mensuales. El primer paso hacia el éxito es escribir la función objetivo value como la suma de los ingresos generados por la producción mensual:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Esto no es un error, al buscar el máximo la función objetivo se minimiza con un signo inverso.
El siguiente paso es prohibir a los empleados hacer horas extras e introducir restricciones sobre el fondo de tiempo laboral:
Lo cual es equivalente 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 restricción formal es que la producción debe ser únicamente positiva:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Y finalmente, la suposición más optimista: debido a nuestros precios bajos y alta calidad, siempre hay una cola de clientes satisfechos. Podemos elegir nosotros mismos los volúmenes de producción mensual, basándonos en la solución de un problema de optimización condicional 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]Redondeamos a números enteros y calculamos la carga mensual de los remeros en el escenario óptimo de producción x = (8, 6, 3) :
- juniors:
8 * 10 + 6 * 20 + 3 * 30 = 290 personas * hora; - mediadores:
8 * 7 + 6 * 15 + 3 * 20 = 206 personas * hora; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 personas * hora.
Conclusión: para que el director reciba su máximo merecido, lo óptimo es hacer al mes 8 landing pages, 6 sitios web medios y 3 tiendas. El senior debe trabajar sin descanso, la carga de los mid-levels será aproximadamente de 2/3, y la de los juniors menos de la mitad.
Conclusión
En el artículo se exponen las principales técnicas para trabajar con el paquete scipy.optimize, utilizadas para resolver problemas de minimización condicional. Personalmente, utilizo scipy solo con fines académicos, por lo que el ejemplo dado tiene un carácter humorístico.
Se puede encontrar mucha teoría y ejemplos en raros en el libro de I.L. Akulich "Programación matemática en ejemplos y tareas". Un uso más hardcore scipy.optimize para construir una estructura 3D a partir de un conjunto de imágenes () se puede ver en .
La fuente principal de información es , los interesados en contribuir a la traducción de esta y otras secciones scipy son bienvenidos a .
Gracias por su participación en la preparación de la publicación.
Fuente: habr.com
