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

Comparison with Static Robust Problem

static = read.table("./static-rap-results.csv", header = FALSE, sep = ',')
colnames(static) <- 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")
static$tag = NULL 

columns_to_keep = c("instance", "Gamma", "deviation", "n_servers", "n_clients", "method", "total_time", "best_bound", "iteration_count")
ccg = all_results %>% filter(method == "CCG", use_heuristic == 0, use_bilevel_separation == 0, use_budgeted_unc == 1) %>% select(all_of(columns_to_keep), p)
static = static %>% select(columns_to_keep)

## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(columns_to_keep)
## 
##   # Now:
##   data %>% select(all_of(columns_to_keep))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ccg = ccg %>% mutate(instance = basename(instance))
static = static %>% mutate(instance = basename(instance))

merged_df = merge(ccg, static,
                   by = c("instance", "Gamma", "deviation"), 
                   all = TRUE,
                   suffixes = c(".ccg", ".static"))
merged_df %>% filter(best_bound.ccg > best_bound.static)

##                   instance Gamma deviation n_servers.ccg n_clients.ccg
## 1  instance_S10_C20__2.txt     6      0.50            10            20
## 2  instance_S10_C20__4.txt     6      0.50            10            20
## 3  instance_S10_C30__0.txt     3      0.25            10            30
## 4  instance_S10_C30__0.txt     6      0.25            10            30
## 5  instance_S10_C30__0.txt     6      0.50            10            30
## 6  instance_S10_C30__0.txt     9      0.25            10            30
## 7  instance_S10_C30__0.txt     9      0.50            10            30
## 8  instance_S10_C30__7.txt     3      0.25            10            30
## 9  instance_S10_C30__9.txt     9      0.25            10            30
## 10 instance_S10_C40__0.txt     4      0.50            10            40
## 11 instance_S10_C40__1.txt     4      0.25            10            40
## 12 instance_S10_C40__1.txt     4      0.50            10            40
## 13 instance_S10_C40__2.txt     4      0.25            10            40
## 14 instance_S10_C40__2.txt     4      0.50            10            40
## 15 instance_S10_C40__2.txt     8      0.50            10            40
## 16 instance_S10_C40__3.txt     4      0.50            10            40
## 17 instance_S10_C40__4.txt     4      0.25            10            40
## 18 instance_S10_C40__4.txt     4      0.50            10            40
## 19 instance_S10_C40__6.txt     4      0.25            10            40
## 20 instance_S10_C40__6.txt     4      0.50            10            40
## 21 instance_S10_C40__7.txt     4      0.25            10            40
## 22 instance_S10_C40__7.txt     4      0.50            10            40
## 23 instance_S10_C40__8.txt     4      0.50            10            40
## 24 instance_S10_C40__9.txt     4      0.25            10            40
## 25 instance_S15_C40__8.txt     4      0.50            15            40
##    method.ccg total_time.ccg best_bound.ccg iteration_count.ccg   p
## 1         CCG           7200          1e+20                   2 0.3
## 2         CCG           7200          1e+20                  54 0.3
## 3         CCG           7200          1e+20                   2 0.1
## 4         CCG           7200          1e+20                   1 0.2
## 5         CCG           7200          1e+20                   2 0.2
## 6         CCG           7200          1e+20                   3 0.3
## 7         CCG           7200          1e+20                   2 0.3
## 8         CCG           7200          1e+20                   3 0.1
## 9         CCG           7200          1e+20                   2 0.3
## 10        CCG           7200          1e+20                   1 0.1
## 11        CCG           7200          1e+20                   1 0.1
## 12        CCG           7200          1e+20                   1 0.1
## 13        CCG           7200          1e+20                   1 0.1
## 14        CCG           7200          1e+20                   2 0.1
## 15        CCG           7200          1e+20                   1 0.2
## 16        CCG           7200          1e+20                   2 0.1
## 17        CCG           7200          1e+20                   1 0.1
## 18        CCG           7200          1e+20                   1 0.1
## 19        CCG           7200          1e+20                   2 0.1
## 20        CCG           7200          1e+20                   2 0.1
## 21        CCG           7200          1e+20                   3 0.1
## 22        CCG           7200          1e+20                   2 0.1
## 23        CCG           7200          1e+20                   3 0.1
## 24        CCG           7200          1e+20                   2 0.1
## 25        CCG           7200          1e+20                   1 0.1
##    n_servers.static n_clients.static method.static total_time.static
## 1                10               20        Static          0.073676
## 2                10               20        Static          0.017307
## 3                10               30        Static          0.017379
## 4                10               30        Static          0.015327
## 5                10               30        Static          0.012198
## 6                10               30        Static          0.017425
## 7                10               30        Static          0.013167
## 8                10               30        Static          0.078152
## 9                10               30        Static          0.075932
## 10               10               40        Static          0.074957
## 11               10               40        Static          0.016686
## 12               10               40        Static          0.017202
## 13               10               40        Static          0.080372
## 14               10               40        Static          0.081235
## 15               10               40        Static          0.081203
## 16               10               40        Static          0.078905
## 17               10               40        Static          0.084900
## 18               10               40        Static          0.086383
## 19               10               40        Static          0.336255
## 20               10               40        Static          0.239991
## 21               10               40        Static          0.068940
## 22               10               40        Static          0.068841
## 23               10               40        Static          0.076108
## 24               10               40        Static          0.063384
## 25               15               40        Static          0.081368
##    best_bound.static iteration_count.static
## 1            1129810                      0
## 2            1000650                      0
## 3            2154810                      0
## 4            2154810                      0
## 5            3088660                      0
## 6            2154810                      0
## 7            3088660                      0
## 8            1273510                      0
## 9            1974460                      0
## 10           3634090                      0
## 11           6488470                      0
## 12           9326280                      0
## 13           3835130                      0
## 14           5507930                      0
## 15           5507930                      0
## 16           4386170                      0
## 17           5312360                      0
## 18           7632320                      0
## 19           3464220                      0
## 20           4973340                      0
## 21           2407820                      0
## 22           3450730                      0
## 23           3408460                      0
## 24           2441350                      0
## 25           3531820                      0

merged_df %>% filter(best_bound.static < 1e-4)

##  [1] instance               Gamma                  deviation             
##  [4] n_servers.ccg          n_clients.ccg          method.ccg            
##  [7] total_time.ccg         best_bound.ccg         iteration_count.ccg   
## [10] p                      n_servers.static       n_clients.static      
## [13] method.static          total_time.static      best_bound.static     
## [16] iteration_count.static
## <0 rows> (or 0-length row.names)

library(tidyr)

# 1. Compute summary
summary_table <- merged_df %>%
  mutate(
    gap = (best_bound.static - best_bound.ccg) / (1e-10 + abs(best_bound.ccg)) * 100,
    time_status = ifelse(total_time.ccg < time_limit, "within_time", "timed_out")
  ) %>%
  group_by(n_servers.ccg, n_clients.ccg, p, deviation, time_status) %>%
  summarise(
    n_instances = length(best_bound.static),
    n_missing_static = sum(is.na(best_bound.static)),
    avg_ccg_time = mean(total_time.ccg),
    avg_static_time = mean(total_time.static),
    avg_gap = ifelse(sum(!is.na(best_bound.static) & best_bound.ccg < 1e20) > 0, 
                 mean(gap[!is.na(best_bound.static) & best_bound.ccg < 1e20]), 
                 NA),
    .groups = "drop"
  )

# 2. Pivot wider
summary_wide <- summary_table %>%
  pivot_wider(
    names_from = time_status,
    values_from = c(n_instances, n_missing_static, avg_ccg_time, avg_static_time, avg_gap),
    names_glue = "{time_status}_{.value}"
  ) %>%
  arrange(n_servers.ccg, n_clients.ccg, p, deviation)

summary_wide = summary_wide %>% select(n_servers.ccg, n_clients.ccg, p, deviation, within_time_n_instances, within_time_n_missing_static, within_time_avg_ccg_time, within_time_avg_static_time, within_time_avg_gap, timed_out_n_instances, timed_out_n_missing_static, timed_out_avg_gap)
# 3. Pretty table with exact column names
kable(summary_wide, "html", digits = 4, escape = FALSE, col.names = c(
  "|V_1|", "|V_2|", "p", "dev",
  "n_instances", "n_static_is_infeasible", "avg_ccg_time", "avg_static_time", "avg_gap (%)",
  "n_instances", "n_static_is_infeasible", "avg_gap (%)"
)) %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  add_header_above(c(" " = 4, "Solved by CCG" = 5, "Time out" = 3))

larger = read.table("./larger-results.csv", header = FALSE, sep = ',')
colnames(larger) <- 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")
larger$tag = NULL 
larger$p = larger$Gamma / larger$n_clients
larger$p = ceiling(larger$p / .05) * .05

ccg = larger %>% filter(method == "CCG", use_heuristic == 1)
static = larger %>% filter(method == "Static")
static$p = NULL

merged_df = merge(ccg, static,
                   by = c("instance", "Gamma", "deviation"), 
                   all = TRUE,
                   suffixes = c(".ccg", ".static"))
merged_df %>% filter(best_bound.ccg > best_bound.static)

##  [1] instance                      Gamma                        
##  [3] deviation                     n_servers.ccg                
##  [5] n_clients.ccg                 time_limit.ccg               
##  [7] method.ccg                    use_heuristic.ccg            
##  [9] total_time.ccg                master_time.ccg              
## [11] separation_time.ccg           best_bound.ccg               
## [13] iteration_count.ccg           fail.ccg                     
## [15] use_bilevel_separation.ccg    use_budgeted_unc.ccg         
## [17] p                             n_servers.static             
## [19] n_clients.static              time_limit.static            
## [21] method.static                 use_heuristic.static         
## [23] total_time.static             master_time.static           
## [25] separation_time.static        best_bound.static            
## [27] iteration_count.static        fail.static                  
## [29] use_bilevel_separation.static use_budgeted_unc.static      
## <0 rows> (or 0-length row.names)

merged_df %>% filter(best_bound.static < 1e-4)

##  [1] instance                      Gamma                        
##  [3] deviation                     n_servers.ccg                
##  [5] n_clients.ccg                 time_limit.ccg               
##  [7] method.ccg                    use_heuristic.ccg            
##  [9] total_time.ccg                master_time.ccg              
## [11] separation_time.ccg           best_bound.ccg               
## [13] iteration_count.ccg           fail.ccg                     
## [15] use_bilevel_separation.ccg    use_budgeted_unc.ccg         
## [17] p                             n_servers.static             
## [19] n_clients.static              time_limit.static            
## [21] method.static                 use_heuristic.static         
## [23] total_time.static             master_time.static           
## [25] separation_time.static        best_bound.static            
## [27] iteration_count.static        fail.static                  
## [29] use_bilevel_separation.static use_budgeted_unc.static      
## <0 rows> (or 0-length row.names)

# 1. Compute summary
summary_table <- merged_df %>%
  mutate(
    gap = (best_bound.static - best_bound.ccg) / (1e-10 + abs(best_bound.ccg)) * 100,
    time_status = "timed_out",#ifelse(total_time.ccg < time_limit, "within_time", "timed_out")
  ) %>%
  group_by(n_servers.ccg, n_clients.ccg, p, deviation, time_status) %>%
  summarise(
    n_instances = length(best_bound.static),
    n_missing_static = sum(is.na(best_bound.static)),
    avg_ccg_time = mean(total_time.ccg),
    avg_static_time = mean(total_time.static),
    avg_gap = ifelse(sum(!is.na(best_bound.static) & best_bound.ccg < 1e20) > 0, 
                 mean(gap[!is.na(best_bound.static) & best_bound.ccg < 1e20]), 
                 NA),
    .groups = "drop"
  )

# 2. Pivot wider
summary_wide <- summary_table %>%
  pivot_wider(
    names_from = time_status,
    values_from = c(n_instances, n_missing_static, avg_ccg_time, avg_static_time, avg_gap),
    names_glue = "{time_status}_{.value}"
  ) %>%
  arrange(n_servers.ccg, n_clients.ccg, p, deviation)

summary_wide = summary_wide %>% select(n_servers.ccg, n_clients.ccg, p, deviation, timed_out_n_instances, timed_out_n_missing_static, timed_out_avg_gap)
# 3. Pretty table with exact column names
kable(summary_wide, "html", digits = 4, escape = FALSE, col.names = c(
  "|V_1|", "|V_2|", "p", "dev",
  "n_instances", "n_static_is_infeasible", "avg_gap (%)"
)) %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  add_header_above(c(" " = 4, "Time out" = 3))