使用库（ompr）将问题建模为使用Miller-Tucker-Zemlin（MTZ）公式的混合整数线性规划（MILP）。
使用库（ompr.roi）和库（ROI.plugin.glpk）来解决问题，GLPK（GNU Linear Programming Kit）解决了大规模MILP问题。
参考：https://dirkschumacher.github.io/ompr/articles/problem-tsp.html
这是一个使用传单的应用程序

library(leaflet)

city.1 <- c(city='Paris',  lat= 48.8566,  lng=2.3522)
city.2 <- c(city='Saumur',  lat= 47.2600, lng=-0.0769)
city.3 <- c(city='Rennes', lat= 48.1147, lng=-1.6794)
city.4 <- c(city='Caen', lat= 49.1800, lng=-0.3700)
city.5 <- c(city='Rouen', lat= 49.4428,  lng=1.0886)
D <- data.frame(rbind(city.1, city.2, city.3, city.4, city.5))
D$lat <- as.numeric(D$lat)
D$lng <- as.numeric(D$lng)

leaflet() %>%
  addTiles() %>% 
  addCircles(data=D, lng = ~lng, lat = ~lat, radius=5000, color='red',
             popup=~city, label=~city)
  
path <- TSP(D)
path

# on the map
  leaflet() %>%
    addTiles() %>% 
    addCircles(data=path, lng = ~lng, lat = ~lat, radius=5000, color='red',
               popup=~city, label=~city) %>%
    addPolylines(data = path, lng = ~lng, lat = ~lat, group = ~group, weight=2)

    
TSP <- function(D) {

  # distance matrix
    M <- matrix(0,nrow=nrow(D), ncol=nrow(D))
    rownames(M) <- D$city -> colnames(M)
    for (r in 1:nrow(D)) {
      pt <- c ( D$lng[r], D$lat[r] )
      E <- distanceTo(pt, D)
      M[r,] <- round(E$dist,0)
    }
#    M
    library(ompr)
    n <- nrow(M) # problem size 
    
    dist_fun <- function(i, j) {
      # a distance extraction function
      vapply(seq_along(i), function(k) M[i[k], j[k]], numeric(1L))
    }
    
    
    model <- MILPModel() %>%
      
      # variable == 1 iff we travel from city i to j
      add_variable(x[i, j], i = 1:n, j = 1:n, type = "integer", lb = 0, ub = 1) %>%
      
      # a helper variable for the MTZ formulation of the tsp
      add_variable(u[i], i = 1:n, lb = 1, ub = n) %>% 
      
      # minimize travel distance
      set_objective(sum_expr(colwise(dist_fun(i, j)) * x[i, j], i = 1:n, j = 1:n), "min") %>%
      
      # you cannot go to the same city
      set_bounds(x[i, i], ub = 0, i = 1:n) %>%
      
      # leave each city
      add_constraint(sum_expr(x[i, j], j = 1:n) == 1, i = 1:n) %>%
      
      # visit each city
      add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = 1:n) %>%
      
      # ensure no subtours (arc constraints)
      add_constraint(u[i] >= 2, i = 2:n) %>% 
      add_constraint(u[i] - u[j] + 1 <= (n - 1) * (1 - x[i, j]), i = 2:n, j = 2:n)
    
#    model    # builds the model
    
    
  # solve the model
    library(ompr.roi)
    library(ROI.plugin.glpk)
    result <- solve_model(model, with_ROI(solver = "glpk", verbose = TRUE))
    
  # extract the solution
    solution <- get_solution(result, x[i, j])
    path <- solution[which(solution$value > 0),]
    
    odd <- function(x) x%%2 != 0;    even <- function(x) x%%2 == 0
    
    from <- 1
    df <- data.frame()
    for (i in 1:nrow(path)) {
      to <- path$j[path$i == from]
#      print(paste("from", dest$city[from], " to", dest$city[to]))
      temp <- data.frame(city=D$city[from], lat=D$lat[from], lng=D$lng[from], group="A")
      if (even(i)) temp$group[1] <- "B" 
      df <- rbind(df, temp)
      from <- to
    }
    temp <- data.frame(city=D$city[from], lat=D$lat[from], lng=D$lng[from], group="A")
    if (odd(i)) temp$group[1] <- "B" 
    df <- rbind(df, temp)
    
    return(df)

}