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
4,705 visitas
perguntada Set 21, 2016 em Ciência da Computação por Saulo (451 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

+3 votos
respondida Set 21, 2016 por Saulo (451 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 (231 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.
comentou Set 28 por paulo (1 ponto)  
x1+⋯+x10 = 1

Essa retrição significa oque ? eu nao consegui entender ela .
comentou Set 28 por Caue (231 pontos)  
Xi é o percentual que cada ativo representa na sua carteira. Se X4 é 0.35, então 35% da sua carteira é o ativo 4.
Assim, a soma do percentual de todos os ativos da sua carteira deve dar 100%, ou seja, 1.
comentou Set 28 por paulo (1 ponto)  
O valor maximo é : 0.00299025974026
não consegui entender oque ele significa neste exemplo, alguem poderia me explica, e também o resultado das variaveis .
comentou Set 29 por Caue (231 pontos)  
“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.”
Se o problema está maximizando o retorno esperado da carteira, então o máximo calculado é exatamente o retorno esperado da carteira.
Sugiro ler tudo novamente com atenção, pois está bem explicado na resposta do Saulo.
comentou Nov 7 por paulo (1 ponto)  
Nesse exemplo, na hora que você faz o modelo e passa para o algoritmo você faz a divisão por 100. teria como explicar o porque dessa divisão? se eu fizer um algoritmo para calcular o simplex do 0, na hora que passar como parâmetro para o calculo eu preciso de passar os valores divididos por 100? E desde já, parabéns pelo post
comentou Nov 7 por Saulo (451 pontos)  
Obrigado pelo retorno.

Porcentagem, por definição, significa "por cento" ou "por cada centena", que, por sua vez, significa dividir o número por 100. Ou seja, 10% na verdade é 10 em 100 partes, 10/100 = 0,1. Você não precisa necessariamente fazer essa divisão por 100. Você consegue reescrever o problema para considerar as porcentagens de 0 a 100 (ao invés de 0 a 1, como eu fiz), e seu ótimo será dado por (x1, x2, ..., x9, x10) = (30, 26.67, 0, 13.37, 0, 0, 30, 0, 0, 0).

Não entendi sua pergunta sobre calcular o simplex do 0.
comentou Nov 7 por paulo (1 ponto)  
Quando eu digo calcular do 0, falo que se eu fizer o algoritmo simplex, sem usar biblioteca alguma, quando eu passar os valores para o simplex, eu teria que passar os valores divididos por 100? ou eu poderia passar os valores diretos?
comentou Nov 7 por Saulo (451 pontos)  
A divisão por 100 é a forma como você estrutura matematicamente o seu problema. Então você não precisa dividir por 100, você poderia passar diretamente os valores, mas você tem que ser coerente na forma como você define cada equação. Se você fizer isso, vai obter o mesmo resultado. Depois tente e veja que dá certo! Se você não conseguir, me avise.
comentou Nov 9 por paulo (1 ponto)  
import numpy as np
from scipy.optimize import linprog

c = [0.37, 0.24, 0.14, 0.30, 0.24, 0.19, 0.28, 0.18, 0.25, 0.24]
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(30*np.ones((10,)), np.zeros((10,)), axis=0)
A_ub = np.append(A_ub, [[2.48, 2.16, 1.95, 2.93, 2.40, 2.00, 2.63, 2.14, 2.73, 2.47]], axis=0)
b_ub = np.append(b_ub, [2.50])

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



O algoritmo ficaria assim? eu não consegui esse mesmo resultado.
comentou Nov 9 por paulo (1 ponto)  
Uma outra duvida também, e se eu quisesse trabalhar com 5 variáveis ou outras condições, oque eu preciso mudar no código, pergunto isso porque eu não manjo muito dessa área de programação e eu tenho que entregar um código na faculdade. Desde já agradeço pelas respostas, você ta me ajudando muito !!!
comentou Nov 11 por Saulo (451 pontos)  
editado Nov 11 por Saulo
Como eu já comentei, você tem que ser coerente usando variáveis percentuais, se você está considerando o percentual de 0 a 100, ou de 0 a 1. Como você está considerando de 0 a 100, sua soma tem que dar 100 (e não 1) na equação 1 da restrição. Na equação 3 da restrição, você está fazendo uma proporção; como x1, ..., x10 variam de 0 a 100, você tem que dividir cada uma dessa variáveis por 100 (ou seja, os coeficientes dessa linha da matriz A tem que levar em consideração essa divisão por 100). Se você fizer essa alteração no seu código, vai dar certo.
comentou Nov 11 por Saulo (451 pontos)  
Você consegue fazer o que você quiser, desde que seja coerente com a estrutura do seu problema matemático e o seu código. Sugiro você estudar melhor o código (ou implementar o seu do zero) para entender melhor. Isso vai te ajudar no futuro, na faculdade e, principalmente, em sua vida profissional. Como sugestão, tente analisar melhor a parte matemática, e veja onde isso está presente no código. Tente fazer as alterações a partir daí.
...