Economia Regional: exemplo de localização Comper

Abstract

We analyse a simple localization problem. Class exercise

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Economia Regional: exemplo de localização Comper. Campo Grande-MS,Brasil: RStudio/Rpubs, 2021. Disponível em http://rpubs.com/amrofi/comper_localization_weber.

1 Introdução

Os primeiros passos são criar ou abrir um diretório de trabalho. Se optar por criar um novo projeto, haverá a possibilidade de criar em uma pasta vazia. Em seguida, sugere-se que coloque os dados nesta pasta, se possível em um arquivo MS Excel e chame a planilha de ‘dados’.Neste caso, a planilha de dados será colocada dentro do código gerado pela função dput(). Os dados básicos foram extraídos de modo parcimonioso, para fins de exercício, a partir do Google Earth (manualmente) e convertido do sistema GMS para o sistema GD pela calculadora do INPE: http://www.dpi.inpe.br/calcula/. O código original veio de notas de aulas de VITON (2014).

2 Dados das coordenadas

Aqui chamaremos os dados embedded.

dados <- structure(list(local = c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J"), 
    scoord = c(-20.4797, -20.4524833333, -20.4916583333, -20.4793916667, -20.4547444444, 
        -20.5157111111, -20.4771222222, -20.4359805556, -20.4656388889, -20.5159), 
    tcoord = c(-54.6105055556, -54.5938805556, -54.6454638889, -54.6407638889, -54.6345361111, 
        -54.6347361111, -54.59175, -54.6048861111, -54.6041, -54.65735), theta = c(1, 
        2, 1, 2, 2, 1, 2, 2, 1, 2), kcost = c(100, 100, 100, 100, 100, 100, 100, 
        100, 100, 100)), row.names = c(NA, -10L), class = c("tbl_df", "tbl", "data.frame"))

O arquivo básico contém 5 variáveis e 10 observações, a saber:

knitr::kable(dados)

local	scoord	tcoord	theta	kcost
A	-20.47970	-54.61051	1	100
B	-20.45248	-54.59388	2	100
C	-20.49166	-54.64546	1	100
D	-20.47939	-54.64076	2	100
E	-20.45474	-54.63454	2	100
F	-20.51571	-54.63474	1	100
G	-20.47712	-54.59175	2	100
H	-20.43598	-54.60489	2	100
I	-20.46564	-54.60410	1	100
J	-20.51590	-54.65735	2	100

em que: local são as lojas identificadas por letras; scoord e tcoord são, respectivamente, latitude e longitude extraídas do Google Earth e convertidas pela calculdadora online do INPE; theta e kcost são as variáveis do modelo, para o peso theta e o custo de transporte associado (shipping rate) \(k\) em unidades monetárias por km e por tonelada. A função a ser minimizada será o custo total de transporte do tipo:

\[ TTC=(k_1.θ_1.d(z_1,f))+(k_2.θ_2.d(z_2,f))+...+(k_n.θ_n.d(z_n,f)) \]

Os thetas são requisitos de insumos (em tons) para produzir 1 ton do produto final n cujo \(θ_n=1\). A medida d é a distância métrica, e \(d(a, b)\) é a distância entre os pontos a e b. Em um cenário padrão, todos os \(k\) seriam iguais.

Neste exemplo, para os dados acima dispostos, a primeira observação é o produto (mercado) em A. A otimização fornecerá o ponto central que minimiza os custos entre estes locais.

O modelo considerará a função distância edist, para uma matriz de coordenadas, e uma função para o custo - cost.

# NAO MEXER NAS FUNCOES
edist <- function(x, cmx) {
    # cmx is an nx2 coordinate matrix
    mx <- matrix(rep(x, each = nrow(cmx)), nrow = nrow(cmx))
    res <- sqrt(rowSums((cmx - mx)^2))
    res
}

cost <- function(x, data) {
    d <- edist(x, data[, c("scoord", "tcoord")])
    c <- sum(d * data[, "theta"] * data[, "kcost"])
    c
}

A função weber2 fará a otimização a partir de um ponto generico de coordenadas (0,0) (o leitor pode alterar se desejar, chamando para um ponto x=c(a,b) inicial). A função wplot fará um plot da solução.

weber2 <- function(start = c(0, 0), data = data) {
    z <- optim(start, cost, gr = NULL, data)
    weber_res <<- z  # save result as a global
    wplot(data, z)
    cat("\nProblem setting (locations can be inputs or final output):\n")
    for (i in 1:nrow(data)) {
        pr <- sprintf("Location %i : (%5.3f,%5.3f)\n   unit requirement %4.2f (tons); trans cost: %5.3f (/ton-mi)", 
            i, data[i, "scoord"], data[i, "tcoord"], data[i, "theta"], data[i, "kcost"])
        cat(pr, "\n")
    }
    # format = f avoids E notation for tiny values, which can arise for problems on a
    # line
    cat("\nSolution (see plot window):\nFactory located at: ", formatC(z$par, digits = 4, 
        format = "f"), "\n")
    cat("Total trans cost: ", formatC(z$value, digits = 5), "\n  ")
}
wplot <- function(data, res = NULL) {
    m <- data[, c("scoord", "tcoord")]
    m <- rbind(m, m[1, ])
    plot(m, type = "o")
    # in the next line, pch=19 is a filled circle
    if (!is.null(res)) {
        points(t(res$par), pch = 19)
    }
}

3 Execução

3.1 Com pesos diferentes

Neste cenário, nossos dados são como obtidos na seção anterior, com theta variando 1 e 2 os diferentes pontos.

weber2(start = c(0, 0), data = dados)

Problem setting (locations can be inputs or final output):
Location 1 : (-20.480,-54.611)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 2 : (-20.452,-54.594)
   unit requirement 2.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 3 : (-20.492,-54.645)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 4 : (-20.479,-54.641)
   unit requirement 2.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 5 : (-20.455,-54.635)
   unit requirement 2.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 6 : (-20.516,-54.635)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 7 : (-20.477,-54.592)
   unit requirement 2.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 8 : (-20.436,-54.605)
   unit requirement 2.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 9 : (-20.466,-54.604)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 10 : (-20.516,-54.657)
   unit requirement 2.00 (tons); trans cost: 100.000 (/ton-mi) 

Solution (see plot window):
Factory located at:  -20.4711 -54.6177 
Total trans cost:  51.229 
  

A solução é portanto o par para latitude -20.4711402 e longitude -54.6177371.

Vamos agora plotar no mapa. Primeiro extrairemos um mapa com os pontos dos mercados.

attach(dados)
# hiper compers
require(RgoogleMaps)
mt = GetMapTiles(latR = scoord, lonR = tcoord, verbose = 1, zoom = 10)

tileserver: http://mt1.google.com/vt/lyrs=m 
nTiles= 1 1 , center= -20.477 -54.622 
1 tiles to download 
tile image format: .png 
NOT downloading existing file 10_356_571.png
sleptTotal= 0 

PlotOnMapTiles(mt, lat = scoord, lon = tcoord, pch = 20, col = c("red", "blue", "green"), 
    cex = 2)

Depois ilustramos como alterar a cor dos pontos.

PlotOnMapTiles(mt, lat = scoord, lon = tcoord, pch = 10, col = c("red"), cex = 2)

# resultado weber_res$par[1] [1] -20.47114 weber_res$par[2] [1] -54.61774

latw = c(weber_res$par[1])
lonw = c(weber_res$par[2])
local_final = rbind(dados$local, "CT")
places = c("HComper", "HComper", "HComper", "HComper", "HComper", "HComper", "HComper", 
    "HComper", "HComper", "HComper", "CT")
mt = GetMapTiles(latR = scoord, lonR = tcoord, verbose = 1, zoom = 12)

tileserver: http://mt1.google.com/vt/lyrs=m 
nTiles= 1 2 , center= -20.477 -54.622 
2 tiles to download 
tile image format: .png 
NOT downloading existing file 12_1426_2285.png
NOT downloading existing file 12_1426_2286.png
sleptTotal= 0 

Por último, colocamos os pontos dos mercados junto com o do CT encontrado.

PlotOnMapTiles(mt, lat = c(scoord, latw), lon = c(tcoord, lonw), pch = c(20, 20, 
    20, 20, 20, 20, 20, 20, 20, 20, 10), col = c("red", "red", "red", "red", "red", 
    "red", "red", "red", "red", "red", "blue"), cex = 3)

3.2 Colocando theta =1 para todos

# colocando theta =1 para todos
dados1 <- structure(list(local = c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J"), 
    scoord = c(-20.4797, -20.4524833333, -20.4916583333, -20.4793916667, -20.4547444444, 
        -20.5157111111, -20.4771222222, -20.4359805556, -20.4656388889, -20.5159), 
    tcoord = c(-54.6105055556, -54.5938805556, -54.6454638889, -54.6407638889, -54.6345361111, 
        -54.6347361111, -54.59175, -54.6048861111, -54.6041, -54.65735), theta = c(1, 
        1, 1, 1, 1, 1, 1, 1, 1, 1), kcost = c(100, 100, 100, 100, 100, 100, 100, 
        100, 100, 100)), row.names = c(NA, -10L), class = c("tbl_df", "tbl", "data.frame"))

Agora os dados tem theta =1 para todos:

knitr::kable(dados1)

local	scoord	tcoord	theta	kcost
A	-20.47970	-54.61051	1	100
B	-20.45248	-54.59388	1	100
C	-20.49166	-54.64546	1	100
D	-20.47939	-54.64076	1	100
E	-20.45474	-54.63454	1	100
F	-20.51571	-54.63474	1	100
G	-20.47712	-54.59175	1	100
H	-20.43598	-54.60489	1	100
I	-20.46564	-54.60410	1	100
J	-20.51590	-54.65735	1	100

weber2(start = c(0, 0), data = dados1)

Problem setting (locations can be inputs or final output):
Location 1 : (-20.480,-54.611)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 2 : (-20.452,-54.594)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 3 : (-20.492,-54.645)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 4 : (-20.479,-54.641)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 5 : (-20.455,-54.635)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 6 : (-20.516,-54.635)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 7 : (-20.477,-54.592)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 8 : (-20.436,-54.605)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 9 : (-20.466,-54.604)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 
Location 10 : (-20.516,-54.657)
   unit requirement 1.00 (tons); trans cost: 100.000 (/ton-mi) 

Solution (see plot window):
Factory located at:  -20.4749 -54.6167 
Total trans cost:  30.822 
  

Agora os resultados mudaram para o CT na latitude -20.4748517 e longitude -54.6167492.
Vamos agora plotar no mapa. Primeiro extrairemos um mapa com os pontos dos mercados. Depois ilustramos como alterar a cor dos pontos. Por último, colocamos os pontos dos mercados junto com o do CT encontrado.

attach(dados1)
# hiper compers
require(RgoogleMaps)
mt = GetMapTiles(latR = scoord, lonR = tcoord, verbose = 1, zoom = 10)

tileserver: http://mt1.google.com/vt/lyrs=m 
nTiles= 1 1 , center= -20.477 -54.622 
1 tiles to download 
tile image format: .png 
NOT downloading existing file 10_356_571.png
sleptTotal= 0 

PlotOnMapTiles(mt, lat = scoord, lon = tcoord, pch = 20, col = c("red", "blue", "green"), 
    cex = 2)

PlotOnMapTiles(mt, lat = scoord, lon = tcoord, pch = 10, col = c("red"), cex = 2)

# resultado -20.4749 -54.6167
latw = c(weber_res$par[1])
lonw = c(weber_res$par[2])
local_final = rbind(dados1$local, "CT")
places = c("HComper", "HComper", "HComper", "HComper", "HComper", "HComper", "HComper", 
    "HComper", "HComper", "HComper", "CT")
mt = GetMapTiles(latR = scoord, lonR = tcoord, verbose = 1, zoom = 12)

tileserver: http://mt1.google.com/vt/lyrs=m 
nTiles= 1 2 , center= -20.477 -54.622 
2 tiles to download 
tile image format: .png 
NOT downloading existing file 12_1426_2285.png
NOT downloading existing file 12_1426_2286.png
sleptTotal= 0 

PlotOnMapTiles(mt, lat = c(scoord, latw), lon = c(tcoord, lonw), pch = c(20, 20, 
    20, 20, 20, 20, 20, 20, 20, 20, 10), col = c("red", "red", "red", "red", "red", 
    "red", "red", "red", "red", "red", "blue"), cex = 3)

4 Usar Google maps para melhor visualização do CT

Neste caso, usaremos o resultado do primeiro exemplo, com pesos diferentes, e procuramos nas coordenadas c(-20.47114,-54.61774) no Google Maps.

Localizacao do CT - Cenario 1. Fonte: Elaboração própria com <https://www.google.com.br/maps/>.

Referências

VITON, Philip A. Using R to Solve Weber Problems. (CRPLAN 6600 slides) Ohio State University, 2014.