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_servers: the number of servers in the instance.
• n_clients: the number of clients 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("./new-results.csv", header = FALSE, sep = ',')
colnames(results) <- c("tag", "instance", "n_servers", "n_clients", "Gamma", "deviation", "time_limit", "method", "use_heuristic", "total_time", "master_time", "separation_time", "best_bound", "iteration_count", "fail")
results$tag = NULL

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

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

results$size = paste0("(", results$n_servers, ",", results$n_clients, "), ", results$p)

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

We also get rid of instances which were too large to be solved by both approaches.

results = results[!(results$n_clients == 40),]
results = results[!(results$n_servers == 25 & results$n_clients == 25),]

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

sum(results$fail)

## [1] 62

We add a tag for unsolved instances.

time_limit = 7200

if (sum(results$fail) > 0) {
  results[results$fail == TRUE,]$total_time = time_limit
}

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

All in all, our result data reads.

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

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()

Summary table

In this section, we create a table summarizing the outcome of our experiments.

We start by computing average computation times over the solved instances.

results_solved = results %>%
  filter(total_time < time_limit) %>%
  group_by(method_extended, n_servers, n_clients, 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()

Then, we also consider unsolved instances: we count the number of such instances and compute the average iteration count.

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

results_fail <- results %>%
  filter(fail == TRUE) %>%
  group_by(method_extended, n_servers, n_clients, p, deviation) %>%
  summarize(
    fail = n(),
    .groups = "drop"
  ) %>%
  ungroup()

We then merge the two tables.

final_result <- merge(results_solved, results_unsolved, by = c("method_extended", "n_servers", "n_clients", "p", "deviation"), all = TRUE)
final_result <- merge(final_result, results_fail, by = c("method_extended", "n_servers", "n_clients", "p", "deviation"), all = TRUE)

Here is our table.

Revision

results_rev = read.table("./results-rev.csv", header = FALSE, sep = ',')
colnames(results_rev) <- c("tag", "instance", "n_servers", "n_clients", "Gamma", "deviation", "time_limit", "method", "use_heuristic", "total_time", "master_time", "separation_time", "best_bound", "iteration_count", "fail", "use_bilevel_separation", "use_budgeted_unc")

results = read.table("./new-results.csv", header = FALSE, sep = ',')
colnames(results) <- c("tag", "instance", "n_servers", "n_clients", "Gamma", "deviation", "time_limit", "method", "use_heuristic", "total_time", "master_time", "separation_time", "best_bound", "iteration_count", "fail")
results$use_bilevel_separation = 0
results$use_budgeted_unc = 1

all_results = rbind(results, results_rev)

all_results = all_results %>% filter(deviation == .25 & method == "CCG")

all_results$p = all_results$Gamma / all_results$n_clients
all_results$p = ceiling(all_results$p / .05) * .05

all_results$size = paste0("(", all_results$n_servers, ",", all_results$n_clients, "), ", all_results$p)

all_results$method_extended = paste0(all_results$method, ", ", ifelse(all_results$use_heuristic, "heuristic, ", ""), ifelse(all_results$use_bilevel_separation, "KKT", "0-1 mapping"))

paged_table(results_rev %>% filter(fail == 1 & total_time < time_limit))

paged_table(results_rev %>% filter(fail == 0 & deviation == .25))

if (sum(all_results$fail) > 0) {
  all_results[all_results$fail == TRUE,]$total_time = time_limit
}

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

budgeted_unc_results = all_results %>% filter(use_budgeted_unc == 1 & method == "CCG")


ggplot(budgeted_unc_results, aes(x = total_time, color = method_extended)) +
  stat_ecdf(geom = "step", linewidth = 1) +
  labs(title = "ECDF Plot", x = "Values", y = "ECDF") +
  theme_minimal()

performance_profile = plot_performance_profile(budgeted_unc_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)

knapsack_unc_results = all_results %>% filter(method == "CCG" & use_budgeted_unc == 0)

ggplot(knapsack_unc_results, aes(x = total_time, color = method_extended)) +
  stat_ecdf(geom = "step", size = 1) +
  labs(title = "ECDF Plot", x = "Values", y = "ECDF") +
  facet_wrap(~p) +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Summary table

In this section, we create a table summarizing the outcome of our experiments.

We start by computing average computation times over the solved instances.

results_solved = knapsack_unc_results %>%
  filter(total_time < time_limit) %>%
  group_by(method_extended, n_servers, n_clients, 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()

Then, we also consider unsolved instances: we count the number of such instances and compute the average iteration count.

results_unsolved <- knapsack_unc_results %>%
  filter(total_time >= time_limit) %>%
  group_by(method_extended, n_servers, n_clients, p, deviation) %>%
  summarize(
    unsolved = n(),
    unsolved_iteration_count = mean(iteration_count),
    .groups = "drop"
  ) %>%
  ungroup()

results_fail <- knapsack_unc_results %>%
  filter(fail == TRUE) %>%
  group_by(method_extended, n_servers, n_clients, p, deviation) %>%
  summarize(
    fail = n(),
    .groups = "drop"
  ) %>%
  ungroup()

We then merge the two tables.

final_result <- merge(results_solved, results_unsolved, by = c("method_extended", "n_servers", "n_clients", "p", "deviation"), all = TRUE)
final_result <- merge(final_result, results_fail, by = c("method_extended", "n_servers", "n_clients", "p", "deviation"), all = TRUE)

Here is our table.