Loading the data

The raw results can be found in the file results.csv with the following columns:

• tag: a tag which should always equal “result” used to grep the result line in the execution log file.
• instance: the name or path of the instance.
• n_facilities: the number of facilities in the instance.
• n_customers: the number of customers in the instance.
• Gamma: the value of \(\Gamma\) which was used.
• deviation: the maximum deviation, in percentage, from the nominal demand.
• time_limit: the time limit which was used when solving the instance.
• method: the name of the method which was used to solve the instance.
• total_time: the total time used to solve the instance.
• master_time: the time spent solving the master problem.
• separation_time: the time spent solving the separation problem.
• best_bound: the best bound found.
• iteration_count: the number of iterations.
• fail: should be empty if the execution of the algorithm went well.

We start by reading the file and by removing the “tag” column.

results = read.table("./results.csv", header = FALSE, sep = ',')
colnames(results) <- c("tag", "instance", "n_facilities", "n_customers", "Gamma", "deviation", "time_limit", "method", "use_heuristic", "total_time", "master_time", "separation_time", "best_bound", "iteration_count", "fail")
results$tag = NULL

We then check that all instances were solved without issue by checking the “fail” column.

sum(results$fail)

## [1] 0

Then, we compute the percentage which \(\Gamma\) represents for the total number of customers

results$p = results$Gamma / results$n_customers
results$p = ceiling(results$p / .05) * .05

results = results %>%
  mutate(
    ratio_demand_capacity = as.integer(sub('.*instance_F\\d+_C\\d+_R(\\d+)__\\d+\\.txt', '\\1', instance))
  )

results$size = paste0("(", results$n_facilities, ",", results$n_customers, "), ", results$p)

results$method_extended = paste0(results$method, ", ", results$use_heuristic)

We add a tag for unsolved instances.

time_limit = 7200

results$unsolved = results$total_time >= time_limit | results$fail

All in all, our result data reads.

Unsolved instances

for (method in unique(results$method_extended)) {
  
  # Sum of unsolved cases for each size
  sum_unsolved = results[results$method_extended == method,] %>% group_by(size) %>% summarise(total_unsolved = sum(unsolved))
  
  # Create a bar plot
  p = ggplot(sum_unsolved, aes(x = size, y = total_unsolved)) +
    geom_bar(stat = "identity") +
    labs(x = "Size", y = "Total Unsolved Cases", title = method) +
    theme_minimal()

  ggsave(paste0("unsolved_", method, ".pdf"), plot = p, width = 10, height = 6)
  
  print(p)
}

Performance profiles and ECDF

We now introduce a function which plots the performance profile of our solvers over a given test set.

add_performance_ratio = function(dataset, 
                                 criterion_column = "total_time",
                                 unsolved_column = "unsolved",
                                 instance_column = "instance",
                                 solver_column = "solver",
                                 output_column = "performance_ratio") {
  
  # Compute best score for each instance
  best = dataset %>%
    group_by(!!sym(instance_column)) %>%
    mutate(best_solver = min(!!sym(criterion_column)))
  
  # Compute performance ratio for each instance and solver
  result = best %>%
    group_by(!!sym(instance_column), !!sym(solver_column)) %>%
    mutate(!!sym(output_column) := !!sym(criterion_column) / best_solver) %>%
    ungroup()
  
  if (sum(result[,unsolved_column]) > 0) {
    result[result[,unsolved_column] == TRUE,output_column] = max(result[,output_column])
  }

  return (result)
}

plot_performance_profile = function(dataset,
                                    criterion_column,
                                    unsolved_column = "unsolved",
                                    instance_column = "instance",
                                    solver_column = "solver"
                                    ) {
  
  dataset_with_performance_ratios = add_performance_ratio(dataset,
                                                          criterion_column = criterion_column,
                                                          instance_column = instance_column,
                                                          solver_column = solver_column,
                                                          unsolved_column = unsolved_column)
  
  solved_dataset_with_performance_ratios = dataset_with_performance_ratios[!dataset_with_performance_ratios[,unsolved_column],]
  
  compute_performance_profile_point = function(method, data) {
    
    performance_ratios = solved_dataset_with_performance_ratios[solved_dataset_with_performance_ratios[,solver_column] == method,]$performance_ratio
    
    unscaled_performance_profile_point = ecdf(performance_ratios)(data)
    
    n_instances = sum(dataset[,solver_column] == method)
    n_solved_instances = sum(dataset[,solver_column] == method & !dataset[,unsolved_column])
    
    return( unscaled_performance_profile_point * n_solved_instances / n_instances )
  }
  
  perf = solved_dataset_with_performance_ratios %>%
    group_by(!!sym(solver_column)) %>%
    mutate(performance_profile_point = compute_performance_profile_point(unique(!!sym(solver_column)), performance_ratio))
  
  result = ggplot(data = perf, aes(x = performance_ratio, y = performance_profile_point, color = !!sym(solver_column))) +
              geom_line()
  
  return (result)
}

performance_profile = plot_performance_profile(results, criterion_column = "total_time", solver_column = "method_extended") +
  labs(x = "Performance ratio", y = "% of instances") +
  scale_y_continuous(limits = c(0, 1)) +
  theme_minimal()
print(performance_profile)

ggsave("performance_profile.pdf", plot = performance_profile)

## Saving 7 x 5 in image

for (ratio in unique(results$ratio_demand_capacity)) {
  
  p = plot_performance_profile(results[results$ratio_demand_capacity == ratio,], criterion_column = "total_time", solver_column = "method_extended") +
    labs(x = "Performance ratio", y = "% of instances", title = paste0("R = ", ratio / 1000)) +
    scale_y_continuous(limits = c(0, 1)) +
    theme_minimal()
  
  print(p)
}

As a complement, we also draw the ECDF.

ggplot(results, aes(x = total_time, color = method_extended)) +
  stat_ecdf(geom = "step") +
  xlab("Total Time") +
  ylab("ECDF") +
  labs(title = "ECDF of Total Time by Method")  +
  scale_x_continuous(breaks = seq(0, max(results$total_time), by = 500), labels = seq(0, max(results$total_time), by = 500)) +
  theme_minimal()

results_solved = results %>%
  filter(total_time < 7200) %>%
  group_by(method_extended, n_facilities, n_customers, p, deviation) %>%
  summarize(
    solved = n(),
    total_time = mean(total_time),
    master_time = mean(master_time),
    separation_time = mean(separation_time),
    solved_iteration_count = mean(iteration_count),
    .groups = "drop"
  ) %>%
  ungroup()

results_unsolved <- results %>%
  filter(total_time >= 7200) %>%
  group_by(method_extended, n_facilities, n_customers, p, deviation) %>%
  summarize(
    unsolved = n(),
    unsolved_iteration_count = mean(iteration_count),
    .groups = "drop"
  ) %>%
  ungroup()

final_result <- merge(results_solved, results_unsolved, by = c("method_extended", "n_facilities", "n_customers", "p", "deviation"), all = TRUE)

final_result[is.na(final_result)] = 0

Number of iterations

for (method in unique(results$method_extended)) {
  
  p = ggplot(results[results$method_extended == method & results$total_time < time_limit,], aes(x = iteration_count)) +
    geom_bar() +
    labs(
      title = method,
      x = "Number of iterations",
      y = "Number of instances"
    ) +
    theme_minimal()
    
  print(p) 
}

Understanding unsolved instances

We first consider those instances with 15 facilities and p = 0.3 since these are the instances which both approaches cannot solve at all.

subset = results[results$n_facilities == 15 & results$p > .25,]

Then, observe how, for the subset of instances such that iteration_count is 1 for at least one of the methods, the computation times are the same (which makes sense because, actually, the same routines are executed - intial master and first separation).

filtered_subset = subset %>%
  group_by(instance) %>%
  filter(any(iteration_count < 1))

ggplot(filtered_subset, aes(x = method, y = master_time)) +
  geom_boxplot() +
  labs(x = "Method", y = "Master time") +
  theme_minimal()

ggplot(filtered_subset, aes(x = method, y = separation_time)) +
  geom_boxplot() +
  labs(x = "Method", y = "Separation time") +
  theme_minimal()

However, for some instances, GBD can do more iterations because the master problem is nicer whereas CCG has a convex MINLP problem which he cannot solve within the time limit.

for (method in unique(subset$method_extended)) {
  
  p = ggplot(subset[subset$method_extended == method,], aes(x = iteration_count)) +
    geom_bar() +
    labs(
      title = method,
      x = "Number of iterations",
      y = "Number of instances"
    ) +
    theme_minimal()
    
  print(p) 
}