R/assign_traffic.R

In cppRouting: Algorithms for Routing and Solving the Traffic Assignment Problem

#' Algorithms for solving the Traffic Assignment Problem (TAP).
#' @description Estimation of the User Equilibrium (UE)
#' @param Graph  An object generated by \link{makegraph} function.
#' @param from A vector of origins
#' @param to A vector of destinations.
#' @param demand A vector describing the flow between each origin-destination pair.
#' @param algorithm character. \code{msa}, \code{fw}, \code{cfw}, \code{bfw} or \code{dial}. Default to \code{bfw}. See details.
#' @param max_gap Numeric. Relative gap to achieve. Default to 0.001.
#' @param max_it Numeric. Maximum number of iterations. Default to \code{.Machine$integer.max}
#' @param aon_method Character.\code{d}, \code{bi}, \code{nba}, \code{cphast} or \code{cbi}. Default to \code{bi}. See details.
#' @param constant numeric. Constant to maintain the heuristic function admissible in NBA* algorithm. Default to 1, when cost is expressed in the same unit than coordinates. See details
#' @param dial_params List. Named list of hyperparameters for \code{dial} algorithm. See details.
#' @param verbose Logical. If \code{TRUE} (default), progression is displayed.
#' @return A \code{list} containing : \itemize{
#'  \item The relative gap achieved
#'  \item Number of iteration
#'  \item A data.frame containing edges attributes, including equilibrated flows, new costs and free-flow travel times.
#'
#'  }
#' @note \code{from}, \code{to} and \code{demand} must be the same length.
#' \code{alpha}, \code{beta} and \code{capacity} must be filled in during network construction. See \link{makegraph}.
#'
#'
#' @details
#' The most well-known assumptions in traffic assignment models are the ones following Wardrop's first principle.
#' Traffic assignment models are used to estimate the traffic flows on a network. These models take as input a matrix of flows that indicate the volume of traffic between origin and destination (O-D) pairs.
#' Unlike All-or-Nothing assignment (see \link{get_aon}), edge congestion is modeled through the \strong{Volume Decay Function (VDF)}.
#' The Volume Decay Function used is the most popular in literature, from the Bureau of Public Roads :
#'
#' \strong{t = t0 * (1 + a * (V/C)^b) }
#' with t = actual travel time (minutes),
#' t0 = free-flow travel time (minutes),
#' a = alpha parameter (unitless),
#' b = beta parameter (unitless),
#' V = volume or flow (veh/hour)
#' C = edge capacity (veh/hour)
#'
#' Traffic Assignment Problem is a convex problem and solving algorithms can be divided into two categories : \itemize{
#' \item link-based : \strong{Method of Successive Average} (\code{msa}) and \strong{Frank-Wolfe variants} (normal : \code{fw}, conjugate : \code{cfw} and bi-conjugate : \code{bfw}).
#' These algorithms uses the descent direction given by AON assignment at each iteration, all links are updated at the same time.
#' \item bush-based : \strong{Algorithm-B} (\code{dial})
#' The problem is decomposed into sub-problems, corresponding to each origin of the OD matrix, that operate on acyclic sub-networks of the original transportation network, called bushes.
#' Link flows are shifted from the longest path to the shortest path recursively within each bush using Newton method.
#' }
#'
#' Link-based algorithms are historically the first algorithms developed for solving the traffic assignment problem. It require low memory and are known to tail in the vicinity of the optimum and usually cannot be used to achieve highly precise solutions.
#' Algorithm B is more recent, and is better suited for achieve the highest precise solution. However, it require more memory and can be time-consuming according the network size and OD matrix size.
#' In \code{cppRouting}, the implementation of algorithm-B allow "batching", i.e. bushes are temporarily stored on disk if memory limit, defined by the user, is exceeded.
#' Please see the package website for practical example and deeper explanations about algorithms. (\url{https://github.com/vlarmet/cppRouting/blob/master/README.md})
#'
#' Convergence criterion can be set by the user using max_gap argument, it is the relative gap which can be written as :
#' \strong{abs(TSTT/SPTT - 1)}
#' with TSTT (Total System Travel Time) = sum(flow * cost),
#' SPTT (Shortest Path Travel Time) = sum(aon * cost)
#'
#' Especially for link-based algorithms (msa, *fw), the larger part of computation time rely on AON assignment. So, choosing the right AON algorithm is crucial for fast execution time.
#' Contracting the network on-the-fly before AON computing can be faster for large network and/or large OD matrix.
#'
#' AON algorithms are : \itemize{
#' \item \code{bi} : bidirectional Dijkstra algorithm
#' \item \code{nba} : bidirectional A* algorithm, nodes coordinates and constant parameter are needed
#' \item \code{d} : Dijkstra algorithm
#' \item \code{cbi} : contraction hierarchies + bidirectional search
#' \item \code{cphast} : contraction hierarchies + phast algorithm
#' }
#' These AON algorithm can be decomposed into two families, depending the sparsity of origin-destination matrix : \itemize{
#' \item recursive pairwise : \code{bi}, \code{nba} and \code{cbi}. Optimal for high sparsity. One-to-one algorithm is called N times, with N being the length of from.
#' \item recursive one-to-many : \code{d} and \code{cphast}. Optimal for dense matrix. One-to-many algorithm is called N times, with N being the number of unique from (or to) nodes
#' }
#'
#' For large instance, it may be appropriate to test different \code{aon_method} for few iterations and choose the fastest one for the final estimation.
#'
#' Hyperparameters for algorithm-b are : \itemize{
#' \item \code{inneriter} : number of time bushes are equilibrated within each iteration. Default to 20
#' \item \code{max_tol} : numerical tolerance. Flow is set to 0 if less than max_tol. Since flow shifting consist of iteratively adding or substracting double types, numerical error can occur and stop convergence.
#' Default to 1e-11.
#' \item \code{tmp_path} : Path for storing bushes during algorithm-B execution. Default using \code{tempdir()}
#' \item \code{max_mem} : Maximum amount of RAM used by algorithm-B in gigabytes. Default to 8.
#'
#' }
#'
#'
#' In New Bidirectional A star algorithm, euclidean distance is used as heuristic function.
#' To understand the importance of constant parameter, see the package description : \url{https://github.com/vlarmet/cppRouting/blob/master/README.md}
#' All algorithms are partly multithreaded (AON assignment).
#' @examples
#' #Choose number of cores used by cppRouting
#' RcppParallel::setThreadOptions(numThreads = 1)
#'
#' #Data describing edges of the graph
#' edges<-data.frame(from_vertex=c(0,0,1,1,2,2,3,4,4),
#'                   to_vertex=c(1,3,2,4,4,5,1,3,5),
#'                   cost=c(9,2,11,3,5,12,4,1,6))
#'
#' # Origin-destination trips
#' trips <- data.frame(from = c(0,0,0,0,1,1,1,1,2,2,2,3,3,4,5,5,5,5,5),
#'                     to = c(1,2,5,3,2,5,2,4,2,5,2,3,5,2,0,0,3,5,1),
#'                     flow = c(10,30,15,5,5,2,3,6,4,15,20,2,3,6,2,1,4,5,3))
#'
#' #Construct graph
#' graph <- makegraph(edges,directed=TRUE, alpha = 0.15, beta = 4, capacity = 5)
#'
#'
#' # Solve traffic assignment problem
#' ## using Bi-conjugate Frank-Wolfe algorithm
#' traffic <- assign_traffic(Graph=graph,
#'                           from=trips$from, to=trips$to, demand = trips$flow,
#'                           algorithm = "bfw")
#' print(traffic$data)
#'
#' ## using algorithm-B
#' traffic2 <- assign_traffic(Graph=graph,
#'                            from=trips$from, to=trips$to, demand = trips$flow,
#'                            algorithm = "dial")
#' print(traffic2$data)
#' @references Wardrop, J. G. (1952). "Some Theoretical Aspects of Road Traffic Research".
#' @references M. Fukushima (1984). "A modified Frank-Wolfe algorithm for solving the traffic assignment problem".
#' @references R. B. Dial (2006). "A path-based user-equilibrium traffic assignment algorithm that obviates path storage and enumeration".
#' @references M. Mitradjieva, P. O. Lindberg (2012).  "The Stiff Is Moving — Conjugate Direction Frank-Wolfe Methods with Applications to Traffic Assignment".


assign_traffic <- function(Graph, from, to, demand, algorithm = "bfw", max_gap = 0.001, max_it = .Machine$integer.max,
                           aon_method = "bi", constant = 1, dial_params = NULL, verbose = TRUE){
  if (length(from)!=length(to) | length(from) != length(demand)) stop("From, to and demand have not the same length")
  demand <- as.numeric(demand)
  if (any(is.na(data.frame(from,to, demand)))) stop("NAs are not allowed in origin/destination trips")
  ind <- which(demand != 0)
  from <- from[ind]
  to <- to[ind]
  demand <- demand[ind]

  if (length(Graph) != 5) stop("Graph object must be generated by makegraph() function")
  if (is.null(Graph$attrib$alpha) | is.null(Graph$attrib$beta)| is.null(Graph$attrib$cap)) stop("alpha, beta and capacity must be defined during graph construction")


  from<-as.character(from)
  to<-as.character(to)

  allnodes<-c(from,to)
  if (sum(allnodes %in% Graph$dict$ref) < length(allnodes)) stop("Some nodes are not in the graph")

  from_id <- Graph$dict$id[match(from,Graph$dict$ref)]
  to_id <- Graph$dict$id[match(to,Graph$dict$ref)]

  if (!algorithm %in% c("msa", "fw", "cfw", "bfw", "dial")) stop("algorithm should be 'msa', 'fw', 'cfw', 'bfw' or 'dial'")
  if (!aon_method %in% c("bi", "nba", "d", "cbi", "cphast")) stop("aon_method should be 'bi', 'nba', 'd', 'cbi' or 'cphast'")

  if (aon_method == "nba" & is.null(Graph$coords)){
    aon_method <- "bi"
    message("nodes coordinates are not provided, using bidirectional Dijkstra")
  }
  phast <- FALSE

  if (aon_method == "d"){
    aon <- ifelse(length(unique(from_id)) <= length(unique(to_id)), 0, 1)
    contract <- FALSE
  }
  if (aon_method == "bi") {
    aon <- 2
    contract <- FALSE
  }
  if (aon_method == "nba") {
    aon <- 3
    contract <- FALSE
  }
  if (aon_method == "cbi") {
    aon <- 2
    contract <- TRUE
  }
  if (aon_method == "cphast"){
    aon <- ifelse(length(unique(from_id)) <= length(unique(to_id)), 0, 1)
    contract <- TRUE
    phast <- TRUE
  }

  vec <- rep(0, nrow(Graph$data))

  if (aon_method == "nba") {
    lat <- Graph$coords$X
    lon <- Graph$coords$Y

    if (constant == 1) warning("Are you sure constant is equal to 1 ?")
  } else{
    lat <- lon <- rep(0, Graph$nbnode)
  }


  if (algorithm %in% c("msa", "fw", "cfw", "bfw")){
    algo <- which(c("msa", "fw", "cfw", "bfw") == algorithm) - 1

    res <- cpptraffic(Graph$data$from, Graph$data$to, Graph$data$dist, vec, vec, Graph$data$dist,
                      Graph$attrib$alpha, Graph$attrib$beta, Graph$attrib$cap, Graph$nbnode,
                      lat, lon, constant, from_id, to_id, demand, max_gap, max_it, algo, aon, contract, phast, verbose)

  }

  if (algorithm == "dial"){
    inneriter <- 20
    max_tol <- 1e-11
    tmp <- paste0(tempdir(), "/cpp/")
    max_mem <- 8

    if (length(dial_params) > 0){
      params <- names(dial_params)
      if ("inneriter" %in% params){
        if (dial_params[["inneriter"]] >= 0) inneriter = dial_params[["inneriter"]]
      }
      if ("max_tol" %in% params){
        if (dial_params[["max_tol"]] >= 0) max_tol = dial_params[["max_tol"]]
      }
      if ("tmp_path" %in% params){
        if (!dir.exists(dial_params[["tmp_path"]])){
          dir.create(dial_params[["tmp_path"]], recursive = TRUE)
        }
        tmp <- dial_params[["tmp_path"]]

        if (!substr(tmp, nchar(tmp), nchar(tmp)) %in% c("/", "\\")) {
          tmp <- paste0(tmp, "/cpp/")
        } else{
          tmp <- paste0(tmp, "cpp/")
        }
      }
      dir.create(tmp, recursive = TRUE, showWarnings = FALSE)

      if ("max_mem" %in% params){
        if (dial_params[["max_mem"]] > 0) max_mem = dial_params[["max_mem"]]
      }
    }
    max_mem <- max_mem * 1e9
    one_bush <- Graph$nbnode * 4 + Graph$nbnode * 8 + nrow(Graph$data) * 4
    batch_size <- floor(max_mem/one_bush)
    n_origin <- ifelse(length(unique(from_id)) <= length(unique(to_id)), length(unique(from_id)), length(unique(to_id)))
    n_batch <- ceiling(n_origin/batch_size)
    if (n_batch == 1) batch_size <- n_origin

    reversed <- ifelse(length(unique(from_id)) <= length(unique(to_id)), FALSE, TRUE)

    if (reversed){
      gfrom <- Graph$data$to
      gto <- Graph$data$from
      dep <- to_id
      arr <- from_id
    } else{
      gfrom <- Graph$data$from
      gto <- Graph$data$to
      dep <- from_id
      arr <- to_id
    }


    res <- cppalgB(gfrom, gto, Graph$data$dist, vec, vec, Graph$data$dist,
                   Graph$attrib$alpha, Graph$attrib$beta, Graph$attrib$cap, Graph$nbnode,
                   lat, lon, constant, dep, arr, demand, max_gap, max_it, aon, batch_size, n_batch, tmp,
                   inneriter, max_tol, contract, phast, verbose)

    if (n_batch > 1){
      file.remove(list.files(tmp, full.names = TRUE))
    }

    if (reversed){
      res <- res[c(2,1,3,4,5,6,7,8,9,10)]
    }

  }

  final <- list()
  final[["gap"]] <- res[[9]]
  final[["iteration"]] <- res[[10]]
  res <- as.data.frame(res[c(1,2,3,4,5,6,7,8)], col.names = c("from", "to", "ftt", "cost", "flow", "capacity", "alpha", "beta"))
  res$from <- Graph$dict$ref[match(res$from, Graph$dict$id)]
  res$to <- Graph$dict$ref[match(res$to, Graph$dict$id)]

  final[["data"]] <- res

  return (final)

}