Ortools RoutingModel not finding best solution to a VRP in a 14-node example

ortools_vrp_4_depots_14_stops.cpp

#include "or_tools_route.h"
#include "ortools/constraint_solver/routing.h"

#include <vector>
#include <cstdio>

using operations_research::RoutingModel;

#define DEPOT_COUNT 4#include "or_tools_route.h"
#include "ortools/constraint_solver/routing.h"
#include "ortools/constraint_solver/routing_flags.h"

#include <vector>
#include <cstdio>

using operations_research::RoutingModel;

#define DEPOT_COUNT 4
#define STOP_COUNT 14
#define SWAP_ROW_0 0
#define SWAP_ROW_1 1

// The first four rows and columns of the distance matrix are the depots.
// In the test, we run with the DM as it is, and obtain an objective of
// 12409m, then swap the places of the depot 0 and depot 1 in the matrix
// and solve again, obtaining an objective of 11567.

// We think the optimal objective should not change with the ordering of the
// depots, and wonder if we are specifying the depots incorrectly.

// We've tested the same code with:
// - a variety of search parameters
// - a symmetric distance matrix (what's below is not)
// - different depot orderings in the depots vector when calling the
//   the RoutingModel constructor.
// None of this seems to make a difference.

long long distance_matrix[STOP_COUNT][STOP_COUNT] = {
  {    0, 1962, 3203, 1222,  714, 3798, 1879, 2729,  354, 2371, 3238, 2712, 2522, 1334 },
  { 1962,    0, 1811, 2350, 1541, 2406, 2436, 1640, 2080, 1282, 1921, 1065, 1082,  861 },
  { 3203, 1811,    0, 3591, 2816,  924, 3569, 2139, 3321, 2088, 1683, 1472, 1489, 2102 },
  { 1234, 2350, 3591,    0, 1640, 4133,  834, 3655, 1038, 3297, 3861, 3231, 3248, 2128 },
  {  714, 1541, 2816, 1640,    0, 3411, 2295, 2015, 1016, 1657, 2524, 1998, 1808,  800 },
  { 3798, 2406,  924, 4126, 3411,    0, 3982, 2409, 3916, 2576, 1843, 1893, 2023, 2697 },
  { 1889, 2436, 3575,  834, 2295, 3988,    0, 3892, 1693, 3534, 3931, 3301, 3318, 2441 },
  { 2729, 1633, 2139, 3655, 2015, 2409, 3885,    0, 3031,  360, 1114,  937,  974, 1944 },
  {  354, 2080, 3321, 1038, 1016, 3916, 1693, 3031,    0, 2673, 3540, 2961, 2824, 1504 },
  { 2371, 1275, 2088, 3297, 1657, 2576, 3527,  360, 2673,    0, 1404,  881,  721, 1586 },
  { 3238, 1921, 1683, 3861, 2524, 1843, 3931, 1114, 3540, 1404,    0, 1000, 1325, 2372 },
  { 2712, 1065, 1472, 3193, 1998, 1893, 3263,  937, 2923,  881, 1000,    0,  469, 1704 },
  { 2522, 1082, 1489, 3210, 1808, 2023, 3280,  974, 2824,  721, 1325,  469,    0, 1721 },
  { 1334,  861, 2102, 2128,  800, 2697, 2441, 1944, 1504, 1586, 2372, 1742, 1737,    0 }
};


class trip_cost_object
{
  public: long long trip_cost(RoutingModel::NodeIndex a, RoutingModel::NodeIndex b)
  {
    return distance_matrix[a.value()][b.value()];
  }
};


void solve(int run_number)
{
  trip_cost_object tco;
  std::vector<std::pair<RoutingModel::NodeIndex, RoutingModel::NodeIndex> > depots;

  for (int depoti = 0; depoti < DEPOT_COUNT; ++depoti)
    depots.push_back(std::make_pair(RoutingModel::NodeIndex(depoti), RoutingModel::NodeIndex(depoti)));

  RoutingModel R(STOP_COUNT, DEPOT_COUNT, depots);
  R.SetArcCostEvaluatorOfAllVehicles(NewPermanentCallback(&tco, &trip_cost_object::trip_cost));

  auto P =
      operations_research::BuildSearchParametersFromFlags();
  P.set_first_solution_strategy(
      operations_research::FirstSolutionStrategy::PARALLEL_CHEAPEST_INSERTION);

  R.CloseModel();

  auto plan = R.SolveWithParameters(P);
  printf("Run %d, Objective: %lld m\n", run_number, plan->ObjectiveValue());

  std::vector<std::vector<RoutingModel::NodeIndex>> routes;
  R.AssignmentToRoutes(*plan, &routes);

  int obj = 0;
  for (int i = 0; i < routes.size(); ++i)
  {
    int prev = i;
    printf("  Depot %d: ", i);
    for (int j = 0; j < routes[i].size(); ++j)
    {
      int next = routes[i][j].value();
      printf("(%4lld m) %2d ", distance_matrix[prev][next], next);
      obj += distance_matrix[prev][next];
      prev = next;
    }
    obj += distance_matrix[prev][i];
    printf("(%4lld m)\n", distance_matrix[prev][i]);
  }
  printf("  Measured objective: %d m\n\n", obj);
}


void swap_depots()
{
  long long temp;
  for (int i = 0; i < STOP_COUNT; ++i)
  {
    temp = distance_matrix[i][SWAP_ROW_0];
    distance_matrix[i][SWAP_ROW_0] = distance_matrix[i][SWAP_ROW_1];
    distance_matrix[i][SWAP_ROW_1] = temp;
  }
  for (int i = 0; i < STOP_COUNT; ++i)
  {
    temp = distance_matrix[SWAP_ROW_0][i];
    distance_matrix[SWAP_ROW_0][i] = distance_matrix[SWAP_ROW_1][i];
    distance_matrix[SWAP_ROW_1][i] = temp;
  }
}


int main()
{
  solve(1);
  swap_depots();
  solve(2);
  swap_depots();
  solve(3);

  return 0;
}

#define STOP_COUNT 14
#define SWAP_ROW_0 0
#define SWAP_ROW_1 1

// The first four rows and columns of the distance matrix are the depots.
// In the test, we run with the DM as it is, and obtain an objective of
// 12409m, then swap the places of the depot 0 and depot 1 in the matrix
// and solve again, obtaining an objective of 11567.

// We think the optimal objective should not change with the ordering of the
// depots, and wonder if we are specifying the depots incorrectly.

// We've tested the same code with:
// - a variety of search parameters
// - a symmetric distance matrix (what's below is not)
// - different depot orderings in the depots vector when calling the
//   the RoutingModel constructor.
// None of this seems to make a difference.

long long distance_matrix[STOP_COUNT][STOP_COUNT] = {
  {    0, 1962, 3203, 1222,  714, 3798, 1879, 2729,  354, 2371, 3238, 2712, 2522, 1334 },
  { 1962,    0, 1811, 2350, 1541, 2406, 2436, 1640, 2080, 1282, 1921, 1065, 1082,  861 },
  { 3203, 1811,    0, 3591, 2816,  924, 3569, 2139, 3321, 2088, 1683, 1472, 1489, 2102 },
  { 1234, 2350, 3591,    0, 1640, 4133,  834, 3655, 1038, 3297, 3861, 3231, 3248, 2128 },
  {  714, 1541, 2816, 1640,    0, 3411, 2295, 2015, 1016, 1657, 2524, 1998, 1808,  800 },
  { 3798, 2406,  924, 4126, 3411,    0, 3982, 2409, 3916, 2576, 1843, 1893, 2023, 2697 },
  { 1889, 2436, 3575,  834, 2295, 3988,    0, 3892, 1693, 3534, 3931, 3301, 3318, 2441 },
  { 2729, 1633, 2139, 3655, 2015, 2409, 3885,    0, 3031,  360, 1114,  937,  974, 1944 },
  {  354, 2080, 3321, 1038, 1016, 3916, 1693, 3031,    0, 2673, 3540, 2961, 2824, 1504 },
  { 2371, 1275, 2088, 3297, 1657, 2576, 3527,  360, 2673,    0, 1404,  881,  721, 1586 },
  { 3238, 1921, 1683, 3861, 2524, 1843, 3931, 1114, 3540, 1404,    0, 1000, 1325, 2372 },
  { 2712, 1065, 1472, 3193, 1998, 1893, 3263,  937, 2923,  881, 1000,    0,  469, 1704 },
  { 2522, 1082, 1489, 3210, 1808, 2023, 3280,  974, 2824,  721, 1325,  469,    0, 1721 },
  { 1334,  861, 2102, 2128,  800, 2697, 2441, 1944, 1504, 1586, 2372, 1742, 1737,    0 }
};


class trip_cost_object
{
  public: long long trip_cost(RoutingModel::NodeIndex a, RoutingModel::NodeIndex b)
  {
    return distance_matrix[a.value()][b.value()];
  }
};


void solve(int run_number)
{
  trip_cost_object tco;
  std::vector<std::pair<RoutingModel::NodeIndex, RoutingModel::NodeIndex> > depots;

  for (int depoti = 0; depoti < DEPOT_COUNT; ++depoti)
    depots.push_back(std::make_pair(RoutingModel::NodeIndex(depoti), RoutingModel::NodeIndex(depoti)));

  RoutingModel R(STOP_COUNT, DEPOT_COUNT, depots);
  R.SetArcCostEvaluatorOfAllVehicles(NewPermanentCallback(&tco, &trip_cost_object::trip_cost));

  auto plan = R.Solve();
  printf("Run %d, Objective: %lld m\n", run_number, plan->ObjectiveValue());

  std::vector<std::vector<RoutingModel::NodeIndex>> routes;
  R.AssignmentToRoutes(*plan, &routes);

  int obj = 0;
  for (int i = 0; i < routes.size(); ++i)
  {
    int prev = i;
    printf("  Depot %d: ", i);
    for (int j = 0; j < routes[i].size(); ++j)
    {
      int next = routes[i][j].value();
      printf("(%4lld m) %2d ", distance_matrix[prev][next], next);
      obj += distance_matrix[prev][next];
      prev = next;
    }
    obj += distance_matrix[prev][i];
    printf("(%4lld m)\n", distance_matrix[prev][i]);
  }
  printf("  Measured objective: %d m\n\n", obj);
}


void swap_depots()
{
  long long temp;
  for (int i = 0; i < STOP_COUNT; ++i)
  {
    temp = distance_matrix[i][SWAP_ROW_0];
    distance_matrix[i][SWAP_ROW_0] = distance_matrix[i][SWAP_ROW_1];
    distance_matrix[i][SWAP_ROW_1] = temp;
  }
  for (int i = 0; i < STOP_COUNT; ++i)
  {
    temp = distance_matrix[SWAP_ROW_0][i];
    distance_matrix[SWAP_ROW_0][i] = distance_matrix[SWAP_ROW_1][i];
    distance_matrix[SWAP_ROW_1][i] = temp;
  }
}


int main()
{
  solve(1);
  swap_depots();
  solve(2);
  swap_depots();
  solve(3);

  return 0;
}