São Paulo e o problema da mochila

Por Julio 12/08/2017

São Paulo é a minha cidade preferida. Não só porque moro aqui, mas também porque é uma cidade cheia de diversidade, boa gastronomia e oportunidades. Para sentir um pouco dessa vibe, recomendo passear na avenida Paulista aos domingos. É sensacional!

Mas a cidade da diversidade só é o que é porque temos muita, muita gente nela. O município tem 12 milhões de habitantes. Esse número é tão grande que temos um paulistano para cada 17 brasileiros! Se São Paulo fosse um país, seria o 77 do mundo, ganhando de países como a Bélgica, Grécia, Portugal, Bolívia e muitas outras.

Outro dia eu estava pensando na seguinte problemática: qual é a área do Brasil ocupada pela população de São Paulo? Ou seja, se pegarmos os municípios com grandes áreas, quanto do país conseguiríamos preencher com 12 milhões de habitantes?

O interessante é que essa questão recai exatamente no problema da mochila, que é um famoso desafio de programação inteira. Depois de estudar profundamente no wikipedia (😂), vi que o problema não é tão trivial como parece.

O problema da mochila

Considere o seguinte contexto: você tem uma mochila com capacidade de 15kg e precisa carregar a combinação de itens com maior valor, com cada item possuindo valores e pesos diferentes.

Knapsack problem. Retirado do [Wikipedia](https://en.wikipedia.org/wiki/Knapsack_problem#/media/File:Knapsack.svg).

Figure 1: Knapsack problem. Retirado do Wikipedia.

Outra forma de pensar nesse problema é com um cardápio de restaurante:

XKCD sobre o knapsack problem.

Figure 2: XKCD sobre o knapsack problem.

Em linguagem matemática, o que temos é a task:

\[ \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} \]

No nosso caso essas letras significam isso aqui:

  • \(n\) é o número de municípios no Brasil (5570).
  • \(v_i\) é a área do município \(i\).
  • \(w_i\) é a população do município \(i\).
  • \(W\) é a população de São Paulo (12 milhões).
  • \(x=(x_1,\dots,x_n)^\top\) é o vetor que seleciona os municípios. Se o município \(i\) faz parte da solução \(x_i=1\) e, caso contrário, \(x_i=0\).

Ou seja, queremos escolher municípios para colocar na mochila tentando maximizar a área, mas o máximo de população que podemos contemplar é 12 milhões.

O problema da mochila é muito interessante pois trata-se de um problema NP-difícil, ou seja, não existe um algoritmo de polinomial capaz de resolvê-lo. Se \(w_i > 0, \forall i\in1,\dots,n\) então a solução pode ser encontrada com um algoritmo pseudo-polinomial.

Forma ad-hoc

Se \(x_i\) pudesse assumir valores entre zero e um (ou seja, se pudéssemos selecionar apenas pedaços de municípios), a solução seria trivial. Bastaria colocar os municípios em ordem decrescente pela razão \(v_i/w_i\) e escolher os municípios ou parte deles até obter \(W\).

Isso indica uma forma sub-ótima de resolver o problema. Chamamos essa solução de ad-hoc. A solução é encontrada assim:

  1. Colocar os municípios em ordem decrescente pela razão \(v_i/w_i\),
  2. Escolher os municípios de maior razão até que a população do próximo município estoure \(W\).
  3. Escolher outros municípios com maior razão na ordem até não ser possível incluir mais nenhum município.

Solução ótima

A solução ótima pode ser encontrada usando a função mknapsack() do pacote adagio. Por exemplo, considere os vetores de pesos w, valores p e máximo cap abaixo.

p <- c(15, 100, 90, 60, 40, 15, 10,  1)
w <- c( 2,  20, 20, 30, 40, 30, 60, 10)
cap <- 102

O vetor-solução é dado por

is <- adagio::mknapsack(p, w, cap)
is$ksack
## [1] 1 1 1 1 0 1 0 0

Dados

As áreas e estimativas das populações dos municípios do Brasil em 2016 foram obtidas do IBGE. A leitura é realizada usando pacotes do tidyverse.

Pacotes:

library(tidyverse)
library(janitor)
library(readxl)
library(httr)

Áreas:

url_area <- paste0(
  'ftp://geoftp.ibge.gov.br',
  '/organizacao_do_territorio/estrutura_territorial/areas_territoriais/2016/',
  'AR_BR_RG_UF_MUN_2016.xls'
) %>% GET(write_disk('area.xls', overwrite = TRUE))

area <- read_excel('area.xls') %>%
  clean_names() %>%
  filter(!is.na(id)) %>% 
  select(uf = nm_uf_sigla, cd_gcmun, muni = nm_mun_2016, area = ar_mun_2016)

População:

url_pop <- paste0(
  'ftp://ftp.ibge.gov.br',
  '/Estimativas_de_Populacao/Estimativas_2016/',
  'estimativa_TCU_2016_20170614.xls'
) %>% GET(write_disk('pop.xls', overwrite = TRUE))

loc <- locale(grouping_mark = '.', decimal_mark = ',')
pop <- read_excel('pop.xls', skip = 1, sheet = 2) %>%
  clean_names() %>%
  filter(!is.na(cod_munic)) %>%
  mutate(pop = parse_number(`população_estimada`, locale = loc),
         cd_gcmun = paste0(cod_uf, cod_munic)) %>%
  select(cd_gcmun, pop)

Para completar, temos os shapefiles dos municípios. Os originais foram baixados nesse link, mas fiz um código para reduzir a resolução e preferi omitir do post. Se quiser ver, entre no repositório que gera esse site. Para ler os shapes usamos o pacote sf:

No final, vamos dar um join das tabelas:

dados <- area %>% 
  left_join(pop, 'cd_gcmun') %>%
  mutate(razao = area / pop) %>%
  filter(!is.na(pop))

Resultados

A solução ad-hoc e ótima são computadas com esse código:

d_solucao <- dados %>% 
  arrange(desc(razao)) %>% # ordena para solucao adhoc funcionar
  mutate(area2 = as.integer(area * 1000), # necessario para mknapsack funcionar
         s_knapsack = adagio::mknapsack(area2, pop, max(pop))$ksack,
         acu = cumsum(pop),
         s_adhoc0 = if_else(acu < max(pop), 1, 0),
         s_adhoc = s_adhoc0) 

Agora, vamos melhorar a solução ad-hoc incluindo os melhores municípios.

id_melhor <- 0
pop_faltam <- with(d_solucao, max(pop) - sum(s_adhoc0 * pop))
while (is.na(id_melhor)) {
  # pega id do melhor municipio a ser incluido
  id_melhor <- with(d_solucao, which(pop <= pop_faltam & s_adhoc == 0)[1])
  if (is.na(id_melhor)) {
    d_solucao$s_adhoc[id_melhor] <- 1
    pop_faltam <- with(d_solucao, max(pop) - sum(s_adhoc * pop))
  }
}

A Tabela 1 mostra os municípios que foram classificados diferentemente nos dois métodos. Note que a solução ótima trocou apenas um município da solução adhoc (Nova Aurora - GO) pelo município de Monte Alegre de Minas - MG.

d_solucao %>% 
  filter(s_adhoc != s_knapsack) %>% 
  select(uf, muni, area, pop, s_adhoc, s_knapsack) %>% 
  knitr::kable(caption = 'Municípios diferentes nas duas soluções.')
Table 1: Municípios diferentes nas duas soluções.
uf muni area pop s_adhoc s_knapsack
GO NOVA AURORA 302.655 2194 1 0
MG MONTE ALEGRE DE MINAS 2595.957 20979 0 1

A Tabela 2 mostra a diferença dos resultados dos dois métodos. A solução ótima fica com apenas 92 pessoas a menos que São Paulo.

Table 2: Diferença dos resultados.
Método Área total População total Diferença para sp
adhoc 5576985.742 12019298 18877
knapsack 5579279.044 12038083 92
São Paulo - 12038175 0

Mapa final

Visualmente, a solução ótima e a solução adhoc são idênticas. Por isso vou mostrar apenas como fica o mapa para a solução ótima.

O resultado aparece na Figura 3. É realmente impressionante ver que aquela regiãozinha vermelha tem a mesma população que toda a região azul do mapa.

Mapa do Brasil final.

Figure 3: Mapa do Brasil final.

É isso! Happy coding ;)

comments powered by Disqus