Problema da mochila
Prof. Steven Dutt-Ross
Outro dia encontrei o meu caderno de uma disciplina de otimização que fiz no doutorado. Ao reler as anotações, encontrei o famoso problema da mochila (Knapsack Problem).
Esse era um problema muito interessante. Lembrei que tive que desenvolver um código em VBA (era o que eu usava na época) para resolver esse problema. Mas agora que eu sei R & Python, pensei em resolver o problema novamente com essas linguagens.
O meu primeiro pensamento foi desenvolver um pacote do R para resolver esse tipo de problema, mas ao pesquisar um pouco na web, percebi que já existiam pacotes para isso. Além disso, encontrei duas publicações em blogs que usam abordagens que eu nunca tinha pensado antes. Uma usa geoprocessamto e a outra usa algoritmo genético.
Assim, nessa publicação gostaria de apresentar o problema da mochila, bem como usar o R para resolver esse tipo de questão. No fim, vou tentar usar o método do Algoritmo Genético para resolver esse mesmo problema de uma forma diferente.
Para quem quiser dar uma olhada nessas duas referências, você pode ver o post do Julio Trecenti, do Curso-R sobre o problema da mochila aqui. Já sobre a integração do Algoritmo Genético com o problema da mochila pode ser vista aqui
O problema da mochila
Considere o seguinte contexto: você tem uma mochila com capacidade máxima de 25kg e precisa carregar a combinação de itens com maior valor possível , com cada item possuindo valores e pesos diferentes. Quais deles vocês escolheria? Provavelmente os item mais leves com valor maior.
Temos três elementos nesse problema: 1- Capacidade Máxima da mochila, 2- Peso dos itens, 3- Valor dos itens.
Desse modo, como exercício aplicado vamos pensar no estado do Rio de Janeiro como uma “mochila”. Neste exemplo aplicado, o objetivo é selecionar os municípios do estado do Rio de janeiro cuja soma seja igual ao PIB do Rio de Janeiro e ao mesmo tempo observar o tamanho (área em km2) do município como critério de seleção. Usando essa analogia temos:
- peso dos itens = área do municipio,
- valor dos itens = PIB,
- capacidade máxima = PIB do rio de janeiro.
Para resolver esse problema, em linguagem matemática, temos a seguinte tarefa:
\[\begin{aligned} & \text{maximizar } \sum_{i=1}^n v_i x_i \\ & \text{sujeito à } \sum_{i=1}^n w_i x_i \leq W, \text{ com } x_i \in\{0,1\}\\ \end{aligned}\]onde:
. n é o número de municípios no estado (91). Foi excluído o
município do Rio.
. vi é a área do município i
. wi é o PIB do município i
. W é o PIB do Rio de Janeiro (472 milhões).
Como funciona a função knapsack do pacote adagio
Para resolver o problema da mochila, podemos usar a função knapsack do pacote adagio. Vamos ver como ele funciona. Olhando o help do pacote, podemos ver que funciona de forma bem simples.
library(adagio)
# Example 1
p <- c(15, 100, 90, 60, 40, 15, 10, 1)
w <- c( 2, 20, 20, 30, 40, 30, 60, 10)
cap <- 102
(is <- knapsack(w, p, cap))$capacity
[1] 102
$profit
[1] 280
$indices
[1] 1 2 3 4 6# [1] 1 2 3 4 6 , capacity 102 and total profit 280
## Example 2
p <- c(70, 20, 39, 37, 7, 5, 10)
w <- c(31, 10, 20, 19, 4, 3, 6)
cap <- 50
(is <- knapsack(w, p, cap))$capacity
[1] 50
$profit
[1] 107
$indices
[1] 1 4# [1] 1 4 , capacity 50 and total profit 107O pacote tem três parâmetros chamados de pesos w,
valores p e máximo cap: 1.
w integer vector of weights.
2. p integer vector of profits.
3. cap maximal capacity of the knapsack, integer
too.
Resolvendo o problema da mochila usando algoritmo genético
Voltando ao problema de uma mochila com capacidade máxima de 25 quilos. O objetivo é maximizar a capacidade de mochila enquanto maximiza seus pontos de survival também. A partir do problema, podemos definir o peso dos itens como a função de restrição, enquanto os pontos survival acumulados dos itens na mochila como a função objetivo.
Para a implementação neste artigo, usamos aqui o pacote GA do R . Primeiro, precisamos inserir os dados e os parâmetros.
#0-1 Knapsack's Problem
library(GA)
item=c('capa de chuva','canivete','água mineral','luvas','saco de dormir','tenda','fogão portátil','comida enlatada','lanches')
weight=c(2,1,6,1,4,9,5,8,3)
survival=c(5,3,15,5,6,18,8,20,8)
data=data.frame(item,weight,survival)
max_weight=25Para ter uma melhor compreensão do algoritmo genético neste problema, suponha que, inicialmente, levemos apenas um canivete, água mineral e lanches em sua mochila.
Podemos escrevê-lo como o “cromossomo”, e como o problema que queremos resolver é um problema de mochila 0–1, então, a partir dessas linhas de códigos abaixo, 1 significa que trouxemos o item, enquanto 0 significa que deixamos o item fora da mochila.
# 1 significa que trouxemos o item, enquanto 0 significa que deixamos o item fora da mochila
chromosomes=c(0,1,1,0,0,0,0,0,1)
data[chromosomes==1,] item weight survival
2 canivete 1 3
3 água mineral 6 15
9 lanches 3 8Em seguida, criamos a função objetivo que queremos otimizar com a função de restrição escrevendo essas linhas de código abaixo. Preste atenção na função ga do pacote GA .
#create the function that we want to optimize
fitness=function(x){
current_survpoint=x%*%data$survival
current_weight=x%*%data$weight
if(current_weight>max_weight)
{
return(0)
}
else
{
return(current_survpoint)
}
}Agora aqui está a parte interessante: o processo de otimização usando o algoritmo genético. Suponha que queremos criar no máximo 30 gerações e 50 indivíduos para o processo de otimização. Para reprodutibilidade, vamos escrever a semente (seed=1234) para guardar a melhor solução.
GA=ga(type='binary',fitness=fitness,nBits=nrow(data),maxiter=30,popSize=50,seed=1234,keepBest=TRUE)
summary(GA)-- Genetic Algorithm -------------------
GA settings:
Type = binary
Population size = 50
Number of generations = 30
Elitism = 2
Crossover probability = 0.8
Mutation probability = 0.1
GA results:
Iterations = 30
Fitness function value = 61
Solution =
x1 x2 x3 x4 x5 x6 x7 x8 x9
[1,] 1 0 1 1 0 0 1 1 1plot(GA)
A partir do resultado acima, podemos ver que o desempenho de cada indivíduo está aumentando a cada geração. Podemos vê-lo a partir do valor de fitness (ou seja, os pontos de sobrevivência neste caso) média e mediana que tende a aumentar a cada geração. Vamos tentar treiná-lo novamente, mas com mais gerações.
GA2=ga(type='binary',fitness=fitness,nBits=nrow(data),maxiter=50,popSize=50,seed=1234,keepBest=TRUE)
GA3=ga(type='binary',fitness=fitness,nBits=nrow(data),maxiter=100,popSize=50,seed=1234,keepBest=TRUE)
plot(GA2)
plot(GA3)
A partir de GA2 e GA3, podemos ver que o resultado da otimização para cada indivíduo é melhor na geração 40-ish e 60-ish, de acordo com a média e a mediana do valor de fitness naquela geração. Também podemos ver que o melhor valor de fitness está aumentando para 62 a partir da 72ª geração.
Como mantemos o melhor resultado em cada processo de otimização, queremos descobrir os itens que podemos trazer para caminhadas com base no melhor resultado da otimização do algoritmo genético. Podemos ver o resumo do GA3 da seguinte forma.
summary(GA3)-- Genetic Algorithm -------------------
GA settings:
Type = binary
Population size = 50
Number of generations = 100
Elitism = 2
Crossover probability = 0.8
Mutation probability = 0.1
GA results:
Iterations = 100
Fitness function value = 62
Solution =
x1 x2 x3 x4 x5 x6 x7 x8 x9
[1,] 1 1 1 1 1 0 0 1 1Do resultado acima, podemos concluir que os itens que podemos incluir na mochila são capa de chuva, canivete, água mineral, luvas, saco de dormir, comida enlatada e lanches. Podemos calcular o peso da mochila, para garantir que a mochila não esteja capacitade superada.
chromosomes_final=c(1,1,1,1,1,0,0,1,1)
cat(GA3@solution%*%data$weight)25Podemos ver que o peso dos itens é o mesmo que a capacidade da mochila!
O problema da mochila para o PIB dos municípios
Agora que conhecemos o problema da mochila e a sua solução com os pacotes adagio e GA , vamos usar essas duas abordagens para o nosso problema. Antes, vamos criar um data.frame com o código abaixo.
dados_pib <-data.frame(
Munic = c("Rio de Janeiro",
"Belford Roxo","Cachoeiras de Macacu","Duque de Caxias",
"Guapimirim","Itaboraí","Itaguaí","Japeri","Magé","Maricá",
"Mesquita","Nilópolis","Niterói","Nova Iguaçu",
"Paracambi","Petrópolis","Queimados","Rio Bonito",
"São Gonçalo","São João de Meriti","Seropédica","Tanguá",
"Aperibé","Bom Jesus do Itabapoana","Cambuci","Italva",
"Itaocara","Itaperuna","Laje do Muriaé","Miracema",
"Natividade","Porciúncula","Santo Antônio de Pádua",
"São José de Ubá","Varre-Sai","Campos dos Goytacazes",
"Carapebus","Cardoso Moreira","Conceição de Macabu","Macaé",
"Quissamã","São Fidélis","São Francisco de Itabapoana",
"São João da Barra","Região Serrana","Bom Jardim",
"Cantagalo","Carmo","Cordeiro","Duas Barras","Macuco",
"Nova Friburgo","Santa Maria Madalena",
"São José do Vale do Rio Preto","São Sebastião do Alto","Sumidouro",
"Teresópolis","Trajano de Moraes","Araruama",
"Armação dos Búzios","Arraial do Cabo","Cabo Frio",
"Casimiro de Abreu","Iguaba Grande","Rio das Ostras",
"São Pedro da Aldeia","Saquarema","Silva Jardim","Barra do Piraí",
"Barra Mansa","Itatiaia","Pinheiral","Piraí","Porto Real",
"Quatis","Resende","Rio Claro","Rio das Flores",
"Valença","Volta Redonda","Areal",
"Comendador Levy Gasparian","Engenheiro Paulo de Frontin","Mendes",
"Miguel Pereira","Paraíba do Sul","Paty do Alferes","Sapucaia",
"Três Rios","Vassouras","Angra dos Reis","Mangaratiba",
"Paraty"),
PIB = c(27179971786,731036095,
105129677,3884394488,99072237,451066428,722525751,132881462,
401686294,3708656972,220341754,273349746,4381265583,
1593313298,88252114,1200444417,345470022,141314651,
1701165961,896614838,360708948,54095355,18376438,
67253563,31279571,31112788,50582119,287135565,13594979,
52961166,27813764,32198816,106386843,15880774,
18130274,2806792113,43919000,27290047,35325470,
1298785681,377269292,69024389,90857102,756448858,1359715951,
62505140,68043504,42078057,37870527,19581447,
15040039,489397947,17516564,46777259,18174228,36455894,
487334211,18941134,343403620,253119881,250113318,
1111851510,200194515,51047591,738470457,219449129,
1051275133,45821739,194697579,498430523,286023242,64062968,
199976589,165709937,25792031,634577701,34775592,
20454076,186028436,942319023,30264340,38585272,23437100,
32800832,57938201,91268018,50910950,62330320,
377489712,101261986,930305905,234037012,266306789)
)
dados_area <-
data.frame(
stringsAsFactors = FALSE,
Munic = c("Campos dos Goytacazes",
"Valença","Macaé","Rio de Janeiro",
"São Francisco de Itabapoana","Itaperuna","Resende","São Fidélis",
"Cachoeiras de Macacu","Silva Jardim","Nova Friburgo",
"Paraty","Rio Claro","Angra dos Reis",
"Santa Maria Madalena","Petrópolis","Teresópolis","Cantagalo",
"Quissamã","Araruama","Santo Antônio de Pádua",
"Bom Jesus do Itabapoana","Trajano de Moraes","Barra do Piraí",
"Paraíba do Sul","Cambuci","Barra Mansa","Sapucaia",
"Vassouras","Cardoso Moreira","Nova Iguaçu",
"Piraí","Rio das Flores","Duque de Caxias",
"Casimiro de Abreu","Rio Bonito","São João da Barra","Itaocara",
"Itaboraí","Cabo Frio","Sumidouro",
"São Sebastião do Alto","Magé","Natividade","Bom Jardim","Duas Barras",
"Mangaratiba","Maricá","Guapimirim","Saquarema",
"Conceição de Macabu","São Pedro da Aldeia",
"Três Rios","Paty do Alferes","Carmo","Carapebus","Miracema",
"Porciúncula","Italva","Miguel Pereira","Quatis",
"Itaguaí","Seropédica","Laje do Muriaé",
"São José de Ubá","São Gonçalo","Itatiaia","Rio das Ostras",
"São José do Vale do Rio Preto","Varre-Sai","Paracambi",
"Volta Redonda","Arraial do Cabo","Tanguá",
"Engenheiro Paulo de Frontin","Niterói","Cordeiro","Areal",
"Comendador Levy Gasparian","Mendes","Aperibé",
"Pinheiral","Japeri","Belford Roxo","Macuco","Queimados",
"Armação dos Búzios","Iguaba Grande","Porto Real",
"Mesquita","São João de Meriti","Nilópolis"),
area = c(403249,130077,121699,
120033,111804,110669,109934,103483,95475,93776,
93543,92430,84680,81321,81096,79114,77334,74721,
71964,63815,60363,59666,59115,58461,57112,55828,
54713,54067,53607,52260,52058,49026,47878,46732,
46292,45946,45240,43318,42996,41358,41341,39721,
39078,38707,38243,37962,36782,36157,35844,35213,
33826,33249,32284,31434,30575,30489,30327,29185,
29119,28793,28483,28261,26519,25353,24969,
24816,24104,22804,22018,20194,19095,18211,15211,
14301,13938,13376,11305,11072,10864,9532,9454,8225,
8170,7899,7836,7593,7098,5098,5089,4117,3522,
1939)
)
library(dplyr)
dados <- dados_area %>% left_join(dados_pib)
remove(dados_pib,dados_area)
dados_sem_rj <- dados %>%filter(Munic!="Rio de Janeiro")library(adagio)
is <- knapsack(dados_sem_rj$area, dados_sem_rj$PIB, 403249)
is$capacity[1] 401483is$profit[1] 23477329307is$indices [1] 30 33 36 39 47 49 61 65 67 71 75 83 85 86 88 89 90 91dados_sem_rj$indices <-rownames(dados_sem_rj)
dados_sem_rj$indices<- as.integer(dados_sem_rj$indices)
selecao <- is$indices
selecao <- data.frame(selecao)
colnames(selecao) <-"indices"
selecao$incluir <- 1
library(dplyr)
dados_sem_rj <- left_join(dados_sem_rj,selecao)
dados_sem_rj$incluir <- ifelse(is.na(dados_sem_rj$incluir),0,1)
dados<-left_join(dados,dados_sem_rj)
dados$incluir<-ifelse(is.na(dados$incluir),2,dados$incluir)De acordo com o pacote Adagio, os “municípios na mochila” são:
dados_sem_rj %>% filter(incluir==1) %>% select(Munic) Munic
1 Nova Iguaçu
2 Duque de Caxias
3 São João da Barra
4 Cabo Frio
5 Maricá
6 Saquarema
7 Itaguaí
8 São Gonçalo
9 Rio das Ostras
10 Volta Redonda
11 Niterói
12 Belford Roxo
13 Queimados
14 Armação dos Búzios
15 Porto Real
16 Mesquita
17 São João de Meriti
18 NilópolisPara melhor visualização, vamos colocar o resultado em um mapa.
library(geobr)
malha <- read_municipality(code_muni=33, year=2010)
colnames(dados)[1]<-"name_muni"
malha2<-malha %>% left_join(dados)
ggplot2::ggplot(malha2) +
ggplot2::geom_sf(
ggplot2::aes(fill = as.factor(incluir)),
size = 0
) +
ggplot2::theme_void() +
ggplot2::theme(legend.position = "bottom") +
ggplot2::labs(fill = "Solução")
Downloading: 770 B
Downloading: 770 B
Downloading: 1.8 kB
Downloading: 1.8 kB
Downloading: 2.1 kB
Downloading: 2.1 kB
Downloading: 2.1 kB
Downloading: 2.1 kB
Downloading: 2.1 kB
Downloading: 2.1 kB
Downloading: 5.5 kB
Downloading: 5.5 kB
Downloading: 5.5 kB
Downloading: 5.5 kB
Downloading: 14 kB
Downloading: 14 kB
Downloading: 30 kB
Downloading: 30 kB
Downloading: 30 kB
Downloading: 30 kB
Downloading: 38 kB
Downloading: 38 kB
Downloading: 54 kB
Downloading: 54 kB
Downloading: 54 kB
Downloading: 54 kB
Downloading: 54 kB
Downloading: 54 kB
Downloading: 70 kB
Downloading: 70 kB
Downloading: 70 kB
Downloading: 70 kB
Downloading: 86 kB
Downloading: 86 kB
Downloading: 86 kB
Downloading: 86 kB
Downloading: 100 kB
Downloading: 100 kB
Downloading: 110 kB
Downloading: 110 kB
Downloading: 130 kB
Downloading: 130 kB
Downloading: 130 kB
Downloading: 130 kB
Downloading: 130 kB
Downloading: 130 kB
Downloading: 140 kB
Downloading: 140 kB
Downloading: 140 kB
Downloading: 140 kB
Downloading: 150 kB
Downloading: 150 kB
Downloading: 170 kB
Downloading: 170 kB
Downloading: 170 kB
Downloading: 170 kB
Downloading: 180 kB
Downloading: 180 kB
Downloading: 190 kB
Downloading: 190 kB
Downloading: 190 kB
Downloading: 190 kB
Downloading: 210 kB
Downloading: 210 kB
Downloading: 210 kB
Downloading: 210 kB
Downloading: 220 kB
Downloading: 220 kB
Downloading: 230 kB
Downloading: 230 kB
Downloading: 230 kB
Downloading: 230 kB 
Agora vamos tentar resolver o problema da mochila com o algoritmo genético. Do mesmo modo que o processo anterior, vamos começar com “Valença”, “Macaé”, e”São Francisco de Itabapoana” na mochila.
chromosomes=c(0,1,1,1,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
dados_sem_rj[chromosomes==1,1][1] "Valença" "Macaé"
[3] "São Francisco de Itabapoana"Do mesmo modo que o passo apresentado no segmento anterior, vamos definar a capacidade máxima e a função fitness.
#create the function that we want to optimize
max_weight =403249
fitness=function(x){
current_pib=x%*%dados_sem_rj$PIB
current_area=x%*%dados_sem_rj$area
if(current_area>max_weight)
{
return(0)
}
else
{
return(current_pib)
}
}Aqui vamos aplicar a função GA para ver o funcionamento nos nosso dados.
GA=ga(type='binary',fitness=fitness,nBits=nrow(dados_sem_rj),maxiter=150,popSize=100,seed=1234,keepBest=TRUE)
summary(GA)-- Genetic Algorithm -------------------
GA settings:
Type = binary
Population size = 100
Number of generations = 150
Elitism = 5
Crossover probability = 0.8
Mutation probability = 0.1
GA results:
Iterations = 150
Fitness function value = 0
Solutions =
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 ... x90 x91
[1,] 1 1 0 0 1 0 0 1 1 0 1 0
[2,] 0 1 0 0 1 0 0 0 0 1 1 1
[3,] 1 1 0 0 1 0 0 0 1 0 1 1
[4,] 1 1 0 1 1 0 1 0 0 1 1 1
[5,] 1 1 0 0 1 0 0 0 1 0 1 0
[6,] 0 1 0 0 1 0 0 0 0 1 1 0
[7,] 0 1 0 0 1 0 0 0 0 1 1 1
[8,] 0 1 0 1 1 1 1 0 1 1 1 1
[9,] 1 1 0 0 1 0 0 0 1 0 1 0
[10,] 0 1 0 1 1 1 1 0 1 1 1 0
...
[93,] 0 1 1 1 1 0 1 0 0 1 1 0
[94,] 1 1 0 0 1 0 1 0 0 1 1 0#plot(GA)Aqui devo dizer que fiz diversas simulações de parametros (e.g. maxiter=1500,popSize=10000). Nada funcionou.
É, não saiu como eu esperava. Onde estou errando?
Vamos tentar com uma amostra de 10 municípios
dados_sem_rj <-dados_sem_rj[2:11,]
chromosomes=c(0,1,1,1,0, 0, 0, 0, 0, 0)
GA=ga(type='binary',fitness=fitness,nBits=nrow(dados_sem_rj),maxiter=150,popSize=100,seed=1234,keepBest=TRUE)
summary(GA)-- Genetic Algorithm -------------------
GA settings:
Type = binary
Population size = 100
Number of generations = 150
Elitism = 5
Crossover probability = 0.8
Mutation probability = 0.1
GA results:
Iterations = 150
Fitness function value = 2422761329
Solution =
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
[1,] 0 1 0 0 1 0 0 0 1 0cat(GA@solution%*%dados_sem_rj$area)325176GA_final <-c(0,1,0,0,1,0,0,0,1,0)Parece que o Algoritmo Genético funciona bem com mochilas pequenas, mas não funcionou bem quando tentei com os 91 municipios. Para ter certeza que está funcionando, vamos comparar os resultados do GA com a abordagem clássica com o pacote Adagio.
GA_final <-c(0,1,0,0,1,0,0,0,1,0)
dados_sem_rj[GA_final==1,] Munic area PIB indices incluir
3 Macaé 121699 1298785681 3 0
6 Resende 109934 634577701 6 0
10 Nova Friburgo 93543 489397947 10 0is2 <- knapsack(dados_sem_rj$area, dados_sem_rj$PIB, 403249)
is2$indices[1] 2 5 9Os dois métodos chegaram ao mesmo resultado (Macaé,Resende,Nova Friburgo) com uma amostra pequena, mas eu preciso investigar depois o motivo do GA não funcionar para os 91 municipios do Rio. Por hoje, fico por aqui.