SciPy (pronúnciase sai pie) é un paquete matemático baseado en numpy que tamén inclúe bibliotecas C e Fortran. SciPy converte a túa sesión interactiva de Python nun ambiente completo de ciencia de datos como MATLAB, IDL, Octave, R ou SciLab.
Neste artigo, analizaremos as técnicas básicas de programación matemática: resolución de problemas de optimización condicional para unha función escalar de varias variables mediante o paquete scipy.optimize. Os algoritmos de optimización sen restricións xa foron discutidos en . Sempre se pode obter axuda máis detallada e actualizada sobre funcións scipy usando o comando help(), Maiús+Tab ou en .
Introdución
A función proporciona unha interface común para resolver problemas de optimización condicionais e non restrinxidos no paquete scipy.optimize. minimize(). Non obstante, sábese que non existe un método universal para resolver todos os problemas, polo que a elección dun método adecuado, como sempre, recae sobre os ombreiros do investigador.
O algoritmo de optimización apropiado especifícase mediante o argumento da función minimize(..., method="").
Para a optimización condicional dunha función de varias variables, están dispoñibles implementacións dos seguintes métodos:
trust-constr— buscar un mínimo local na rexión de confianza. , ;SLSQP— programación cuadrática secuencial con restricións, método newtoniano para resolver o sistema de Lagrange. .TNC- Newton truncado Constrinxido, número limitado de iteracións, bo para funcións non lineais cun gran número de variables independentes. .L-BFGS-B— un método do equipo Broyden–Fletcher–Goldfarb–Shanno, implementado cun consumo de memoria reducido debido á carga parcial de vectores da matriz de Hesse. , .COBYLA— MARE Optimización restrinxida por aproximación lineal, optimización restrinxida con aproximación lineal (sen cálculo de gradientes). .
Dependendo do método elixido, as condicións e restricións para resolver o problema fíxanse de forma diferente:
- obxecto de clase
Boundspara os métodos L-BFGS-B, TNC, SLSQP, trust-constr; - a lista
(min, max)para os mesmos métodos L-BFGS-B, TNC, SLSQP, trust-constr; - un obxecto ou unha lista de obxectos
LinearConstraint,NonlinearConstraintpara COBYLA, SLSQP, métodos trust-constr; - dicionario ou lista de dicionarios
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}para métodos COBYLA, SLSQP.
Esquema do artigo:
1) Considere o uso dun algoritmo de optimización condicional na rexión de confianza (method="trust-constr") con restricións especificadas como obxectos Bounds, LinearConstraint, NonlinearConstraint ;
2) Considere a programación secuencial usando o método dos mínimos cadrados (método = "SLSQP") con restricións especificadas en forma de dicionario {'type', 'fun', 'jac', 'args'};
3) Analiza un exemplo de optimización de produtos fabricados utilizando o exemplo dun estudo web.
Método de optimización condicional = "trust-constr"
Implementación do método trust-constr baseado en para problemas con restricións da forma de igualdade e así para problemas con restricións en forma de desigualdades. Ambos métodos son implementados por algoritmos para atopar un mínimo local na rexión de confianza e son moi axeitados para problemas a gran escala.
Formulación matemática do problema de atopar un mínimo en forma xeral:
Para restricións de igualdade estritas, o límite inferior é igual ao límite superior .
Para unha restrición unidireccional, establécese o límite superior ou inferior np.inf co signo correspondente.
Sexa necesario atopar o mínimo dunha función de Rosenbrock coñecida de dúas variables:
Neste caso, establécense as seguintes restricións no seu dominio de definición:
No noso caso, hai unha solución única no punto , para o que só son válidas as restricións primeira e cuarta.
Imos repasar as restricións de abaixo a arriba e ver como podemos escribilas en scipy.
Restricións и imos definilo usando o obxecto Límites.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restricións и Escribimos en forma lineal:
Definamos estas restricións como un obxecto LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])E finalmente a restrición non lineal en forma matricial:
Definimos a matriz jacobiana para esta restrición e unha combinación lineal da matriz hessiana cun vector arbitrario :
Agora podemos definir unha restrición non lineal como un obxecto 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 o tamaño é grande, as matrices tamén se poden especificar en forma escasa:
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)ou como obxecto 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)Ao calcular a matriz de Hesse require moito esforzo, podes usar unha clase . As seguintes estratexias están dispoñibles: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())O Hessiano tamén se pode calcular usando diferenzas finitas:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')A matriz xacobiana para restricións tamén se pode calcular usando diferenzas finitas. Non obstante, neste caso non se pode calcular a matriz de Hesse usando diferenzas finitas. O Hessian debe definirse como unha función ou mediante a clase HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())A solución ao problema de optimización é a seguinte:
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 é necesario, a función para calcular o Hessian pódese definir usando a 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)ou o produto do Hessiano e un vector arbitrario a través do 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, pódense aproximar as derivadas primeira e segunda da función que se está optimizando. Por exemplo, o hessiano pódese aproximar usando a función SR1 (aproximación cuasi newtoniana). O gradiente pódese aproximar mediante diferenzas 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"
O método SLSQP está deseñado para resolver problemas de minimización dunha función da forma:
Onde и — conxuntos de índices de expresións que describen restricións en forma de igualdades ou desigualdades. — conxuntos de límites inferiores e superiores para o dominio de definición da función.
As restricións lineais e non lineais descríbense en forma de dicionarios con teclas 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])
}A busca do mínimo realízase do seguinte xeito:
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 ]Exemplo de optimización
En relación coa transición á quinta estrutura tecnolóxica, vexamos a optimización da produción co exemplo dun estudo web, que nos aporta uns ingresos pequenos pero estables. Imaxinémonos como o director dunha galera que produce tres tipos de produtos:
- x0 - venda de páxinas de destino, a partir de 10 tr.
- x1 - sitios web corporativos, a partir de 20 tr.
- x2 - tendas en liña, a partir de 30 tr.
O noso amable equipo de traballo está formado por catro xuvenís, dous medios e un sénior. O seu fondo de tempo de traballo mensual:
- Xuño:
4 * 150 = 600 чел * час, - medios:
2 * 150 = 300 чел * час, - señor:
150 чел * час.
Permitir que o primeiro menor dispoñible pase (0, 1, 2) horas no desenvolvemento e implantación dun sitio de tipo (x10, x20, x30), medio - (7, 15, 20), senior - (5, 10, 15). ) horas do mellor momento da túa vida.
Como calquera director normal, queremos maximizar os beneficios mensuais. O primeiro paso para o éxito é escribir a función obxectivo value como a cantidade de ingresos dos produtos producidos ao mes:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Este non é un erro; ao buscar o máximo, a función obxectivo minimízase co signo contrario.
O seguinte paso é prohibir aos nosos empregados o exceso de traballo e introducir restricións no horario de traballo:
O que é 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]])
}Unha restrición formal é que a produción do produto só debe ser positiva:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])E, finalmente, a hipótese máis optimista é que, debido ao baixo prezo e á alta calidade, unha fila de clientes satisfeitos está constantemente facendo cola para nós. Podemos escoller os volumes de produción mensuais nós mesmos, baseándonos en resolver o problema de optimización restrinxido 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 enteiros e calculemos a carga mensual de remeiros cunha distribución óptima de produtos x = (8, 6, 3) :
- Xuño:
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 o director reciba o seu merecido máximo, é óptimo crear 8 landing pages, 6 sitios de tamaño medio e 3 tendas ao mes. Neste caso, o sénior debe arar sen levantar a vista da máquina, a carga dos medios será de aproximadamente 2/3, os xuvenís menos da metade.
Conclusión
O artigo describe as técnicas básicas para traballar co paquete scipy.optimize, usado para resolver problemas de minimización condicional. Persoalmente uso scipy con fins puramente académicos, polo que o exemplo dado é de carácter tan humorístico.
Pódense atopar unha gran cantidade de teoría e exemplos virtuais, por exemplo, no libro de I.L. Akulich "Programación matemática en exemplos e problemas". Aplicación máis hardcore scipy.optimize para construír unha estrutura 3D a partir dun conxunto de imaxes () pódese ver en .
A principal fonte de información é aqueles que desexen contribuír á tradución desta e doutras seccións scipy Benvido a .
Grazas para a participación na elaboración da publicación.
Fonte: www.habr.com