Proyecto Final: Simulacion 1

Paula Cazali

14000060

Cargamos las librerias necesarias para el proyecto:

library(readr)
library(tidyverse)
library(httr)
library(jsonlite)
library(ggmap)
library(gepaf)
library(dplyr)

Se usara un API key y y el API de google para obtener las direcciones desde cada punto.

api_key <- "api_key_goes_here"
register_google(key=api_key)

api_key <- "&key=api_key_goes_here"
url_part_1 <- "https://maps.googleapis.com/maps/api/distancematrix/json?units=metric&origins="
url_part_2 <- "&destinations="

Cargamos la informacion de los hospitales con un .csv y creamos un dataset con las coordenadas unidas

# hosp_coordinates <- read_csv("Universidad Galileo/Maestria/Simulacion 1/Proyecto Final/Proyecto Final/hosp_coordinates.csv")
hosp_coordinates <- read_csv("hosp_coordinates.csv")

Parsed with column specification:
cols(
  id = col_double(),
  lat = col_double(),
  lon = col_double(),
  hospital = col_character()
)

hosp_coordinates <- hosp_coordinates %>%
  unite(coordinate, c(lat,lon), sep = ",")
hosp_coordinates

Funcion que calcula la distancia desde un origen a multiples destinos. Devuelve un vector con las distancias

get_distance <- function(origin_id, coordinates, distance = T){
  destination <- coordinates[!coordinates$id == origin_id,]
  origin <- coordinates[coordinates$id == origin_id,]
  url <- paste0(url_part_1,
                origin$coordinate,
                url_part_2,
                paste0(destination$coordinate,collapse = "|"),
                api_key,collapse = "")
  
  response <- GET(url)
  resp_json <- fromJSON(content(response, as = "text"))
  resp_df <- resp_json[[3]][[1]][[1]]
  if(distance) return(resp_df[[1]][2]$value)
  else return(resp_df[[2]][2]$value)
}

En este for se hara una matriz que tiene las distancias desde cada ciudad hasta cada ciudad, se usa la funcion get_distance()

dist_matrix <- c()
for(hosp in 1:nrow(hosp_coordinates)){
  hosp_coordinates$id
  distance_vec <- get_distance(hosp, hosp_coordinates)
  distance_vec <- append(distance_vec, 0, after = hosp - 1)
  print(distance_vec)
  dist_matrix <- c(dist_matrix, distance_vec)
  #print(hosp)
}

 [1]     0  5543  7458  9543  9776  9682  4553  4553  3344  4908 11895  7853 10778
 [1]  5211     0  4661  6188  6421  5502  5237  5237  4711  5593 11143  8564  7069
 [1] 5777 4213    0 2311 2544 5488 7660 7660 4901 8016 9013 8384 5399
 [1]  7022  5275  2374     0   233  4016  8906  8906  6147  9261  9532 12396  3915
 [1]  7549  5056  2043  1017     0  3783  9432  9432  6674  9788  9727 12591  3791
 [1]  8936  5016  4412  4405  4638     0  9012  9012  7863  9368 13736 11651  3323
 [1]  5102  4259  7085  9164  9398  8398     0     0  4602   721 11516  9512  9493
 [1]  5102  4259  7085  9164  9398  8398     0     0  4602   721 11516  9512  9493
 [1]  3081  5555  5785  7087 13063  9129  6151  6151     0  6507 10439  4772 15958
 [1]  4816  3973  6799  8878  9111  8112  1198  1198  4316     0 11230  9226  9207
 [1]  9098 10742  9007  9606  9512 12983 10981 10981  7604 11336     0  9235 12407
 [1]  7101 10738  8965 13441 13347 12940 10977 10977  5710 11333 10723     0 16242
 [1]  9563  5625  5039  4265  4498  3137  9185  9185  8490  9541 12150 15014     0

dist_matrix <- matrix(dist_matrix, ncol = nrow(hosp_coordinates), byrow = TRUE)

Algoritmo Genetico para obtener la ruta mas corta

Definicion de funciones

# distancia total de toda la ruta
distance_route <- function(route,distance){
  sum_distance<-0
  for(i in 1:length(route[-1])){
    out <- distance[route[i],route[i+1] ]
    sum_distance <- sum_distance + out
  }
  return(sum_distance)
}

# genera una ruta aleatoria
sample_route <- function(..., cities){
  c(1,sample(2:cities),1) %>% return()
}

# Selecciona los padres en base a la distancia (roulette)
select_mating_parents <- function(...,pop_size, roullete, population){
  sum_fit_p <-
    sample(1:sum(roullete$rank), size = 1 )
  pindex <-
    roullete %>% 
    filter(cumsum_rank < sum_fit_p) %>%
    nrow()
  p1 <- roullete[pindex + 1,] %>% pull(parent)
  sum_fit_p <-
    sample(1:sum(roullete$rank),size = 1 )
  pindex <-
    roullete %>% 
    filter(cumsum_rank < sum_fit_p) %>%
    nrow()
  p2 <- roullete[pindex + 1,] %>% pull(parent)
  return(population[c(p1,p2)])
}

# Funcion para hacer el crossover con dos padres, se usa material genetico de uno y del otro.
position_based_crossover <- function(parents, n = 3){
  p1 <- parents[[1]]
  p2 <- parents[[2]]
  length_parent <- length(p1)
  p1 <- p1[2:(length_parent-1)]
  p2 <- p2[2:(length_parent-1)]
  index <- sample(1:(length_parent-2), size = n, replace = FALSE)
  child1 <- p1
  fill_values <- setdiff(p2, p1[index])
  child1[-index] <- fill_values
  return(c(1, child1, 1)) 
}

#Cambia de posicion dos ciudades, una por la otra.
swap <- function(vec, i, j){
  out <- vec
  out[i] <- vec[j]
  out[j] <- vec[i]
  return(out)
}

# Funcion para mutar los hijos
mutate_child <- function(child, rate = 0.03){
  if(runif(1) < rate){
    index <- 2:(length(child) - 1)
    i <- sample(index, size = 1)
    index <- setdiff(index,i)
    j <- sample(index, size = 1)
    child <- swap(child, i, j)
    #print(1)
  }
  return(child)
}

Obteniendo la ruta mas corta

population <- lapply(1:100, sample_route, cities = nrow(hosp_coordinates))

# FOR PARA OBTENER LA RUTA MAS CORTA.
for(i in 1:100){
  pop_fitness <- lapply(population, distance_route, distance = dist_matrix)
  
  roullete <-
    tibble(parent = 1:100, fitness = pop_fitness %>% unlist()) %>%
    arrange(desc(fitness))
  
  roullete$rank <- 1:nrow(roullete)
  
  roullete <-
    roullete %>% 
    mutate(cumsum_rank = cumsum(rank))
  
  mating_parents <-
    lapply(1:100, select_mating_parents,
           pop_size = 100,
           roullete = roullete, 
           population = population)
  
  children <- lapply(mating_parents, position_based_crossover, n = 5)
  
  population <- lapply(children, mutate_child, rate = 0.1)
  #print(pop_fitness[[which.min(pop_fitness %>% unlist())]])
}

min_route <- population[which.min(pop_fitness %>% unlist())]
min_route[[1]]

 [1]  1 12 11  3  4  5 13  6  2  8  7 10  9  1

Graficar la ruta mas corta

url_route1 <- "https://maps.googleapis.com/maps/api/directions/json?origin="
url_route2 <- "&destination="

# OBTENER LOS POLYLINES
poly_lines <- data.frame()

for(hosp in 1:length(min_route[[1]])){
  if(hosp == length(min_route[[1]])){
    origin <- hosp_coordinates[min_route[[1]][hosp],]
    destination <- hosp_coordinates[min_route[[1]][1],]
    #print(hosp)
  }else{
    #print(hosp)
    origin <- hosp_coordinates[min_route[[1]][hosp],]
    destination <- hosp_coordinates[min_route[[1]][hosp+1],]
  }
  url <- paste0(url_route1,
                origin$coordinate,
                url_route2,
                paste0(destination$coordinate),
                api_key,collapse = "")
  
  response <- GET(url)
  resp_json <- fromJSON(content(response, as = "text"))
  ruta <- resp_json$routes$overview_polyline$points
  ruta <- decodePolyline(ruta)
  aux_df <- ruta
  poly_lines <- rbind(poly_lines, aux_df)
}

Graficar la ruta mas corta

qmap('Estacion Tivoli',zoom = 13, maptype="roadmap") +
  geom_point(aes(x = lon, y = lat), 
             data = poly_lines, 
             colour = "red" )

Source : https://maps.googleapis.com/maps/api/staticmap?center=Estacion%20Tivoli&zoom=13&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx
Source : https://maps.googleapis.com/maps/api/geocode/json?address=Estacion+Tivoli&key=xxx

Algoritmo Simulated Annealing

# Generamos una ruta al azar 
simulated_route <- sample_route(cities = nrow(hosp_coordinates))

# Distancia que hay que recorrer en esa ruta
dist_actual <- distance_route(simulated_route, dist_matrix)

temperature <- 20
delta_temp <- -0.01

for (i in seq(from = temperature, to = 0, by = delta_temp)){
  dist_anterior <- dist_actual
  # quitaremos el punto que esta fijo en la ruta: el Aeropuerto (punto de inicio y fin)
  simulated_route <- simulated_route[2:(length(simulated_route)-1)]
  #cambio en las ciudades
  swap_index <- sample(1:(length(simulated_route)), 2, replace = FALSE)
  simulated_route_2 <- c(1,swap(simulated_route, swap_index[1],swap_index[2]),1)
  dist_actual <- distance_route(simulated_route_2, dist_matrix)
  
  delta_dist <- dist_anterior - dist_actual
  if(delta_dist > 0){
    simulated_route <- simulated_route_2
    #print("cambio")
  } else if(exp(delta_dist/temperature) > runif(1, 0, 1)){
    simulated_route <- simulated_route_2
    #print("probabilidad")
  }else {
    simulated_route <- c(1, simulated_route, 1)
    #print("No cambio")
  } 
}

distancia_obtenida <- distance_route(simulated_route, dist_matrix)
distancia_obtenida

[1] 66264

Obtener los polylines

poly_lines_SA <- data.frame()

for(hosp in 1:length(simulated_route)){
  if(hosp == length(simulated_route)){
    origin <- hosp_coordinates[simulated_route[hosp],]
    destination <- hosp_coordinates[simulated_route[1],]
    #print(hosp)
  }else{
    #print(hosp)
    origin <- hosp_coordinates[simulated_route[hosp],]
    destination <- hosp_coordinates[simulated_route[hosp+1],]
  }
  url <- paste0(url_route1,
                origin$coordinate,
                url_route2,
                paste0(destination$coordinate),
                api_key,collapse = "")
  
  response <- GET(url)
  resp_json <- fromJSON(content(response, as = "text"))
  ruta_SA <- resp_json$routes$overview_polyline$points
  ruta_SA <- decodePolyline(ruta_SA)
  aux_df <- ruta_SA
  poly_lines_SA <- rbind(poly_lines_SA, aux_df)
}

Graficar la ruta encontrada por el algoritmo Simulated Annealing

qmap('Estacion Tivoli',zoom = 13, maptype="roadmap") +
  geom_point(aes(x = lon, y = lat), 
             data = poly_lines_SA, 
             colour = "blue" )

Source : https://maps.googleapis.com/maps/api/staticmap?center=Estacion%20Tivoli&zoom=13&size=640x640&scale=2&maptype=roadmap&language=en-EN&key=xxx
Source : https://maps.googleapis.com/maps/api/geocode/json?address=Estacion+Tivoli&key=xxx