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:
import pyomo.environ as pyo
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.neighbors import DistanceMetricdf_dados = pd.read_csv("lojas_1.csv", delimiter=";")
df_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:
dist = DistanceMetric.get_metric('euclidean')
dij = df_dados[['X', 'Y']].to_numpy()
dij = dist.pairwise(dij)
pd.DataFrame(dij)
Visualizando a dispersão dos bairros:
plt.figure(figsize=(12, 6))
plot = sns.scatterplot(data=df_dados, x="X", y="Y", s=100)
for index, row in df_dados.iterrows():
plot.text(row['X']+3, row['Y'], index, horizontalalignment='left')
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:
modelo = pyo.ConcreteModel()O próximo passo é a criação dos índices de linhas e colunas e a definição dos parâmetros do modelo:
# Índices:
modelo.M = range(len(dij))
modelo.N = range(len(dij))# Parâmetros
modelo.d = pyo.Param(modelo.M, modelo.N, initialize=lambda modelo, i, j: dij[i][j])
p = 3
Agora definimos as variáveis de decisão:
modelo.y = pyo.Var(modelo.N, within=pyo.Binary)
modelo.x = pyo.Var(modelo.M, modelo.N, within=pyo.Binary)Em seguida criamos a função objetivo:
def f_obj(modelo):
return sum(modelo.x[i,j] * modelo.d[i,j] for i in modelo.M for j in modelo.N)modelo.objetivo = pyo.Objective(rule=f_obj, sense=pyo.minimize)
Agora vamos às restrições:
modelo.restricao_a = pyo.Constraint(expr=sum(modelo.y[j] for j in modelo.N) == p)
modelo.restricao_b = pyo.ConstraintList()for i in modelo.M:
modelo.restricao_b.add(sum(modelo.x[i,j] for j in modelo.N) == 1)
modelo.restricao_c = pyo.ConstraintList()for i in modelo.M:
for j in modelo.N:
modelo.restricao_c.add(modelo.x[i,j] <= modelo.y[j])
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:
resultado = pyo.SolverFactory('glpk').solve(modelo)
print(resultado)
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:
modelo.y.pprint()
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:
list_y = list(modelo.y.keys())
[j for j in list_y if modelo.y[j]() == 1]
Eu particularmente prefiro incluir no conjunto de dados uma coluna indicando se o bairro foi escolhido ou não como mediana.
dados_modelo = df_dados.copy()
dados_modelo['Mediana'] = [modelo.y[i]() for i in list_y]
dados_modelo
Da mesma forma podemos imprimir as variáveis x_ij, mas não é muito conveniente analisar as alocações desta maneira:
modelo.x.pprint()Podemos imprimir somente as alocações realizadas:
list_x = list(modelo.x.keys())
alocacoes = [i for i in list_x if modelo.x[i]() == 1]
alocacoes.sort(key=lambda x:x[0])
print(alocacoes)
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:
medianas = [alocacao[1] for alocacao in alocacoes]
dados_modelo['Alocacao'] = medianas
dados_modelo
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:
dados_modelo['Distancia'] = [dij[alocacao[0], alocacao[1]] for alocacao in alocacoes]
dados_modelo
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:
plt.figure(figsize=(12, 6))
plot = sns.scatterplot(data=dados_modelo,
x="X", y="Y",
hue="Alocacao",
alpha=.7,
s=100,
palette="tab10")
highlights = dados_modelo[dados_modelo.index == dados_modelo.Alocacao]
for index, row in highlights.iterrows():
plot.text(row['X']+3, row['Y'], row['Alocacao'], horizontalalignment='left')
Podemos somar a distância total por mediana:
dados_modelo_resumo = dados_modelo.copy()
resumo = dados_modelo_resumo.groupby('Alocacao', as_index=False).agg({"Distancia": "sum"})
resumo
E também conferir se a distância total está conforme a saída do modelo:
resumo.Distancia.sum()
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:
df_dados = pd.read_csv("lojas_2.csv", delimiter=";")
print(df_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:
# Declaração do modelo:
modelo2 = pyo.ConcreteModel()# Matriz de distâncias:
dist = DistanceMetric.get_metric('euclidean')
dij = df_dados[['X', 'Y']].to_numpy()
dij = dist.pairwise(dij)
# Índices:
modelo2.M = range(len(dij))
modelo2.N = range(len(dij))# Parâmetros:
modelo2.d = pyo.Param(modelo2.M, modelo2.N, initialize=lambda modelo, i, j: dij[i][j])
p = 3
Mas agora temos mais parâmetros para adicionar no modelo. Precisamos acrescentar as demandas e as capacidades:
# Vetores de demanda e capacidade
Ci = list(df_dados.Demanda)
Kj = list(df_dados.Capacidade)modelo2.C = pyo.Param(modelo2.M, initialize=lambda modelo2, i: Ci[i])
modelo2.K = pyo.Param(modelo2.N, initialize=lambda modelo2, j: Kj[j])
As variáveis de decisão, função objetivo e os conjuntos de restrições A e B também permanecem da mesma forma:
# Variáveis de decisão:
modelo2.y = pyo.Var(modelo2.N, within=pyo.Binary)
modelo2.x = pyo.Var(modelo2.M, modelo2.N, within=pyo.Binary)# Função objetivo:
def f_obj(modelo):
return sum(modelo.C[i] * modelo.x[i,j] * modelo.d[i,j] for i in modelo.M for j in modelo.N)modelo2.obj = pyo.Objective(rule=f_obj, sense=pyo.minimize)# Sujeito a:modelo2.restricao_a = pyo.Constraint(expr=sum(modelo2.y[j] for j in modelo2.N) == p)modelo2.restricao_b = pyo.ConstraintList()
for i in modelo2.M:
modelo2.restricao_b.add(sum(modelo2.x[i,j] for j in modelo2.N) == 1.0)
Já o conjunto de restrições C, conforme já vimos, foi alterado. Portanto, precisamos efetuar as modificações no Pyomo:
modelo2.restricao_c = pyo.ConstraintList()
for j in modelo2.N:
modelo2.restricao_c.add(sum(modelo2.C[i] * modelo2.x[i,j] for i in modelo2.M) <= modelo2.K[j] * modelo2.y[j])Pronto, podemos executar o solver:
resultado2 = pyo.SolverFactory('glpk').solve(modelo2)
print(resultado2)
Agora incluímos os dados da solução no conjunto de dados:
# Medianas
list_y = list(modelo2.y.keys())
dados_modelo2 = df_dados.copy()
dados_modelo2['Mediana'] = [modelo2.y[i]() for i in list_y]# Alocações:
list_x = list(modelo2.x.keys())
alocacoes = [i for i in list_x if modelo2.x[i]() == 1]
alocacoes.sort(key=lambda x:x[0])
medianas = [alocacao[1] for alocacao in alocacoes]
dados_modelo2['Alocacao'] = medianas# Distância
dados_modelo2['Distancia'] = [dij[alocacao[0], alocacao[1]] for alocacao in alocacoes]dados_modelo2
Para facilitar a nossa análise, vamos adicionar o cálculo da distância total no conjunto de dados:
dados_modelo2['Distancia_total'] = dados_modelo2['Distancia'] * dados_modelo2['Demanda']
print(dados_modelo2)
E agora podemos visualizar o gráfico com as alocações:
plt.figure(figsize=(12, 6))
plot = sns.scatterplot(data=dados_modelo2,
x="X", y="Y",
hue="Alocacao",
size='Demanda',
sizes=(50, 250),
alpha=.7,
palette="tab10")
highlights = dados_modelo2[dados_modelo2.index == dados_modelo2.Alocacao]
for index, row in highlights.iterrows():
plot.text(row['X']+3, row['Y'], row['Alocacao'], horizontalalignment='left')
Ao comparar com o primeiro exemplo, percebemos que a loja do bairro 0 foi substituída pela do bairro 8.
dados_modelo2_resumo = dados_modelo2.copy()
resumo2 = dados_modelo2_resumo.groupby('Alocacao', as_index=False).agg({"Demanda": "sum",
"Distancia": "sum",
"Distancia_total": "sum"})
sum_dist2 = resumo2.Distancia.sum()
sum_dist_total2 = resumo2.Distancia_total.sum()
print("Distância:", sum_dist2)
print("Distância_total:", sum_dist_total2)
resumo2
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:
df_dados = pd.read_csv("lojas_3.csv", delimiter=";")
df_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:
# Declaração do modelo:
modelo3 = pyo.ConcreteModel()# Matriz de distâncias:
dist = DistanceMetric.get_metric('euclidean')
dij = df_dados[['X', 'Y']].to_numpy()
dij = dist.pairwise(dij)
# Índices:
modelo3.M = range(len(dij))
modelo3.N = range(len(dij))# Parâmetros:# Vetores de demanda e capacidade
Ci = list(df_dados.Demanda)
Kj = list(df_dados.Capacidade)modelo3.d = pyo.Param(modelo3.M, modelo3.N, initialize=lambda modelo, i, j: dij[i][j])
modelo3.C = pyo.Param(modelo3.M, initialize=lambda modelo3, i: Ci[i])
modelo3.K = pyo.Param(modelo3.N, initialize=lambda modelo3, j: Kj[j])
Precisamos incluir os novos parâmetros:
# Vetor de custos:
Sj = list(df_dados.Custo)modelo3.S = pyo.Param(modelo3.N, initialize=lambda modelo, j: Sj[j])
B = 50000
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:
# Variáveis de decisão:
modelo3.y = pyo.Var(modelo3.N, within=pyo.Binary)
modelo3.x = pyo.Var(modelo3.M, modelo3.N, within=pyo.Binary)# Função objetivo:
def f_obj(modelo):
return sum(modelo.C[i] * modelo.x[i,j] * modelo.d[i,j] for i in modelo.M for j in modelo.N)modelo3.obj = pyo.Objective(rule=f_obj, sense=pyo.minimize)# Sujeito a:modelo3.restricao_b = pyo.ConstraintList()
for i in modelo3.M:
modelo3.restricao_b.add(sum(modelo3.x[i,j] for j in modelo3.N) == 1.0)
modelo3.restricao_c = pyo.ConstraintList()
for j in modelo3.N:
modelo3.restricao_c.add(sum(modelo3.C[i] * modelo3.x[i,j] for i in modelo3.M) <= modelo3.K[j] * modelo3.y[j])
No entanto, é necessário alterar o conjunto de restrições A:
modelo3.restricao_a = pyo.Constraint(expr=sum(modelo3.S[j] * modelo3.y[j] for j in modelo3.N) <= B)Pronto, podemos executar o solver:
resultado3 = pyo.SolverFactory('glpk').solve(modelo3)
print(resultado3)
Vamos incluir os dados da solução no dataset:
# Medianas
list_y = list(modelo3.y.keys())
dados_modelo3 = df_dados.copy()
dados_modelo3['Mediana'] = [modelo3.y[i]() for i in list_y]# Alocações:
list_x = list(modelo3.x.keys())
alocacoes = [i for i in list_x if modelo3.x[i]() == 1]
alocacoes.sort(key=lambda x:x[0])
medianas = [alocacao[1] for alocacao in alocacoes]
dados_modelo3['Alocacao'] = medianas# Distâncias
dados_modelo3['Distancia'] = [dij[alocacao[0], alocacao[1]] for alocacao in alocacoes]
dados_modelo3['Distancia_total'] = dados_modelo3['Distancia'] * dados_modelo3['Demanda']dados_modelo3
E agora visualizar as alocações no gráfico de dispersão:
plt.figure(figsize=(12, 6))
plot = sns.scatterplot(data=dados_modelo3,
x="X", y="Y",
hue="Alocacao",
size='Demanda',
sizes=(50, 250),
alpha=.7,
palette="tab10")
highlights = dados_modelo3[dados_modelo3.index == dados_modelo3.Alocacao]
for index, row in highlights.iterrows():
plot.text(row['X']+3, row['Y'], row['Alocacao'], horizontalalignment='left')
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:
dados_modelo3_resumo = dados_modelo3.copy()
dados_modelo3_resumo['Custo'] = dados_modelo3_resumo['Custo'] * dados_modelo3_resumo['Mediana']
resumo3 = dados_modelo3_resumo.groupby('Alocacao', as_index=False).agg({"Demanda": "sum",
"Distancia": "sum",
"Distancia_total": "sum",
"Custo": "sum"})
sum_dist3 = resumo3.Distancia.sum()
sum_dist_total3 = resumo3.Distancia_total.sum()
sum_custo3 = resumo3.Custo.sum()
print("Distância:", sum_dist3)
print("Distância_total:", sum_dist_total3)
print("Custo:", sum_custo3)
resumo3
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:
# Declaração do modelo:
modelo4 = pyo.ConcreteModel()# Matriz de distâncias:
dist = DistanceMetric.get_metric('euclidean')
dij = df_dados[['X', 'Y']].to_numpy()
dij = dist.pairwise(dij)
# Índices:
modelo4.M = range(len(dij))
modelo4.N = range(len(dij))# Parâmetros:# Vetores de demanda, capacidade e custo
Ci = list(df_dados.Demanda)
Kj = list(df_dados.Capacidade)
Sj = list(df_dados.Custo)modelo4.d = pyo.Param(modelo4.M, modelo4.N, initialize=lambda modelo, i, j: dij[i][j])
modelo4.C = pyo.Param(modelo4.M, initialize=lambda modelo4, i: Ci[i])
modelo4.K = pyo.Param(modelo4.N, initialize=lambda modelo4, j: Kj[j])
modelo4.S = pyo.Param(modelo4.N, initialize=lambda modelo4, j: Sj[j])
B = 70000
Em seguida declaramos o parâmetro T:
T = [3,5,7]As variáveis de decisão, função objetivo e as restrições também são exatamente as mesmas do modelo anterior:
# Variáveis de decisão:
modelo4.y = pyo.Var(modelo4.N, within=pyo.Binary)
modelo4.x = pyo.Var(modelo4.M, modelo4.N, within=pyo.Binary)# Função objetivo:
def f_obj(modelo):
return sum(modelo.C[i] * modelo.x[i,j] * modelo.d[i,j] for i in modelo.M for j in modelo.N)modelo4.obj = pyo.Objective(rule=f_obj, sense=pyo.minimize)# Sujeito a:modelo4.restricao_a = pyo.Constraint(expr=sum(modelo4.S[j] * modelo4.y[j] for j in modelo4.N) <= B)modelo4.restricao_b = pyo.ConstraintList()
for i in modelo4.M:
modelo4.restricao_b.add(sum(modelo4.x[i,j] for j in modelo4.N) == 1.0)
modelo4.restricao_c = pyo.ConstraintList()
for j in modelo4.N:
modelo4.restricao_c.add(sum(modelo4.C[i] * modelo4.x[i,j] for i in modelo4.M) <= modelo4.K[j] * modelo4.y[j])
Por fim incluímos a nova restrição:
modelo4.restricao_d = pyo.Constraint(expr=sum(modelo4.y[j] for j in modelo4.N if j in T) <= 1)Tudo pronto, podemos executar o solver e analisar os dados da saída do modelo:
resultado4 = pyo.SolverFactory('glpk').solve(modelo4)
print(resultado4)
# Medianas
list_y = list(modelo4.y.keys())
dados_modelo4 = df_dados.copy()
dados_modelo4['Mediana'] = [modelo4.y[i]() for i in list_y]# Alocações:
list_x = list(modelo4.x.keys())
alocacoes = [i for i in list_x if modelo4.x[i]() == 1]
alocacoes.sort(key=lambda x:x[0])
medianas = [alocacao[1] for alocacao in alocacoes]
dados_modelo4['Alocacao'] = medianas# Distâncias
dados_modelo4['Distancia'] = [dij[alocacao[0], alocacao[1]] for alocacao in alocacoes]
dados_modelo4['Distancia_total'] = dados_modelo4['Distancia'] * dados_modelo4['Demanda']dados_modelo4
plt.figure(figsize=(12, 6))
plot = sns.scatterplot(data=dados_modelo4,
x="X", y="Y",
hue="Alocacao",
size='Demanda',
sizes=(50, 250),
alpha=.7,
palette="tab10")
highlights = dados_modelo4[dados_modelo4.index == dados_modelo4.Alocacao]
for index, row in highlights.iterrows():
plot.text(row['X']+3, row['Y'], row['Alocacao'], horizontalalignment='left')
dados_modelo4_resumo = dados_modelo4.copy()
dados_modelo4_resumo['Custo'] = dados_modelo4_resumo['Custo'] * dados_modelo4_resumo['Mediana']
resumo4 = dados_modelo4_resumo.groupby('Alocacao', as_index=False).agg({"Demanda": "sum",
"Distancia": "sum",
"Distancia_total": "sum",
"Custo": "sum"})
sum_dist4 = resumo4.Distancia.sum()
sum_dist_total4 = resumo4.Distancia_total.sum()
sum_custo4 = resumo4.Custo.sum()
print("Distância:", sum_dist4)
print("Distância_total:", sum_dist_total4)
print("Custo:", sum_custo4)
resumo4
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.
