# In dirkschumacher/romp: Model and Solve Mixed Integer Linear Programs

## Introduction

In this article we will model the minimum graph coloring problem. The goal: color a map with as few colors as possible while no two adjacent regions having the same color.

First, let's load some useful packages needed for (spatial) data processing.

```library(rgeos)
library(rgdal)
library(maptools)
library(dplyr)
```

Then we read in the 50 states of the US.

```# data from https://github.com/datasets/geo-boundaries-us-110m
# reference data from here: https://github.com/nvkelso/natural-earth-vector/tree/master/110m_cultural
```

Next step is to create an adjacency list to determine neighboring states.

```# this gives as an adjancy list
neighbors <- spdep::poly2nb(map_data)

# a helper function that determines if two nodes are adjacent
purrr::map2_lgl(i, j, ~ .y %in% neighbors[[.x]])
}
```

```is_adjacent(1, 2)
```

## Optimization model

Next, we will model the problem with `ompr` as a mixed integer linear program that tries to find a coloring with as few colors as possible.

```library(ompr)
```
```n <- nrow(map_data@data) # number of nodes aka states
max_colors <- 4 # 4 should be enough. But you increase this number.

# based on the formulation from here
# http://wwwhome.math.utwente.nl/~uetzm/do/IP-FKS.pdf
model <- MIPModel() %>%

# 1 iff node i has color k
add_variable(x[i, k], type = "binary", i = 1:n, k = 1:max_colors) %>%

# 1 iff color k is used
add_variable(y[k], type = "binary", k = 1:max_colors) %>%

# minimize colors
# multiply by k for symmetrie breaking (signifcant diff. in solution time)
set_objective(sum_expr(k * y[k], k = 1:max_colors), sense = "min") %>%

# each node is colored
add_constraint(sum_expr(x[i, k], k = 1:max_colors) == 1, i = 1:n) %>%

# if a color k is used, set y[k] to 1
add_constraint(x[i, k] <= y[k], i = 1:n, k = 1:max_colors) %>%

# no adjacent nodes have the same color
add_constraint(x[i, k] + x[j, k] <= 1, i = 1:n, j = 1:n, k = 1:max_colors, is_adjacent(i, j))
```
```model
```

The last constraint uses filter expressions:

```add_constraint(x[i, k] + x[j, k] <= 1, i = 1:n, j = 1:n, k = 1:max_colors, is_adjacent(i, j))
```

This code generates a constraint for each combination of i, j and k with the additional constraint, that `is_adjacent(i, j)` needs to evaluate to `TRUE`. Using filter expressions let's you design a very specific set of constraints with minimal code. This principle can be used with all functions that acceppt quantifiers: `sum_expr`, `set_bounds`, `add_constraint` and `add_variable`.

## Solve it

Now since we have built the model, we will solve it using the excellent `ROI` package and the `GLPK` solver. `ompr.roi` is a solver that takes an `ompr` package and translates it to `ROI`: this makes it possible to use range of solvers out of the box.

```library(ROI)
library(ROI.plugin.glpk)
library(ompr.roi)
```
```result <- solve_model(model, with_ROI("glpk", presolve = TRUE, verbose = TRUE))
```
```result
```

Yay. We found the minimal coloring with `r objective_value(result)` colors.

## Plot the result

Last step is to plot the result. First we will get the colors from the optimal solution.

```assigned_colors <- get_solution(result, x[i, k]) %>%
filter(value > 0.9) %>%
arrange(i)
```
```head(assigned_colors, 5)
```

Then we need to prepare the data for ggplot and join the colors to the data.

```library(ggplot2)
color_data <- map_data@data
color_data\$color <- assigned_colors\$k
plot_data_fort <- fortify(map_data, region = "adm1_code") %>%
by = c("id" = "adm1_code")) %>%
mutate(color = factor(color))
```

Now we have everything to plot it:

```ggplot(plot_data_fort, aes(x = long, y = lat, group = group)) +
geom_polygon(aes(fill = color)) +
coord_quickmap() +
viridis::scale_fill_viridis(discrete = TRUE, option = "D")
```

## Feedback

Do you have any questions, ideas, comments? Or did you find a mistake? Let's discuss on Github.

dirkschumacher/romp documentation built on Oct. 6, 2018, 4:36 a.m.