SciPy (pronuncia-se sai pie) é um pacote matemático baseado em numpy que também inclui bibliotecas C e Fortran. SciPy transforma sua sessão interativa de Python em um ambiente completo de ciência de dados como MATLAB, IDL, Octave, R ou SciLab.
Neste artigo, veremos as técnicas básicas de programação matemática - resolvendo problemas de otimização condicional para uma função escalar de diversas variáveis usando o pacote scipy.optimize. Algoritmos de otimização irrestrita já foram discutidos em . Ajuda mais detalhada e atualizada sobre funções scipy sempre pode ser obtida usando o comando help(), Shift+Tab ou em .
Introdução
Uma interface comum para resolver problemas de otimização condicional e irrestrita no pacote scipy.optimize é fornecida pela função minimize(). Porém, sabe-se que não existe um método universal para resolver todos os problemas, portanto a escolha de um método adequado, como sempre, recai sobre os ombros do pesquisador.
O algoritmo de otimização apropriado é especificado usando o argumento da função minimize(..., method="").
Para otimização condicional de uma função de diversas variáveis, estão disponíveis implementações dos seguintes métodos:
trust-constr— procurar um mínimo local na região de confiança. , ;SLSQP— programação quadrática sequencial com restrições, método newtoniano para resolução do sistema de Lagrange. .TNC- Newton truncado Número restrito e limitado de iterações, bom para funções não lineares com um grande número de variáveis independentes. .L-BFGS-B— um método da equipe Broyden–Fletcher–Goldfarb–Shanno, implementado com consumo de memória reduzido devido ao carregamento parcial de vetores da matriz Hessiana. , .COBYLA— Otimização restrita MARE por aproximação linear, otimização restrita com aproximação linear (sem cálculo de gradiente). .
Dependendo do método escolhido, as condições e restrições para resolver o problema são definidas de forma diferente:
- objeto de classe
Boundspara 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; - um objeto ou uma lista de objetos
LinearConstraint,NonlinearConstraintpara métodos COBYLA, SLSQP, trust-constr; - dicionário ou lista de dicionários
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}para métodos COBYLA, SLSQP.
Esboço do artigo:
1) Considere o uso de um algoritmo de otimização condicional na região de confiança (método=”trust-constr”) com restrições especificadas como objetos Bounds, LinearConstraint, NonlinearConstraint ;
2) Considere a programação sequencial usando o método dos mínimos quadrados (método = "SLSQP") com restrições especificadas na forma de um dicionário {'type', 'fun', 'jac', 'args'};
3) Analisar um exemplo de otimização de produtos manufaturados utilizando o exemplo de um web studio.
Método de otimização condicional = "trust-constr"
Implementação do método trust-constr baseado em para problemas com restrições da forma de igualdade e sobre para problemas com restrições na forma de desigualdades. Ambos os métodos são implementados por algoritmos para encontrar um mínimo local na região de confiança e são adequados para problemas de grande escala.
Formulação matemática do problema de encontrar um mínimo na forma geral:
Para restrições de igualdade estritas, o limite inferior é igual ao limite superior .
Para uma restrição unidirecional, o limite superior ou inferior é definido np.inf com o sinal correspondente.
Seja necessário encontrar o mínimo de uma função de Rosenbrock conhecida de duas variáveis:
Neste caso, são impostas as seguintes restrições ao seu domínio de definição:
No nosso caso, existe uma solução única no ponto , para os quais apenas a primeira e a quarta restrições são válidas.
Vamos examinar as restrições de baixo para cima e ver como podemos escrevê-las em scipy.
Restrições и vamos defini-lo usando o objeto Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restrições и Vamos escrever na forma linear:
Vamos definir essas restrições como um objeto 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 restrição não linear em forma matricial:
Definimos a matriz Jacobiana para esta restrição e uma combinação linear da matriz Hessiana com um vetor arbitrário :
Agora podemos definir uma restrição não linear como um 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)Se o tamanho for grande, as matrizes também poderão ser especificadas na forma esparsa:
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 um 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)Ao calcular a matriz Hessiana requer muito esforço, você pode usar uma classe . As seguintes estratégias estão disponíveis: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())O Hessiano também pode ser calculado usando diferenças finitas:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')A matriz Jacobiana para restrições também pode ser calculada usando diferenças finitas. Contudo, neste caso a matriz Hessiana não pode ser calculada usando diferenças finitas. O Hessian deve ser definido como uma função ou usando a classe HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())A solução para o problema de otimização é assim:
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 necessário, a função de cálculo do Hessian pode ser definida usando a 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)ou o produto do Hessiano e um vetor arbitrário atravé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, a primeira e a segunda derivadas da função que está sendo otimizada podem ser aproximadas. Por exemplo, o Hessiano pode ser aproximado usando a função SR1 (aproximação quase newtoniana). O gradiente pode ser aproximado por diferenças 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 otimização condicional = "SLSQP"
O método SLSQP foi projetado para resolver problemas de minimização de uma função na forma:
Onde и — conjuntos de índices de expressões que descrevem restrições na forma de igualdades ou desigualdades. — conjuntos de limites inferiores e superiores para o domínio de definição da função.
Restrições lineares e não lineares são descritas na forma de dicionários com chaves 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 pelo mínimo é realizada da seguinte forma:
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 otimização
No que diz respeito à transição para a quinta estrutura tecnológica, vejamos a otimização da produção a partir do exemplo de um web studio, que nos traz um rendimento pequeno mas estável. Imaginemo-nos como diretores de uma cozinha que produz três tipos de produtos:
- x0 - venda de landing pages, a partir de 10 tr.
- x1 - sites corporativos, a partir de 20 tr.
- x2 - lojas online, a partir de 30 tr.
Nossa simpática equipe de trabalho inclui quatro juniores, dois intermediários e um sênior. Seu fundo mensal de tempo de trabalho:
- Junho:
4 * 150 = 600 чел * час, - meio:
2 * 150 = 300 чел * час, - senhor:
150 чел * час.
Deixe o primeiro júnior disponível gastar (0, 1, 2) horas no desenvolvimento e implantação de um site do tipo (x10, x20, x30), médio - (7, 15, 20), sênior - (5, 10, 15 ) horas do melhor momento da sua vida.
Como qualquer diretor normal, queremos maximizar os lucros mensais. O primeiro passo para o sucesso é escrever a função objetivo value como o valor da receita dos produtos produzidos por mês:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Isto não é um erro; ao procurar o máximo, a função objetivo é minimizada com o sinal oposto.
O próximo passo é proibir nossos funcionários de trabalhar demais e introduzir restrições ao horário de trabalho:
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]])
}Uma restrição formal é que a saída do produto deve ser apenas positiva:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])E, finalmente, a suposição mais otimista é que, devido ao baixo preço e à alta qualidade, uma fila de clientes satisfeitos se forma constantemente para nós. Podemos escolher nós mesmos os volumes de produção mensais, com base na resolução do problema de otimização restrita com 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]Vamos arredondar livremente para números inteiros e calcular a carga mensal dos remadores com uma distribuição ideal de produtos x = (8, 6, 3) :
- Junho:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - meio:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senhor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Conclusão: para que o diretor receba o máximo merecido, o ideal é criar 8 landing pages, 6 sites de médio porte e 3 lojas por mês. Neste caso, o mais velho deve arar sem tirar os olhos da máquina, a carga dos médios será de aproximadamente 2/3, os mais novos menos da metade.
Conclusão
O artigo descreve as técnicas básicas para trabalhar com o pacote scipy.optimize, usado para resolver problemas de minimização condicional. Pessoalmente eu uso scipy puramente para fins acadêmicos, razão pela qual o exemplo dado é de natureza tão cômica.
Muita teoria e exemplos virtuais podem ser encontrados, por exemplo, no livro de I. L. Akulich “Programação matemática em exemplos e problemas”. Aplicação mais hardcore scipy.optimize construir uma estrutura 3D a partir de um conjunto de imagens () pode ser visualizado em .
A principal fonte de informação é aqueles que desejam contribuir para a tradução desta e de outras seções scipy Bem-vindo ao .
Obrigado pela participação na preparação da publicação.
Fonte: habr.com