Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como rodar o algoritmo de SIMPLEX para resolver um problema de programação linear em Python?

+1 voto
2,043 visitas
perguntada Set 21, 2016 em Ciência da Computação por Saulo (426 pontos)  

Dê um exemplo interessante de como rodar o algoritmo de SIMPLEX para resolver um problema de programação linear em Python.

Compartilhe

1 Resposta

+2 votos
respondida Set 21, 2016 por Saulo (426 pontos)  
editado Set 21, 2016 por Saulo

Existem diversas formas para se utilizar o método Simplex em Python:

(1) Implementar o algoritmo. Por exemplo, veja os posts de:

(2) Usar bibliotecas que possuem o algoritmo já implementado, tais como:

  • SciPy.
    Para a minimização de problemas sem restrição, pode-se usar a função fmin, que usa o algoritmo simplex "downhill", também chamado de método de Nelder-Mead, ou a função minimize com o parâmetro '' method='nelder-mead' ''. Para a minimização de problemas com restrição, pode-se usar a função linprog.

  • PuLP.
    A documentação da biblioteca, Optimisation Concepts - Linear Programing, nos mostra que é possível utilizar os algoritmos "Primal Simplex Method", "Dual Simplex Method" e "Interior Point Method" usando a função cplex. Veja o exemplo "A Blending Problem".

Alguns exemplos interessantes de problemas que podem ser resolvidos usando o algoritmo SIMPLEX podem ser obtidos em:

  • BEILBY, M. H. "Economics and Operational Research". Academic Press, 1976.
  • ARENALES, M.; ARMENTANO, V.; MORABITO, R.; YANASSE, H. "Pesquisa Operacional". Elsevier, 2007.
  • BELFIORE, P.; FÁVERO, L. P. "Pesquisa Operacional para Cursos de Engenharia". Elsevier, 2013.

As referências acima também explicam o método Simplex. Vamos utilizar um problema proposto por Belfiore e Fávero (2013), que está detalhado nas páginas 42 a 45 do livro.

Um investidor, que opera diariamente no home broker, deseja selecionar uma nova carteira de investimentos que maximiza seu retorno esperado com um nível de risco assumido. A partir de uma análise histórica das principais ações do Ibovespa, foi selecionado um conjunto de 10 ações que poderiam compor a carteira. O quadro a seguir mostra o retorno médio e o desvio padrão obtido a partir do histórico do retorno diário de cada uma das ações no período 14/1/2009 a 13/1/2010:

\[ \small \begin{array}{ccllcc} ~ & \text{Variável} & \text{Ação} & \text{Código} & \text{Retorno médio} & \text{Desvio-padrão} \\ 1 & x_{1} & \text{Banco Brasil ON} & \text{BBAS3} & 0,37\% & 2,48\% \\ 2 & x_{2} & \text{Bradesco PN} & \text{BBDC4} & 0,24\% & 2,16\% \\ 3 & x_{3} & \text{Eletrobrás PNB} & \text{ELET6} & 0,14\% & 1,95\% \\ 4 & x_{4} & \text{Gerdau PN} & \text{GGBR4} & 0,30\% & 2,93\% \\ 5 & x_{5} & \text{Itausa PN} & \text{ITSA4} & 0,24\% & 2,40\% \\ 6 & x_{6} & \text{Petrobras PN} & \text{PETR4} & 0,19\% & 2,00\% \\ 7 & x_{7} & \text{Sid Nacional ON} & \text{CSNA3} & 0,28\% & 2,63\% \\ 8 & x_{8} & \text{Telemar PN} & \text{TNLP4} & 0,18\% & 2,14\% \\ 9 & x_{9} & \text{Usiminas PNA} & \text{USIM5} & 0,25\% & 2,73\% \\ 10 & x_{10} & \text{Vale PNA} & \text{VALE5} & 0,24\% & 2,47\% \\ \end{array} \]

A fim de diversificar a carteira e, consequentemente, minimizar o risco do portfólio, deseja-se investir, no máximo, 30% em cada uma das ações. Além disso, o risco da carteira, medido por meio do desvio-padrão, não poderia ultrapassar o valor de 2,5%.

Sendo \(x_i\) a porcentagem da ação \(i\), \(i=1,\ldots,10\), o problema em questão pode ser estruturado como:

\[ \small \begin{array}{rl} \displaystyle\max_{x_1,\ldots,x_{10}} & 0,0037\,x_1 + 0,0024\,x_2 + 0,0014\,x_3 + 0,0030\,x_4 + 0,0024\,x_5 \\ ~ & + 0,0019\,x_6 + 0,0028\,x_7 + 0,0018\,x_8 + 0,0025\,x_9 + 0,0024\,x_{10} \\ \text{s.a.} & ~ \\ ~ & \begin{array}{rcl} x_1 + \cdots + x_{10} & = & 1 \\ x_1, \ldots, x_{10} & \leq & 0,30 \\ 0,0248\,x_1 + 0,0216\,x_2 + \cdots + 0,0273\,x_9 + 0.0247\,x_{10} & \leq & 0,0250 \\ x_1, \ldots, x_{10} & \geq & 0 \\ \end{array} \end{array} \]

Para resolver o problema, vamos utilizar a função linprog da biblioteca SciPy, que resolve o seguinte problema com restrição:

\[ \begin{array}{l} \displaystyle\min_{x}{~~~~c^T\,x} \\ \begin{array}{rcl} \text{s.a.} ~ A_{ub} & \leq & b_{ub} \\ A_{eq} & = & b_{eq} \\ \end{array} \end{array} \]

O código a seguir, escrito em Python 2.7, nos fornece o resultado desejado:

import numpy as np
from scipy.optimize import linprog

c = [0.0037, 0.0024, 0.0014, 0.0030, 0.0024, 0.0019, 0.0028, 0.0018, 0.0025, 0.0024]
c = np.multiply(-1.0, c)

A_eq = np.ones((1,10))
b_eq = np.array([1.0])

A_ub = np.append(np.eye(10), -1.0*np.eye(10), axis=0)
b_ub = np.append(0.30*np.ones((10,)), np.zeros((10,)), axis=0)
A_ub = np.append(A_ub, [[0.0248, 0.0216, 0.0195, 0.0293, 0.0240, 0.0200, 0.0263, 0.0214, 0.0273, 0.0247]], axis=0)
b_ub = np.append(b_ub, [0.0250])

res = linprog(c, A_eq=A_eq, b_eq=b_eq, A_ub=A_ub, b_ub=b_ub, method='simplex', options = {'disp': True, 'maxiter': 1000})

print "Valor otimo: ", -res.fun
print "X:"
for k, xk in enumerate(res.x):
    print "x_{", str(k+1), "} = ", xk

A simulação nos fornece que o valor máximo é igual a \(0.00299025974026\) e também que:
\[ \begin{array}{rcl} x_{ 1 } & = & 0.3 \\ x_{ 2 } & = & 0.266233766234 \\ x_{ 3 } & = & 0.0 \\ x_{ 4 } & = & 0.133766233766 \\ x_{ 5 } & = & 0.0 \\ x_{ 6 } & = & 0.0 \\ x_{ 7 } & = & 0.3 \\ x_{ 8 } & = & 0.0 \\ x_{ 9 } & = & 0.0 \\ x_{ 10 } & = & 0.0 \\ \end{array} \]

comentou Set 23, 2016 por Caue (226 pontos)  
Exemplo muito interessante e relacionado ao curso.

É válido destacar que apenas 4 ativos (de 10) foram selecionados pelo modelo proposto. Isso é ruim do ponto de vista de diversificação da carteira. Um dos motivos é que não foi considerada a correlação entre esses ativos, mas apenas o desvio padrão de cada um deles.
...