Resolvendo o problema das p-Medianas com Python e Pyomo

O Pyomo é um pacote baseado em Python para formulação, resolução e análise de modelos 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á, e aproveito, inclusive, para agradecê-los por disponibilizarem este ótimo material, talvez o melhor que temos em língua portuguesa.

Para demonstrar um pouco sobre o funcionamento do Pyomo, neste artigo será abordada a modelagem e a resolução do problema das p-Medianas com algumas variações.

Para baixar o notebook com os códigos e os conjuntos de dados utilizados, acesse o repositório no GitHub.

Exemplo 1

O problema clássico das p-Medianas tem como objetivo definir a localização de p instalações (medianas) em uma rede de n nós, minimizando a soma das distâncias entre cada nó e a instalação mais próxima. Neste exemplo consideraremos o seguinte caso: uma rede de lojas de departamento decidiu construir 3 novas instalações em uma determinada cidade. Entre os 10 bairros disponíveis, onde construir as novas lojas, de modo que a soma das distâncias entre cada bairro e a loja mais próxima seja mínima?

Aproveitaremos o conjunto de dados utilizado pelo professor Gustavo Loch, da Universidade Federal do Paraná, nas aulas que ele abordou o problema das p-Medianas com C# e Gurobi.

Iniciamos carregando as bibliotecas e o conjunto de dados:

Imagine que cada linha do conjunto de dados representa um dos bairros candidatos a receberem as novas lojas, sendo X e Y as coordenadas do centroide de cada bairro. Se o objetivo é minimizar as distâncias entre cada bairro e a loja mais próxima, vamos precisar, obviamente, de uma matriz de distâncias. Para o exemplo optei por considerar distância euclidiana:

Visualizando a dispersão dos bairros:

O modelo

Antes de partirmos para o Pyomo é importante definirmos o modelo matemático exato.

Uma vez que a matriz é composta por m linhas e n colunas, definimos os índices da seguinte forma:

Em seguida precisamos estabelecer os parâmetros do modelo, ou seja, os dados já conhecidos previamente. No exemplo em questão os parâmetros são a matriz de distâncias e o número de novas lojas a serem construídas.

E agora definimos as variáveis de decisão:

Vamos entender melhor o que são estas duas variáveis de decisão. A variável y_j representa todas as possíveis localizações para as p novas lojas. Caso o bairro j seja escolhido para receber uma nova loja, a variávei y_j recebe o valor 1. Caso contrário, recebe 0. Por exemplo, se os bairros 2, 5 e 8 fossem escolhidos para receberem as novas lojas, y_j receberia os seguintes valores:

Já a variável x_ij representa as possíveis alocações de cada bairro à loja mais próxima. Caso o bairro i seja alocado à loja j, a variável x_ij recebe o valor 1. Caso contrário, recebe 0. A ilustração abaixo representa os bairros 0, 1 e 2 alocados à loja do bairro 2, os bairros 3, 5 e 9 alocados à loja do bairro 5 e os bairros 4, 6, 7 e 8 alocados à loja do bairro 8.

A função objetivo fica definida da seguinte forma:

Aqui basicamente estamos somando o produto entre distância e alocação com o objetivo de encontrar as variáveis x_ij que minimizam esta soma.

Definida a função objetivo, precisamos estabelecer as restrições do problema. O primeiro conjunto de restrições garante que exatamente p bairros sejam escolhidos para receberem as novas lojas:

Observando a matriz x_ij do exemplo acima, percebemos que a soma de cada linha resulta em 1. Isto se deve ao conjunto de restrições abaixo, que tem como objetivo garantir que um bairro seja alocado a uma única loja:

Já o terceiro conjunto de restrições estabelece o “diálogo” entre y_j e x_ij. O objetivo aqui é garantir que os bairros só sejam alocados aos nós escolhidos como mediana:

E, por fim, as restrições abaixo garantem que as variáveis de decisão sejam binárias:

O modelo completo fica da seguinte forma:

Modelagem e resolução através do Pyomo

Agora que temos o modelo formulado, podemos partir para o Pyomo. Nesta etapa basicamente iremos reescrever, através da sintaxe do Python e do Pyomo, o modelo que acabemos de formular. Como as bibliotecas já foram carregadas, iniciamos declarando o modelo:

O próximo passo é a criação dos índices de linhas e colunas e a definição dos parâmetros do modelo:

Agora definimos as variáveis de decisão:

Em seguida criamos a função objetivo:

Agora vamos às restrições:

As restrições abaixo já foram definidas com as variáveis de decisão, através do argumento within=pyo.Binary :

Tudo certo, agora vamos executar o solver. Para o exemplo escolhi utilizar o GLPK:

Os dados da saída do modelo são intuitivos. Através deles podemos avaliar se o modelo chegou a uma solução factível, o tempo que ele levou para alcançar a solução, o valor da função objetivo (que no neste caso é a distância total entre os bairros e as novas lojas), a quantidade de restrições e etc. Para verificar quais bairros foram escolhidos como mediana podemos utilizar o comando abaixo:

Aqui percebemos que os bairros escolhidos para receber as novas lojas foram 0, 1 e 9. No entanto, esta forma de exibir os dados não é conveniente quando trabalhamos com problemas maiores. Uma alternativa é a seguinte:

Eu particularmente prefiro incluir no conjunto de dados uma coluna indicando se o bairro foi escolhido ou não como mediana.

Da mesma forma podemos imprimir as variáveis x_ij, mas não é muito conveniente analisar as alocações desta maneira:

Podemos imprimir somente as alocações realizadas:

Cada tupla representa a alocação do bairro i à mediana j. Aqui obervamos que os bairros 0 e 8 foram alocados à loja do bairro 0, os bairros 1, 2, 4 e 6 à loja do bairro 1 e os bairros 3, 5 e 9 à loja do bairro 9. Novamente opto por inserir estas informações no conjunto de dados:

Agora conseguimos visualizar quais bairros foram escolhidos como mediana e a qual loja cada bairro foi alocado. Adicionalmente podemos incluir no conjunto de dados as distâncias entre cada bairro e suas respectivas medianas:

Pegando como exemplo a 5º linha, o valor de 106,042444 representa a distância entre o bairro 4 e a loja do bairro 1, ao qual ele foi alocado.

Vamos agora visualizar o gráfico com as alocações:

Podemos somar a distância total por mediana:

E também conferir se a distância total está conforme a saída do modelo:

Exemplo 2

Imagine agora que a tal rede de lojas deseja acrescentar mais alguns dados no modelo para melhorar a tomada de decisão. Antes de tudo vamos carregar os novos dados:

A coluna Demanda representa a quantidade de clientes por bairro e a Capacidade, a capacidade de loja atender a demanda, caso ela seja construída. O objetivo agora é definir onde construir as 3 novas lojas com base não só nas distâncias, mas também na demanda de cada bairro, respeitando a capacidade das lojas.

O modelo

Precisamos realizar as modificações no modelo matemático. A primeira delas é a inclusão dos novos parâmetros:

O parâmetro C_i representa a demanda por bairro, enquanto K_j representa capacidade da loja, caso ela seja construída.

Outra modificação necessária é na função objetivo:

Observe que agora a função objetivo pondera as distâncias pela demanda de cada bairro e o objetivo é encontrar as variáveis x_ij que minimizam a soma dos produtos.

Uma outra modificação se faz necessária. Precisamos garantir que a soma das demandas do conjunto de bairros alocados à loja j não supere a capacidade da respectiva loja. Portanto, a restrição

deve ser substituída por

O modelo completo fica da seguinte forma:

Modelagem e resolução através do Pyomo

A declaração do modelo, os índices e os parâmetros utilizados no modelo anterior são exatamente os mesmos:

Mas agora temos mais parâmetros para adicionar no modelo. Precisamos acrescentar as demandas e as capacidades:

As variáveis de decisão, função objetivo e os conjuntos de restrições A e B também permanecem da mesma forma:

Já o conjunto de restrições C, conforme já vimos, foi alterado. Portanto, precisamos efetuar as modificações no Pyomo:

Pronto, podemos executar o solver:

Agora incluímos os dados da solução no conjunto de dados:

Para facilitar a nossa análise, vamos adicionar o cálculo da distância total no conjunto de dados:

E agora podemos visualizar o gráfico com as alocações:

Ao comparar com o primeiro exemplo, percebemos que a loja do bairro 0 foi substituída pela do bairro 8.

Observamos que a alteração da solução com relação ao primeiro exemplo não afetou a soma das distâncias, mas é certo que a nova solução agora garante que as lojas possuem capacidade suficiente para atender as demandas.

Exemplo 3

Agora vamos mudar um pouco o problema. Nos exemplos anteriores tínhamos definida a quantidade de novas lojas a serem construídas. Agora imagine que ao invés disso, nós temos um orçamento disponível de 50.000 dinheiros para construir p novas lojas e precisamos decidir quantas lojas construir e onde construir cada uma delas. Vamos carregar o novo conjunto de dados:

Observe que agora temos a coluna Custo, que indica o custo para construir cada uma das lojas, e que este custo varia em função do bairro. O objetivo agora, então, é decidir a quantidade de novas lojas e onde cada uma delas deverá ser construída, respeitando orçamento e minimizando as distâncias ponderadas pela demanda.

O modelo

Vamos realizar algumas modificações no modelo matemático anterior. Primeiro precisamos eliminar o parâmetro p e acrescentar B, que representa o orçamento disponível para construção das novas lojas:

Também devemos incluir no modelo o parâmetro que representa o custo para construção das novas lojas:

O conjunto de restrições

Deve ser eliminado do modelo, para dar lugar ao conjunto

que garante que a soma dos custos das novas lojas não supere o orçamento disponível.

O modelo completo fica da seguinte forma:

Modelagem e resolução através do Pyomo

A declaração do modelo, os índices e os parâmetros utilizados no exemplo anterior são basicamente os mesmos, só não podemos esquecer de remover o parâmetro p:

Precisamos incluir os novos parâmetros:

As variáveis de decisão, função objetivo e os conjuntos de restrições B e C também seguem exatamente como no exemplo anterior:

No entanto, é necessário alterar o conjunto de restrições A:

Pronto, podemos executar o solver:

Vamos incluir os dados da solução no dataset:

E agora visualizar as alocações no gráfico de dispersão:

Agora observamos que dado o orçamento de 50.000 unidades monetárias, o número de novas lojas que minimiza a soma das distâncias ponderadas pela demanda é 5, alocadas nos bairros 1, 4, 5, 7 e 8.

Vamos analisar os dados resumidos:

Aqui verificamos que o custo total da instalação das 5 novas lojas é de 48.214 unidades monetárias, ou seja, ainda sobram 1.786 dos 50.000 disponíveis. Também podemos notar que a distância total é de 17.603, significantemente menor que valor de 37.881 apresentado pelo modelo do exemplo anterior. Sob o ponto de vista das distâncias, o novo modelo apresenta uma melhora considerável, no entanto, esta redução possui uma relação direta com o número de novas lojas.

No exemplo anterior o modelo prescreveu a construção de 3 novas lojas nos bairros 1, 8 e 9 e se calcularmos o custo da construção (13.991 + 8.884 + 12.634 = 35.509) observamos que o modelo anterior, sob o ponto de vista econômico, é mais interessante que o atual.

Será que existe um meio termo? Se custo para a construção das 3 lojas do modelo anterior foi de 35k e o custo para a construção das 5 lojas do modelo atual ficou em 48k, podemos reduzir o orçamento disponível e tentar encontrar uma solução intermediária.

Considerando um orçamento de 42.000, chegamos à seguinte solução:

Limitando o orçamento em 42.000 observamos que as lojas 5 e 7 da solução considerando 50.000 dão lugar à loja 9, o que, consequentemente, reduz o custo de construção das novas lojas para 41.116 e eleva a soma das distâncias ponderadas para 24.414.

O objetivo do modelo, da maneira como ele foi formulado, é minimizar as distâncias. Logo, quanto maior o orçamento disponível, maior a tendência do modelo prescrever a construção de mais lojas, umas vez que quanto mais lojas, menor a soma das distâncias.

Vamos rodar o modelo novamente, agora considerando um orçamento total de 70,000:

Com um orçamento maior, a solução que minimiza a soma das distâncias prescreve a construção de 6 novas lojas, com um custo total de 62.139.

Exemplo 4

Percebemos na última variação do modelo anterior que as lojas 3, 5 e 7 são relativamente próximas. Supondo que um dos requisitos do negócio é que somente uma destas 3 lojas sejam construídas, ou nenhuma delas, vamos então alterar o nosso modelo para atender esta necessidade.

Recapitulando, o objetivo agora é construir p novas lojas com um orçamento de 70.000, entretanto, não é permitido que as lojas 3, 5 e 7 sejam construídas simultaneamente.

O modelo

Vamos aproveitar o modelo anterior, porém, precisamos acrescentar um novo parâmetro para representar o grupo das lojas 3, 5 e 7.

Em seguida devemos acrescentar a restrição que garante que no máximo uma das lojas do grupo T seja construída.

O modelo matemático completo fica definido da seguinte forma:

Modelagem e resolução através do Pyomo

Vamos aproveitar os índices e os parâmetros do modelo anterior:

Em seguida declaramos o parâmetro T:

As variáveis de decisão, função objetivo e as restrições também são exatamente as mesmas do modelo anterior:

Por fim incluímos a nova restrição:

Tudo pronto, podemos executar o solver e analisar os dados da saída do modelo:

Podemos notar que do grupo das lojas 3,5 e 7, somente a loja 5 foi escolhida como mediana. Em compensação, o custo de instalação das lojas aumentou de 62.139 para 65.630 e a função objetivo foi de 12.315 para 13.876, em comparação com o modelo sem esta nova restrição. No entanto, como não estamos lidando com problemas reais, fica impossível afirmar qual é o melhor modelo, uma vez que esta decisão depende dos requisitos e condições do negócio.

Conclusões

O objetivo deste artigo foi demonstrar o funcionamento do Pyomo através do problema das p-Medianas. As demonstrações realizadas deixaram evidente que os recursos do pacote e da linguagem Python tornam relativamente simples a modelagem de problemas de Programação Linear. É importante que o usuário esteja bem familiarizado com a linguagem Python para ser capaz de modelar problemas menos convencionais, principalmente pelo fato de que o Pyomo ainda não possui um material tão rico disponível para consulta na internet, quando comparado com outras ferramentas. A boa notícia é que Python, sem dúvidas, é uma das linguagens mais fáceis de se aprender, o que torna possível qualquer usuário, independentemente da sua formação, com um pouco de esforço aprender o suficiente em um curto espaço de tempo.

Para quem está começando no mundo da Programação Matemática, ficou clara também a importância de ter o modelo muito bem formulado e compreendido, para depois transformar tudo em código. Inverter a ordem das coisas pode trazer grandes complicações para a resolução do problema. Também foi constatado que não há limites para adicionar complexidade nos modelos, mas que iniciar com um modelo mais simples e complicar gradativamente, se necessário, pode ser o melhor caminho.

Em situações reais, utilizando as melhores práticas da programação é possível desenvolver soluções extremamente eficientes com Python e Pyomo, entretanto, por motivos didáticos, a abordagem na resolução dos problemas apresentados neste artigo ignorou em alguns sentidos tais boas práticas.

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

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.

LOCH, G. V. Problema das p-Medianas utilizando C# e Gurobi. 2020. Disponível em: <https://www.youtube.com/watch?v=dk0LHa_UVKc&list=PLTv9E12xrX8Om5Brh3PDXA_XN7amnS8Lc>. Acesso em: 12 de jan. 2021.

Pyomo. Pyomo Documentation 5.7.2. 2017. Disponível em: <https://pyomo.readthedocs.io/en/stable/>. Acesso em: 12 de jan. 2021.

Sandia National Laboratories. Pyomo. 2011. Disponível em: <https://www.osti.gov/servlets/purl/1376827>. Acesso em: 12 de jan. 2021.

Cientista de dados @ Grupo Boticário, interesse em modelos lineares generalizados e programação matemática, músico e lifelong learner.

Cientista de dados @ Grupo Boticário, interesse em modelos lineares generalizados e programação matemática, músico e lifelong learner.