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,592 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 4 dias atrás 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 2 dias atrás 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 2 dias atrás 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 23 horas atrás por Saulo (451 pontos)  
editado 23 horas atrás 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 23 horas atrás 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í.
...