SciPy, optimización con condiciones.

SciPy (pronunciado sai pie) es un paquete matemático basado en numpy que también incluye bibliotecas C y Fortran. Con SciPy, una sesión interactiva de Python se convierte en un entorno de procesamiento de datos completo como MATLAB, IDL, Octave, R o SciLab.

En este artículo, consideraremos 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 incondicional ya se consideran en el último artículo . Siempre se puede obtener ayuda más detallada y actualizada sobre funciones scipy utilizando el comando help (), Shift + Tab, o en la documentación oficial .

Introduccion

La interfaz minimize() los problemas de optimización condicional e incondicional en el paquete scipy.optimize. Sin embargo, se sabe que no existe un método universal para resolver todos los problemas, por lo tanto, como siempre, la elección de un método adecuado corresponde al investigador.

Se especifica un algoritmo de optimización adecuado utilizando el argumento para la función minimize(..., method="") .

Para la optimización condicional de la función de varias variables, están disponibles los siguientes métodos:

• trust-constr : busca un mínimo local en el área de confianza. Artículo de Wiki, artículo de Habr ;
• SLSQP : programación cuadrática secuencial con restricciones, el método newtoniano de resolución del sistema Lagrange. Artículo de Wiki .
• TNC - Newton Truncated Constrained, un número limitado de iteraciones, es bueno para funciones no lineales con una gran cantidad de variables independientes. Artículo de Wiki .
• L-BFGS-B : un método de los cuatro Broyden - Fletcher - Goldfarb - Shanno, implementado con un consumo de memoria reducido debido a la carga parcial de vectores de la matriz de Hesse. Artículo de Wiki , artículo de Habr .
• COBYLA - MARES Optimización restringida por aproximación lineal, optimización limitada con aproximación lineal (sin cálculo de gradiente). Artículo de Wiki .

Dependiendo del método elegido, las condiciones y restricciones para resolver el problema se establecen de manera diferente:

• un objeto de la clase Bounds para los métodos L-BFGS-B, TNC, SLSQP, trust-constr;
• la lista (min, max) para los mismos métodos L-BFGS-B, TNC, SLSQP, trust-constr;
• objeto o lista de objetos LinearConstraint , NonlinearConstraint para los métodos COBYLA, SLSQP, trust-constr;
• diccionario o lista de diccionarios {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} para los métodos COBYLA, SLSQP.

El bosquejo del artículo:
1) Considere la aplicación del algoritmo de optimización condicional en el área de confianza (método = "trust-constr") con las restricciones especificadas en forma de objetos Bounds , LinearConstraint , NonlinearConstraint ;
2) Considere la programación secuencial de mínimos cuadrados (método = "SLSQP") con las restricciones especificadas en forma de diccionario {'type', 'fun', 'jac', 'args'} ;
3) Desmontar un ejemplo de optimización de productos en el ejemplo de un estudio web.

Método de optimización condicional = "trust-constr"

La implementación del trust-constr basa en EQSQP para problemas con restricciones de la forma de igualdad y en TRIP para problemas con restricciones en forma de desigualdades. Ambos métodos se implementan mediante algoritmos para encontrar un mínimo local en el área de confianza y son muy adecuados para tareas a gran escala.

La formulación matemática del problema de búsqueda mínima en la forma general:

$$

$$\ min_x f (x)$$

$$

$$c ^ l \ leq c (x) \ leq c ^ u,$$

$$

$$x ^ l \ leq x \ leq x ^ u$$

Para restricciones de igualdad estrictas, el límite inferior se establece igual al superior $$$c ^ l_j = c ^ u_j$ .
Para una restricción unidireccional, np.inf establece el límite superior o inferior con el signo correspondiente.
Sea necesario encontrar el mínimo de la función conocida de Rosenbrock de dos variables:

$$

$$\ min_ {x_0, x_1} 100 (x_0 - x_1 ^ 2) ^ 2 + (1-x_0) ^ 2$$

En este caso, las siguientes restricciones se establecen en su dominio de definición:

$$

$$x_0 ^ 2 + x_1 \ leq 1$$

$$

$$x_0 ^ 2 - x_1 \ leq 1$$

$$

$$2x_0 + x_1 = 1$$

$$

$$x_0 + 2x_1 \ leq 1$$

$$

$$0 \ leq x_0 \ leq 1$$

$$

$$-0.5 \ leq x_1 \ leq 2.0$$

En nuestro caso, hay una solución única en este punto $$$[x_0, x_1] = [0.4149, 0.1701]$ para el cual solo las restricciones primera y cuarta son válidas.
Repasemos las restricciones de abajo hacia arriba y consideremos cómo se pueden escribir en scipy.
Limitaciones $$$0 \ leq x_0 \ leq 1$ y $$$-0.5 \ leq x_1 \ leq 2.0$ definido usando el objeto Bounds.

 from scipy.optimize import Bounds bounds = Bounds ([0, -0.5], [1.0, 2.0]) 

Limitaciones $$$x_0 + 2 x_1 \ leq 1$ y $$$2 x_0 + x_1 = 1$ escribimos en forma lineal:

$$

$$\ begin {bmatrix} - \ infty \\ 1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 & 2 \\ 2 & 1 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \ end {bmatrix } \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$$

Defina 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, una restricción no lineal en forma matricial:

$$

$$c (x) = \ begin {bmatrix} x_0 ^ 2 + x_1 \\ x_0 ^ 2 - x_1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$$

Definimos la matriz de Jacobi para esta restricción y la combinación lineal de la matriz de Hesse con un vector arbitrario. $$$v$ :

$$

$$J (x) = \ begin {bmatrix} 2x_0 & 1 \\ 2x_0 & -1 \ end {bmatrix}$$

$$

$$H (x, v) = \ sum \ limits_ {i = 0} ^ 1 v_i \ nabla ^ 2 c_i (x) = v_0 \ begin {bmatrix} 2 & 0 \\ 2 & 0 \ end {bmatrix} + v_1 \ begin {bmatrix} 2 y 0 \\ 2 y 0 \ end {bmatrix}$$

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) 

Cuando el cálculo de la matriz de Hesse $$$H (x, v)$ costoso, puede usar la clase HessianUpdateStrategy . 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 arpillera también se puede calcular usando diferencias finitas:

 nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point') 

La matriz de Jacobi para restricciones también se puede calcular utilizando diferencias finitas. Sin embargo, en este caso, la matriz de Hesse no se puede calcular utilizando diferencias finitas. Hessian 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 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` 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 de cálculo de Hesse 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 arpillera 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 derivadas primera y segunda de la función optimizada se pueden calcular aproximadamente. Por ejemplo, la arpillera se puede aproximar usando la función SR1 (aproximación cuasi-newtoniana). El gradiente puede ser aproximado por 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 el problema de minimizar la función en forma de:

$$

$$\ min_x f (x)$$

$$

$$c_j (x) = 0, j \ in \ mathcal {E}$$

$$

$$c_j (x) \ geq 0, j \ in \ mathcal {I}$$

$$

$$lb_i \ leq x_i \ leq ub_i, i = 1, ..., N$$

Donde $$$\ mathcal {E}$ y $$$\ mathcal {I}$ - conjuntos de índices de expresión que describen restricciones en forma de igualdades o desigualdades. $$$[lb_i, ub_i]$ - conjuntos de límites inferior y superior para el dominio de definición de función.

Las restricciones lineales y no lineales se describen en forma de diccionarios con las teclas type , 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 mínima 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 al quinto modo tecnológico, consideremos la optimización de la producción utilizando el ejemplo de un estudio web, que nos brinda un ingreso pequeño pero estable. Permítanos presentarnos como el director de una galería, que produce tres tipos de productos:

• x0 - venta de desembarques, desde 10 tr
• x1 - sitios corporativos, desde 20 tr
• x2 - tiendas en línea, desde 30 tr

Nuestro equipo de trabajo amigable incluye cuatro Juns, dos Middle y uno Senior. Fondo de su tiempo de trabajo durante un mes:

• Junio: 4 * 150 = 600 * ,
• medios: 2 * 150 = 300 * ,
• Senior: 150 *

Suponga que el primer junior que gasta en desarrollar y desplegar un sitio de tipo (x0, x1, x2) debe pasar (10, 20, 30) horas, medio - (7, 15, 20), senior - (5, 10, 15) horas de lo mejor tiempo de su vida.

Como cualquier director normal, queremos maximizar nuestras ganancias mensuales. El primer paso para el éxito es anotar el value función objetivo como la suma de los ingresos de 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 el signo opuesto.

El siguiente paso es prohibir el procesamiento a nuestros empleados e introducir restricciones en el fondo del tiempo de trabajo:

$$

$$\ begin {bmatrix} 10 y 20 y 30 \\ 7 y 15 y 20 \\ 5 y 10 y 15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ leq \ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix}$$

Que es equivalente:

$$

$$\ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix} - \ begin {bmatrix} 10 y 20 y 30 \\ 7 y 15 y 20 \\ 5 y 10 y 15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ geq 0$$

 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]]) } 

Restricción formal: el resultado solo debe ser positivo:

 bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf]) 

Y, por último, la suposición más optimista: debido al bajo precio y la alta calidad, una línea de clientes satisfechos se alinea constantemente para nosotros. Podemos elegir volúmenes de producción mensuales nosotros mismos, en función de la solución del 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 la piruleta al entero más cercano y calculamos la carga mensual de remeros con una distribución óptima de producción x = (8, 6, 3) :

• Junio: 8 * 10 + 6 * 20 + 3 * 30 = 290 * ;
• medios: 8 * 7 + 6 * 15 + 3 * 20 = 206 * ;
• Senior: 8 * 5 + 6 * 10 + 3 * 15 = 145 * .

Conclusión: para que el director obtenga su merecido máximo, es óptimo realizar 8 aterrizajes, 6 sitios medianos y 3 tiendas por mes. Al mismo tiempo, el señor debe arar sin separarse de la máquina, la carga del medio es de aproximadamente 2/3, menos de la mitad de junio.

Conclusión

El artículo describe las técnicas básicas para trabajar con el paquete scipy.optimize que se utilizan para resolver problemas de minimización condicional. Personalmente, uso scipy únicamente con fines académicos, por lo que este ejemplo es muy cómico.

Se pueden encontrar muchos ejemplos de teoría y winrar, por ejemplo, en el libro de I. L. Akulich "Programación matemática en ejemplos y problemas". Se puede encontrar una aplicación más hardcore de scipy.optimize para construir una estructura 3D para un conjunto de imágenes ( un artículo en el centro ) en el scipy-cookbook .

La principal fuente de información es docs.scipy.org. Si desea traducir esta y otras secciones de scipy a traducción, bienvenido a GitHub .

Gracias a mefistófeas por contribuir a esta publicación.