SciPy, otimização com condições

SciPy (pronuncia-se sai pie) é um pacote matemático baseado em numpy que também inclui as bibliotecas C e Fortran. Com o SciPy, uma sessão interativa em Python se transforma em um ambiente completo de processamento de dados, como MATLAB, IDL, Octave, R ou SciLab.

Neste artigo, consideraremos as técnicas básicas de programação matemática - resolvendo problemas de otimização condicional para uma função escalar de várias variáveis ​​usando o pacote scipy.optimize. Algoritmos de otimização incondicional já são considerados no último artigo . Ajuda mais detalhada e atualizada sobre as funções do scipy sempre pode ser obtida usando o comando help (), Shift + Tab, ou na documentação oficial .

1. Introdução

A interface geral para resolver problemas de otimização condicional e incondicional no pacote scipy.optimize é fornecida pela função minimize() . No entanto, sabe-se que não existe um método universal para solucionar todos os problemas, portanto, como sempre, a escolha de um método adequado cabe ao pesquisador.

Um algoritmo de otimização adequado é especificado usando o argumento para a função minimize(..., method="") .

Para otimização condicional da função de várias variáveis, os seguintes métodos estão disponíveis:

• trust-constr - procure um mínimo local na área de confiança. Artigo Wiki, artigo Habr ;
• SLSQP - programação quadrática seqüencial com restrições, o método newtoniano de resolver o sistema de Lagrange. Artigo da Wiki .
• TNC - Newton restrito truncado, um número limitado de iterações, é bom para funções não lineares com um grande número de variáveis ​​independentes. Artigo da Wiki .
• L-BFGS-B - um método dos quatro Broyden - Fletcher - Goldfarb - Shanno, implementado com consumo de memória reduzido devido ao carregamento parcial de vetores da matriz Hessiana. Artigo Wiki , artigo Habr .
• COBYLA - MARES Otimização restrita por aproximação linear, otimização limitada com aproximação linear (sem cálculo de gradiente). Artigo da Wiki .

Dependendo do método escolhido, as condições e restrições para resolver o problema são definidas de forma diferente:

• um objeto da classe Bounds para 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;
• objeto ou lista de objetos LinearConstraint , NonlinearConstraint para métodos COBYLA, SLSQP, trust-constr;
• lista de dicionários ou dicionários {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} para os métodos COBYLA, SLSQP.

O esboço do artigo:
1) Considere a aplicação do algoritmo de otimização condicional na área de confiança (método = "trust-constr") com as restrições especificadas na forma de objetos Bounds , LinearConstraint , NonlinearConstraint ;
2) Considere a programação seqüencial de mínimos quadrados (método = "SLSQP") com as restrições especificadas na forma de um dicionário {'type', 'fun', 'jac', 'args'} ;
3) Desmontar um exemplo de otimização de produtos no exemplo de um web studio.

Método de otimização condicional = "trust-constr"

A implementação do trust-constr baseia-se no EQSQP para problemas com restrições da forma de igualdade e no TRIP 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 área de confiança e são adequados para tarefas de grande escala.

A formulação matemática do problema mínimo de pesquisa na forma geral:

$$

$$\ min_x f (x)$$

$$

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

$$

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

Para restrições estritas de igualdade, o limite inferior é definido igual ao superior $$$c ^ l_j = c ^ u_j$.
Para uma restrição unidirecional, um limite superior ou inferior é definido por np.inf com o sinal correspondente.
Seja necessário encontrar o mínimo da função conhecida de Rosenbrock de duas variáveis:

$$

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

Nesse caso, as seguintes restrições são definidas em seu domínio de definição:

$$

$$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$$

No nosso caso, existe uma solução única no momento $$$[x_0, x_1] = [0,4149, 0,1701]$para o qual apenas as primeira e quarta restrições são válidas.
Vamos examinar as restrições de baixo para cima e considerar como elas podem ser escritas em scipy.
Limitações $$$0 \ leq x_0 \ leq 1$e $$$-0,5 \ leq x_1 \ leq 2,0$definido usando o objeto Bounds.

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

Limitações $$$x_0 + 2 x_1 \ leq 1$e $$$2 x_0 + x_1 = 1$nós escrevemos na forma linear:

$$

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

Defina 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, uma restrição não linear na forma de matriz:

$$

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

Definimos a matriz de Jacobi para essa restrição e a combinação linear da matriz de Hessian com um vetor arbitrário $$$v$:

$$

$$J (x) = \ begin {bmatrix} 2x_0 e 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 e 0 \\ 2 & 0 \ end {bmatrix}$$

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 de 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) 

Quando o cálculo da matriz de Hessian $$$H (x, v)$caro, você pode usar a classe HessianUpdateStrategy . As seguintes estratégias estão disponíveis: BFGS e SR1 .

 from scipy.optimize import BFGS nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS()) 

O Hessian 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 de Jacobi para restrições também pode ser calculada usando diferenças finitas. No entanto, nesse 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 é 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 necessário, a função de cálculo 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 Hessian 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, as primeira e segunda derivadas da função otimizada podem ser calculadas aproximadamente. Por exemplo, o Hessian pode ser aproximado usando a função SR1 (aproximação quasi-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 o problema de minimizar a função na 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$$

Onde $$$\ mathcal {E}$e $$$\ mathcal {I}$- conjuntos de índices de expressão que descrevem restrições na forma de igualdades ou desigualdades. $$$[lb_i, ub_i]$- conjuntos de limites inferior e superior para o domínio de definição de função.

Restrições lineares e não lineares são descritas na forma de dicionários com o type teclas, fun e 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 pesquisa mínima é 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

Em conexão com a transição para o quinto modo tecnológico, consideremos a otimização da produção usando o exemplo de um estúdio na web, o que gera uma renda pequena, mas estável. Vamos nos apresentar como diretor de uma galeria, que produz três tipos de produtos:

• x0 - venda de desembarques, a partir de 10 tr
• x1 - sites corporativos, a partir de 20 tr
• x2 - lojas online, a partir de 30 tr

Nossa equipe de trabalho amigável inclui quatro Juns, dois Middle e um Senior. Fundo do seu tempo de trabalho durante um mês:

• Junho: 4 * 150 = 600 * ,
• meio: 2 * 150 = 300 * ,
• Senior: 150 *

Suponha que o primeiro júnior a gastar no desenvolvimento e implantação de um site do tipo (x0, x1, x2) deva gastar (10, 20, 30) horas, meio - (7, 15, 20), sênior - (5, 10, 15) horas das melhores época de sua vida.

Como qualquer diretor normal, queremos maximizar nosso lucro mensal. O primeiro passo para o sucesso é anotar o value função objetivo como a soma da renda da produção mensal:

 def value(x): return - 10*x[0] - 20*x[1] - 30*x[2] 

Isso não é um erro: ao procurar o máximo, a função objetivo é minimizada com o sinal oposto.

O próximo passo é proibir o processamento para nossos funcionários e introduzir restrições ao fundo de tempo de trabalho:

$$

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

O que é equivalente:

$$

$$\ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix} - \ begin {bmatrix} 10 e 20 e 30 \\ 7 e 15 e 20 \\ 5 e 10 e 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]]) } 

Restrição formal - o resultado deve ser apenas positivo:

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

E, finalmente, a suposição mais otimista - por causa do preço baixo e da alta qualidade, uma linha de clientes satisfeitos está constantemente se alinhando para nós. Nós mesmos podemos escolher os volumes mensais de produção, com base na solução do problema de otimização condicional 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] 

Arredondamos o pirulito para o número inteiro mais próximo e calculamos a carga mensal de remadores com uma distribuição ideal da produção x = (8, 6, 3) :

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

Conclusão: para que o diretor obtenha seu merecido máximo, é ideal realizar 8 desembarques, 6 locais médios e 3 lojas por mês. Ao mesmo tempo, o senhor deve arar sem se afastar da máquina, a carga do meio é de cerca de 2/3, menos da metade de junho.

Conclusão

O artigo descreve as técnicas básicas para trabalhar com o pacote scipy.optimize que são usadas para resolver problemas de minimização condicional. Pessoalmente, eu uso o scipy puramente para fins acadêmicos, então este exemplo é tão cômico.

Muitos exemplos de teoria e winrar podem ser encontrados, por exemplo, no livro de I. L. Akulich "Programação matemática em exemplos e problemas". Uma aplicação mais scipy.optimize do scipy.optimize para criar uma estrutura 3D para um conjunto de imagens ( um artigo no hub ) pode ser encontrada no livro de receitas scipy-cookbook .

A principal fonte de informação é docs.scipy.org.Se você deseja traduzir esta e outras seções do scipy em tradução, bem-vindo ao GitHub .

Agradecemos a mephistopheies por contribuir com esta publicação.