Introduction

Flow Free is a game where given a grid and a number of (origin, destination) pairs you have to design flows that connect those OD paris and cover most of the edges in that grid. See Wikipedia for more. Flow Free is a game by Big Duck Games, LLC.

This is another test for the ompr package to discover bugs and improve the api.

Load some libraries

library(ompr)
library(ompr.roi)
library(ROI.plugin.glpk)
library(magrittr)
library(dplyr)
library(ggplot2)

Input data

start_positions <- tibble::tribble(
  ~start, ~end, ~color,
  1,   16, 1,
  2,   3, 2,
  4,   12, 3
)

This means, we have three flows with three different OD pairs.

Build the model

This model is just a prototype. It might be able to formulate everything a more clever.

# some helper parameters for the model
n <- 16
sn <- sqrt(n)
start_positions <- mutate(start_positions, x1 = (start - 1) %% sn + 1,
                           x2 = (end - 1) %% sn + 1,
                           y1 = ceiling(start / sn),
                           y2 = ceiling(end / sn))
s <- function(color) as.integer(start_positions[start_positions$color == color, "start"])
e <- function(color) as.integer(start_positions[start_positions$color == color, "end"])
m <- nrow(start_positions)

model <- MIPModel() %>%
  # 1 iff edge i,j is covered by flow color k
  add_variable(x[i, j, k], i = 1:n, j = 1:n, k = 1:m,
               type = "integer", lb = 0, ub = 0) %>%
  
  # a helper variable to model a flow to have connected flows
  add_variable(u[i], i = 1:n, lb = 0, ub = n) %>%
  
  # maximize flow
  set_objective(sum_exp(u[s(k)], k = 1:m)) %>%

  # no flow to itself
  add_constraint(x[i, i, k] == 0, i = 1:n, k = 1:m) %>%

  # only one direction allowed
  add_constraint(x[i, j, k] + x[j, i, k] <= 1, i = 1:n, j = 1:n, k = 1:m) %>%

  # only one flow per node
  add_constraint(sum_exp(x[i, j, k], k = 1:m) <= 1, i = 1:n, j = 1:n) %>%

  # at most one flow can start from a node
  add_constraint(sum_exp(x[i, j, k], k = 1:m, j = 1:n) <= 1, i = 1:n) %>%

  # at most one flow can go to a node
  add_constraint(sum_exp(x[i, j, k], k = 1:m, i = 1:n) <= 1, j = 1:n) %>%

  # if we travel an arc, we need to consume one flow
  add_constraint(u[i] - u[j] >= 1 - n * (1 - x[i, j, k]), i = 1:n, j = 1:n, k = 1:m) %>%
  add_constraint(u[i] - u[j] <= 1 + n * (1 - x[i, j, k]), i = 1:n, j = 1:n, k = 1:m) %>%
  add_constraint(u[e(k)] == 0, k = 1:m) %>% # every end point is a sink

  # set the start/end node
  add_constraint(sum_exp(x[s(k), j, k], j = 1:n) == 1, k = 1:m) %>%
  add_constraint(sum_exp(x[i, s(k), k], i = 1:n) == 0, k = 1:m) %>%
  add_constraint(sum_exp(x[i, e(k), k], i = 1:n) == 1, k = 1:m) %>%
  add_constraint(sum_exp(x[e(k), j, k], j = 1:n) == 0, k = 1:m)

# disable diagonals
for (i in 1:n) {
  for (j in 1:n) {
    x1 <- (j - 1) %% sn + 1
    y1 <- ceiling(j / sn)
    x2 <- (i - 1) %% sn + 1
    y2 <- ceiling(i / sn)
    d <- sqrt((x1 - x2)^2 + (y1 - y2)^2)
    if (d == 1) { # horizontal
      model <- set_bounds(model, x[i, j, k], ub = 1, i = i, j = j, k = 1:m)
    }
  }
}

# now add flow constraints to all non start/edges
for (k in 1:m) {
  start <- s(k)
  end <- e(k)
  i_wo_se <- setdiff(1:n, c(start, end))
  model <- add_constraint(model, sum_exp(x[j, i, k], j = 1:n) -
                            sum_exp(x[i, j, k], j = 1:n) == 0, i = i_wo_se)
}
model
## Mixed linear integer optimization problem
## Variables:
##   Continuous: 16 
##   Integer: 768 
##   Binary: 0 
## Search direction: maximize 
## Constraints: 2697

Solve the model

solution <- solve_model(model, with_ROI("glpk", verbose = TRUE)) %>%
  get_solution(x[i, j, k]) %>%
  filter(value > 0.8)
## <SOLVER MSG>  ----
## GLPK Simplex Optimizer, v4.60
## 2697 rows, 784 columns, 9711 non-zeros
##       0: obj =  -0.000000000e+00 inf =   6.000e+00 (6)
##      59: obj =   2.000000000e+00 inf =   0.000e+00 (0)
## *    86: obj =   3.400000000e+01 inf =   0.000e+00 (0)
## OPTIMAL LP SOLUTION FOUND
## GLPK Integer Optimizer, v4.60
## 2697 rows, 784 columns, 9711 non-zeros
## 768 integer variables, 144 of which are binary
## Integer optimization begins...
## +    86: mip =     not found yet <=              +inf        (1; 0)
## +   113: >>>>>   1.300000000e+01 <=   3.400000000e+01 161.5% (6; 0)
## +   275: mip =   1.300000000e+01 <=     tree is empty   0.0% (0; 39)
## INTEGER OPTIMAL SOLUTION FOUND
## <!SOLVER MSG> ----
solution
##    variable  i  j k value
## 1         x  1  5 1     1
## 2         x  5  6 1     1
## 3         x 10  9 1     1
## 4         x  6 10 1     1
## 5         x  9 13 1     1
## 6         x 13 14 1     1
## 7         x 14 15 1     1
## 8         x 15 16 1     1
## 9         x  2  3 2     1
## 10        x  8  7 3     1
## 11        x  4  8 3     1
## 12        x  7 11 3     1
## 13        x 11 12 3     1

Plot the solution

# some helper data for the plot
grid_data <- expand.grid(x = 1:sn, y = 1:sn) %>%
  left_join(start_positions, by = c("x" = "x1", "y" = "y1")) %>%
  left_join(start_positions, by = c("x" = "x2", "y" = "y2")) %>%
  mutate(color = if_else(is.na(color.x), color.y, color.x),
         is_start = !is.na(start.x),
         is_end = !is.na(end.y),
         node_type = if_else(is_start, "start", if_else(is_end, "end", "regular"))) %>%
  select(x, y, color, node_type)

solution <- solution %>%
  mutate(x1 = (i - 1) %% sn + 1,
         x2 = (j - 1) %% sn + 1,
         y1 = ceiling(i / sn),
         y2 = ceiling(j / sn))

non_flow_grid <- filter(grid_data, node_type == "regular") %>% 
  anti_join(solution, by = c("x" = "x1", "y" = "y1")) %>% 
  anti_join(solution, by = c("x" = "x2", "y" = "y2"))

flow_grid <- grid_data %>% 
  anti_join(non_flow_grid, by = c("x", "y"))

Start parameters

ggplot(NULL, aes(x = x, y = y)) +
  geom_point(data = grid_data, color = "gray", size = 4) +
  geom_point(data = start_positions, aes(x = x1, y = y1, color = factor(color)), size = 4) +
  geom_point(data = start_positions, aes(x = x2, y = y2, color = factor(color)), size = 4) +
  theme(text = element_text(color = "#ffffff"),
        rect = element_rect(fill = "#000000", color = "#000000"),
        plot.background = element_rect(fill = "#000000", color = "#000000"),
        panel.background = element_rect(fill = "#000000", color = "#000000"),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none") +
  ggtitle("Flow Free with R and GLPK")

Solution

ggplot(NULL, aes(x = x, y = y)) +
  geom_point(data = non_flow_grid, color = "gray", size = 4) +
  geom_segment(data = solution, aes(x = x1, y = y1,
                                    xend = x2, yend = y2, color = factor(k)), size = 3) +
  geom_point(data = solution, aes(x = x1, y = y1, color = factor(k)), size = 4) +
  geom_point(data = solution, aes(x = x2, y = y2, color = factor(k)), size = 4) +
  theme(text = element_text(color = "#ffffff"),
        rect = element_rect(fill = "#000000", color = "#000000"),
        plot.background = element_rect(fill = "#000000", color = "#000000"),
        panel.background = element_rect(fill = "#000000", color = "#000000"),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none") +
  ggtitle("Flow Free with R and GLPK")

Lessons learned for OMPR