SciPy (pronunciado sai pie) es un paquete matemático basado en numpy que también incluye bibliotecas C y Fortran. SciPy convierte su sesión interactiva de Python en un entorno completo de ciencia de datos como MATLAB, IDL, Octave, R o SciLab.
En este artículo, veremos las técnicas básicas de programación matemática: resolver problemas de optimización condicional para una función escalar de varias variables utilizando el paquete scipy.optimize. Los algoritmos de optimización sin restricciones ya se han discutido en . Siempre se puede obtener ayuda más detallada y actualizada sobre las funciones de scipy usando el comando help(), Shift+Tab o en .
introducción
La función proporciona una interfaz común para resolver problemas de optimización condicional y sin restricciones en el paquete scipy.optimize. minimize(). Sin embargo, se sabe que no existe un método universal para resolver todos los problemas, por lo que la elección de un método adecuado, como siempre, recae sobre los hombros del investigador.
El algoritmo de optimización apropiado se especifica utilizando el argumento de la función. minimize(..., method="").
Para la optimización condicional de una función de varias variables, se encuentran disponibles implementaciones de los siguientes métodos:
trust-constr— buscar un mínimo local en la región de confianza. , ;SLSQP— programación cuadrática secuencial con restricciones, método newtoniano para resolver el sistema de Lagrange. .TNC- Newton truncado Restringido, número limitado de iteraciones, bueno para funciones no lineales con una gran cantidad de variables independientes. .L-BFGS-B— un método del equipo Broyden–Fletcher–Goldfarb–Shanno, implementado con un consumo de memoria reducido debido a la carga parcial de vectores de la matriz de Hesse. , .COBYLA— MARE Optimización restringida por aproximación lineal, optimización restringida con aproximación lineal (sin cálculo de gradiente). .
Dependiendo del método elegido, las condiciones y restricciones para resolver el problema se establecen de forma diferente:
- objeto de clase
Boundspara métodos L-BFGS-B, TNC, SLSQP, confianza-constr; - lista
(min, max)para los mismos métodos L-BFGS-B, TNC, SLSQP, trust-constr; - un objeto o una lista de objetos
LinearConstraint,NonlinearConstraintpara COBYLA, SLSQP, métodos de control de confianza; - diccionario o lista de diccionarios
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}para COBYLA, métodos SLSQP.
Esquema del artículo:
1) Considere el uso de un algoritmo de optimización condicional en la región de confianza (método=”trust-constr”) con restricciones especificadas como objetos Bounds, LinearConstraint, NonlinearConstraint ;
2) Considere la programación secuencial utilizando el método de mínimos cuadrados (método = "SLSQP") con restricciones especificadas en forma de diccionario. {'type', 'fun', 'jac', 'args'};
3) Analizar un ejemplo de optimización de productos fabricados utilizando el ejemplo de un estudio web.
Método de optimización condicional="trust-constr"
Implementación del método. trust-constr Residencia en para problemas con restricciones de la forma de igualdad y en para problemas con restricciones en forma de desigualdades. Ambos métodos se implementan mediante algoritmos para encontrar un mínimo local en la región de confianza y son muy adecuados para problemas a gran escala.
Formulación matemática del problema de encontrar un mínimo en forma general:
Para restricciones de igualdad estrictas, el límite inferior se establece igual al límite superior .
Para una restricción unidireccional, el límite superior o inferior se establece np.inf con el signo correspondiente.
Sea necesario encontrar el mínimo de una función de Rosenbrock conocida de dos variables:
En este caso, se establecen las siguientes restricciones en su dominio de definición:
En nuestro caso, existe una solución única en el punto , para lo cual sólo son válidas las restricciones primera y cuarta.
Repasemos las restricciones de abajo hacia arriba y veamos cómo podemos escribirlas en scipy.
restricciones и Definámoslo usando el objeto Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])restricciones и Escribámoslo 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 una combinación lineal de la matriz de Hesse con un vector arbitrario :
Ahora podemos definir una 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 también se pueden especificar 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)Al calcular la matriz de Hesse requiere mucho esfuerzo, puedes usar una clase . Las siguientes estrategias están disponibles: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())El hessiano también se puede calcular mediante diferencias finitas:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')La matriz jacobiana para restricciones también se puede calcular utilizando diferencias finitas. Sin embargo, en este caso la matriz de Hesse no se puede calcular mediante diferencias finitas. El hessiano debe definirse como una función o utilizando la clase HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())La solución al problema de optimización se ve así:
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 es necesario, la función para calcular el hessiano se puede definir usando 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 del hessiano 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, se pueden aproximar la primera y segunda derivadas de la función que se está optimizando. Por ejemplo, el hessiano se puede aproximar usando la función SR1 (aproximación cuasi-newtoniana). El gradiente se puede aproximar mediante 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)Método de optimización condicional="SLSQP"
El método SLSQP está diseñado para resolver problemas de minimización de una función en la forma:
Donde и — conjuntos de índices de expresiones que describen restricciones en forma de igualdades o desigualdades. — conjuntos de límites inferior y superior para el dominio de definición de la función.
Las restricciones lineales y no lineales se describen en forma de diccionarios con claves. 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 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)Optimization terminated successfully. (Exit mode 0)
Current function value: 0.34271757499419825
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
[0.41494475 0.1701105 ]Ejemplo de optimización
En relación con la transición a la quinta estructura tecnológica, consideremos la optimización de la producción en el ejemplo de un estudio web, que nos aporta unos ingresos pequeños pero estables. Imaginémonos como el director de una galera que elabora tres tipos de productos:
- x0 - venta de páginas de destino, desde 10 tr.
- x1 - webs corporativas, desde 20 tr.
- x2 - tiendas online, desde 30 tr.
Nuestro amigable equipo de trabajo incluye cuatro juniors, dos middles y un senior. Su fondo mensual de tiempo de trabajo:
- Junios:
4 * 150 = 600 чел * час, - medios:
2 * 150 = 300 чел * час, - señor:
150 чел * час.
Deje que el primer junior disponible dedique (0, 1, 2) horas al desarrollo e implementación de un sitio del tipo (x10, x20, x30), medio - (7, 15, 20), senior - (5, 10, 15 ) horas del mejor momento de tu vida.
Como cualquier director normal, queremos maximizar las ganancias mensuales. El primer paso hacia el éxito es escribir la función objetivo. value como la cantidad de ingresos de los productos producidos por mes:
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 el signo opuesto.
El siguiente paso es prohibir a nuestros empleados trabajar demasiado e introducir restricciones en las horas de trabajo:
Qué es 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 restricción formal es que la producción del producto sólo debe ser positiva:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Y finalmente, la suposición más optimista es que, debido al bajo precio y la alta calidad, una cola de clientes satisfechos hace cola constantemente para recibirnos. Podemos elegir nosotros mismos los volúmenes de producción mensual, basándonos en resolver el problema de optimización restringida 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]Redondeemos libremente a números enteros y calculemos la carga mensual de remeros con una distribución óptima de productos. x = (8, 6, 3) :
- Junios:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - medios:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - señor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Conclusión: para que el director reciba su merecido máximo, lo óptimo es crear 8 páginas de destino, 6 sitios medianos y 3 tiendas por mes. En este caso, el mayor debe arar sin levantar la vista de la máquina, la carga de los medianos será aproximadamente 2/3, los jóvenes menos de la mitad.
Conclusión
El artículo describe las técnicas básicas para trabajar con el paquete. scipy.optimize, utilizado para resolver problemas de minimización condicional. Personalmente uso scipy con fines puramente académicos, por lo que el ejemplo dado tiene un carácter tan cómico.
Se pueden encontrar muchas teorías y ejemplos virtuales, por ejemplo, en el libro de I. L. Akulich "Programación matemática en ejemplos y problemas". Aplicación más dura scipy.optimize para construir una estructura 3D a partir de un conjunto de imágenes () se puede ver en .
La principal fuente de información es aquellos que deseen contribuir a la traducción de esta y otras secciones scipy Bienvenido a .
Gracias por participar en la preparación de la publicación.
Fuente: habr.com