Modelando um problema de Programação Linear e resolvendo com Python e Pyomo
De forma resumida, um problema de otimização é quando desejamos encontrar um conjunto de valores, dentro de um espaço restrito, com o objetivo de minimizar ou maximizar uma função. Quando o problema de otimização possui função objetivo e restrições lineares, temos um problema de Programação Linear.
A Programação Linear possui inúmeras aplicações nas mais variadas áreas, tais como Marketing, Supply Chain, Saúde, Logística, Segurança Pública, entre outras. Neste artigo o objetivo é, através de um case didático, demonstrar o processo de modelagem e a resolução com Python e Pyomo.
O Pyomo é um pacote baseado em Python para lidar com problemas de otimização. Um modelo escrito em Pyomo pode ser resolvido através de diversos solvers, entre eles CPLEX, Gurobi e GLPK. entretanto, fica por conta do pacote realizar as conversões para os formatos específicos de cada solver, dispensando o usuário de realizar alterações significativas no código.
Para quem tem pouca ou nenhuma experiência com Python e Pyomo, recomendo este manual, desenvolvido pelo Claudemir Woche e o professor Anselmo Pitombeira, da Universidade Federal do Ceará. Talvez este seja o melhor material sobre Pyomo em língua portuguesa.
Para baixar o notebook com todos os códigos, acesse o repositório no GitHub.
O problema do mix de produção
Vamos utilizar o exemplo da página 26 do livro Pesquisa Operacional para Cursos de Engenharia (Belfiore e Fávero, 2013).
A empresa Venix de brinquedos está revendo seu planejamento de produção de carrinhos e triciclos. O lucro líquido por unidade de carrinho e triciclo produzido é de R$ 12,00 e R$ 60,00, respectivamente.
As matérias-primas e os insumos necessários para a fabricação de cada um dos produtos são terceirizados, cabendo à empresa os processos de usinagem, pintura e montagem. O processo de usinagem requer 15 minutos de mão de obra especializada por unidade de carrinho e 30 minutos por unidade de triciclo produzida.
O processo de pintura requer 6 minutos de mão de obra especializada por unidade de carrinho e 45 minutos por unidade de triciclo produzida.
Já o processo de montagem necessita de 6 minutos e 24 minutos para uma unidade de carrinho e de triciclo produzida, respectivamente.
O tempo disponível por semana é de 36, 22 e 15 horas para os processos de usinagem, pintura e montagem, respectivamente.
A empresa quer determinar quanto produzir de cada produto por semana, respeitando as limitações de recursos, de forma a maximizar o lucro líquido semanal.
Formulação do modelo matemático
Antes de partir para o código, é importante ter o modelo matemático bem definido.
Variáveis de decisão
Já sabemos o quanto a empresa lucra por unidades vendidas de carrinhos (R$ 12,00) e triciclos (R$ 60,00), agora precisamos saber o quanto produzir de cada um deles, para que o lucro semanal seja o maior possível. Se precisamos decidir o quanto produzir de cada produto, as quantidades a serem produzidas são, de fato, a resposta para o nosso problema. Portanto, temos aqui as variáveis de decisão do modelo:
Função objetivo
O objetivo do problema é maximizar o lucro líquido semanal da empresa, e para que isto aconteça, precisamos decidir o quanto fabricar de carrinhos (𝑥1) e triciclos (𝑥2). O lucro semanal, portanto, é calculado da seguinte forma:
A função objetivo (𝑧) é quem relaciona o objetivo do problema com as variáveis de decisão. Portanto, fica definida por:
A notação abaixo, no entanto, é mais comum:
Restrições
É necessário "informar" para o modelo todas as condições que o problema impõe, sejam limitações de recursos, requisitos de negócio e etc. Portanto, é aqui que incluímos as restrições do modelo.
Neste problema consideraremos que as matérias-primas e os insumos necessários para a fabricação são ilimitados, portanto, são informações desprezadas pelo modelo. No entanto, existem limitações de recursos para os processos de usinagem, pintura e montagem. A tabela abaixo resume o tempo de mão de obra especializada necessário em cada processo, por produto:
Restrições de disponibilidade de mão de obra
Para o processo de usinagem, o tempo disponível é de 36 horas. Ou seja, dadas as limitações de recursos (máquinas, equipamentos, colaboradores e etc.), o processo de usinagem deve consumir no máximo 36 horas. Formulamos esta restrição da seguinte forma:
Quanto ao processo de pintura, o tempo disponível é de 22 horas. Logo:
Por fim, para o processo de pintura o tempo disponível é de 15 horas. Portanto, temos:
Restrições de não negatividade
Modelamos as restrições de não negatividade para que o modelo não encontre valores negativos para as variáveis de decisão, uma vez que seria impossível fabricar quantidades negativas de carrinhos e triciclos. Portanto, são elas:
O modelo matemático completo fica definido abaixo:
Representação matricial
O case em estudo possui 2 variáveis de decisão e 5 restrições. Você deve estar se perguntando: e se eu tiver um modelo com 200 variáveis de decisão e 1000 restrições, vou ter que escrever tudo isso? A resposta é não! Existe uma forma genérica de representar o modelo, mas para chegar lá é importante compreender o modelo na forma matricial.
Importante: se você não tem familiaridade com operações com matrizes, recomendo fortemente que assista estas ótimas aulas do professor Fernando Grings ou visite este ótimo material do professor Wagner Bonat.
Função objetivo
Para expressarmos a função objetivo na forma matricial, vamos primeiro atribuir as variáveis de decisão a um vetor coluna, que nada mais é que uma matriz 𝑚×1 (𝑚 linhas e 1 coluna):
Em seguida, atribuimos os coeficientes da função objetivo a um vetor linha (transposta do vetor coluna):
Portanto, dadas as características da multiplicação de matrizes, a função objetivo fica definida por:
Restrições
Para representarmos matricialmente as restrições de disponibilidade de mão de obra, começaremos atribuindo os coeficientes das restrições a uma matriz 𝑚×𝑛, em que 𝑚 é o número de processos de fabricação (linhas), no caso 3 (usinagem, pintura e montagem) e 𝑛 é o número de produtos (colunas), que no caso são 2 (carrinhos e triciclos):
Agora, atribuimos o lado direito das inequações (RHS), que são as horas disponíveis para execução de cada um dos processos, a um vetor coluna:
Logo, as restrições de disponibilidade de mão de obra ficam definidas por:
E, por fim, definimos matricialmente as restrições de não negatividade:
O modelo completo na forma matricial fica da seguinte forma:
Representação genérica
Agora que conseguimos enxergar matricialmente o modelo matemático, fica muito mais fácil de compreendermos a sua forma genérica. A forma genérica, além de tudo, vai tornar o processo de “tradução” para o modelo computacional muito mais fácil.
Índices
Pensando no modelo matricial, precisamos representar os índices das linhas e colunas das matrizes. Para facilitar, vamos começar analisando a matriz 𝐀, que é a mais completa do modelo.
A representação genérica da matriz 𝐀 fica da seguinte forma:
Sendo 𝑎𝑖𝑗 o elemento da linha 𝑖 e coluna 𝑗 da matriz.
Observe que o índice das linhas varia de 1 a 3, podendo ser representado como 𝑖=1,2,3. Seguindo a mesma ideia, o índice das colunas pode ser representado como 𝑗=1,2. Generalizando, temos:
Note que no problema em estudo, 𝑚 é o número de linhas da matriz, ou seja, o número de processos de fabricação. Portanto, 𝑖=1,2,3 representa os processos de usinagem, pintura e montagem, respectivamente.
Seguindo na mesma linha de raciocínio, 𝑛 é o número de colunas da matriz, ou seja, o número de produtos que a empresa fabrica. Portanto, 𝑗=1,2 representa os produtos carrinho e triciclo, respectivamente.
Ainda olhando para as restrições, temos o vetor coluna 𝐛 representado de forma genérica:
Observe que cada linha corresponde a um processo de fabricação, representados pelo índice 𝑖=1,2,3.
Agora olhando para a função objetivo, vamos analisar o vetor linha 𝐜𝑇:
Observe que cada elemento do vetor corresponde a um produto dos produtos fabricados, que são representados pelo índice 𝑗=1,2. O mesmo raciocínio serve para o vetor coluna 𝐱.
Portanto, para o problema em estudo, os índices 𝑖 e 𝑗 são suficientes.
Tradicionalmente o índice 𝑖 está associado às linhas da matriz, e o 𝑗, às colunas. Outro ponto de atenção é que normalmente 𝑚 está para 𝑖, assim como 𝑛 está para 𝑗. Existem variações de uma fonte para outra, no entanto, é importante que estas relações estejam bem definidas no nosso modelo, para que a gente não faça confusão.
Parâmetros
Agora que conhecemos os índices do modelo, podemos declarar os parâmetros, que são todos os dados de entrada do modelo. Ou seja, dados que já conhecemos previamente. Não confunda com variáveis de decisão!
Começando pelos coeficientes da função objetivo, temos:
Agora olhando para as restrições, temos:
Variáveis de decisão
A representação genérica das variáveis de decisão fica da seguinte forma:
Função objetivo
A representação genérica da função objetivo expressa a soma dos produtos entre os parâmetros 𝑐𝑗 e as variáveis de decisão 𝑥𝑗, ou seja:
Portanto, temos:
Restrições
A representação genérica das restrições de disponibilidade de mão de obra expressa a inequação composta pela soma dos produtos entre os parâmetros 𝑎𝑖𝑗 e as variáveis de decisão 𝑥𝑗, e o parâmetro 𝑏𝑖. Ou seja:
Logo, temos:
Por fim, as restrições de não negatividade:
A formulação genérica completa fica da seguinte forma:
Modelo computacional e resolução
Agora que temos bem claro o modelo genérico, o trabalho de “tradução” para o Pyomo fica muito mais fácil.
Vamos começar definindo os dados de entrada do modelo. Como estamos lidando com um exemplo meramente didático, vou produzir estes dados direto no código. Mas aqui poderiamos ler estes dados de alguma fonte externa.
Lembre-se que se estamos falando de dados de entrada, estamos falando dos parâmetros do modelo, que são eles:
c = [12,60]a = [[0.25, 0.5],
[0.1, 0.75],
[0.1, 0.4]]b = [36, 22, 15]
Agora vamos importar o Pyomo:
import pyomo.environ as pyoA próxima etapa é declarar o modelo:
modelo = pyo.ConcreteModel()Agora vamos declarar os índices no modelo. Sabendo que a função len(a) retorna o número de linhas do parâmetro 𝑎𝑖𝑗, e a função len(a[0][:]), o número de colunas, podemos declarar os índices da seguinte forma:
modelo.m = range(len(a))modelo.n = range(len(a[0][:]))
O próximo passo é declarar os parâmetros no modelo. Nesta etapa é importante ter bastante atenção com os índices!
modelo.c = pyo.Param(modelo.n, initialize=lambda modelo, j: c[j])modelo.a = pyo.Param(modelo.m, modelo.n, initialize=lambda modelo, i, j: a[i][j])modelo.b = pyo.Param(modelo.m, initialize=lambda modelo, i: b[i])
Agora declaramos as variáveis de decisão. Atenção com o índice utilizado!
modelo.x = pyo.Var(modelo.n, within=pyo.NonNegativeReals)Chegou a vez da função objetivo. Observe que a função do Pyomo summation() facilita a nossa vida, já que faz o trabalho de construir a soma dos produtos:
def f_obj(modelo):
return pyo.summation(modelo.c, modelo.x)modelo.z = pyo.Objective(rule=f_obj, sense=pyo.maximize)
E, por fim, as restrições. Observe que na declaração das variáveis de decisão já foram incluídas as restrições de não negatividade (within=pyo.NonNegativeReals). Portanto, faltam apenas as retrições de disponibilidade de mão de obra:
def f_constr(modelo, i):
return sum(modelo.a[i,j] * modelo.x[j] for j in modelo.n) <= modelo.b[i]modelo.restr_horas = pyo.Constraint(modelo.m, rule=f_constr)
Observe uma relação importante: Dentro da função sum(), o trecho for j in modelo.n representa o índice do somatório, ao passo que dentro da função Constraint(), o argumento modelo.m representa o ∀𝑖 (para todo 𝑖).
Agora que temos o modelo computacional completamente declarado, vamos para a resolução:
resultado = pyo.SolverFactory('gurobi', solver_io="python").solve(modelo)Para visualizar o valor das variáveis de decisão:
l = list(modelo.x.keys())
for i in l:
print("x" + str(i+1),'=', modelo.x[i]())
E para visualizar o valor da função objetivo:
modelo.z()
Portanto, interpretamos a saída do modelo da seguinte forma: para alcançar o lucro máximo semanal de R$ 2.040,00, a empresa deve fabricar 70 unidades do produto 1 (carrinho) e 20 unidades do produto 2 (triciclo).
Observe que isto não é uma previsão. A saída do modelo é uma prescrição do que a empresa deve fazer para atingir o objetivo de maximizar o lucro. Portanto, sob a perspectiva da ciência de dados, a otimização está dentro do contexto da análise prescritiva.
Conclusões
O objetivo deste artigo foi demonstrar o processo de modelagem e a resolução de um case didático de Programação Linear, com o uso do Pyomo, uma biblioteca baseada em Python para lidar com problemas de otimização. A principal mensagem aqui é que a plena compreensão do problema e a formulação correta do modelo matemático são as principais etapas do processo. Se houverem falhas nestas etapas, a ferramenta utilizada para resolver o problema, seja o Pyomo ou qualquer outra, não irá fazer milagres, já que garbage in, garbage out.
Também foi demonstado que a boa compreensão da forma genérica do modelo, com índices e parâmetros, torna a construção do modelo computacional muito mas fácil. Mas para compreender integralmente a forma genérica, é importante compreender a forma matricial, que por sua vez exige conhecimento em álgebra linear, principalmente no contexo de matrizes.
Por fim, mesmo o Pyomo não estando entre as mais populares linguagens de modelagem de otimização, o fato do Python ser uma das linguagens de programação mais utilizadas no mundo torna o seu uso bastante atrativo.
Referências
BELFIORE, P; FÁVERO, L. P. Pesquisa Operacional Para Cursos de Engenharia. [s.l.] Elsevier Brasil, 2013.
BONAT, W. H. Tópicos em Matemática para Cientistas de Dados. 2021. Disponível em: http://leg.ufpr.br/~wagner/TMCD/. Acesso em: 24 de abr. 2021.
CARVALHO, C. W. V.; PITOMBEIRA NETO, A. R. Manual de uso da biblioteca Pyomo para Programação Matemática. 2020. Disponível em: http://www.opl.ufc.br/files/Manual_Pyomo_2020.pdf. Acesso em: 12 de jan. 2021.
GRINGS, F. Matriz transposta, matriz inversa, determinantes e sistemas lineares, regra de cramer, escalonamento. 2020. Disponível em: https://www.youtube.com/playlist?list=PLE6qFDd4x9w9mmuAHJuBspzGI9mQ_sNVB. Acesso em: 24 de abr. 2021.
SCARPIN, C. T. Programação Linear PPGMNE. 2020. Disponível em: https://www.youtube.com/watch?v=Zspg7bl4unU&list=PLY8sfsQRSw792jf1y1KDvRYbL5QMB5CoG. Acesso em: 24 de abr. 2021.
