Objectives

We want to allocate resources to zones while minimising the total cost (e.g. travel time, travel distance, expenditure, etc.) for the entire system.

Allocation Costs

I <- 10
J <- 5

my_costs <- matrix(runif(I*J), nrow = I, ncol = J)
rownames(my_costs) <- sprintf("R_%s", seq.int(from = 1,length.out = I))
colnames(my_costs) <- sprintf("Z_%s", seq.int(from = 1,length.out = J))

my_costs
##             Z_1        Z_2       Z_3       Z_4        Z_5
## R_1  0.31300760 0.13139723 0.4963993 0.7840266 0.73679765
## R_2  0.53842965 0.87764273 0.7447344 0.0672037 0.91149671
## R_3  0.25143509 0.90748595 0.1146217 0.5751526 0.53631711
## R_4  0.95993728 0.36790132 0.2225167 0.2831948 0.01454969
## R_5  0.30849284 0.74401005 0.9107487 0.3551488 0.52337242
## R_6  0.39425429 0.40168052 0.3757401 0.3075921 0.11172232
## R_7  0.04196323 0.07738462 0.7931153 0.5933098 0.76521537
## R_8  0.72600738 0.46235135 0.4789284 0.8689144 0.18859216
## R_9  0.68708525 0.88440675 0.7547065 0.2911243 0.36670425
## R_10 0.99598999 0.31697881 0.8720274 0.6002004 0.48280044

This system has \(I=10\) resources and \(J=5\) zones, so that the resources are denoted as \(\{R_1, R_2, R_3...R_I\}\) while the zones are denoted as \(\{Z_1, Z_2, Z_3...Z_J\}\).

A \(I \times J\) matrix denoted as \(c\) represents the allocation cost for each resource to each zone. The value \(c_{i,j}\) represents the cost of assigning resource \(R_i\) to zone \(Z_j\).

In this example, in order to assign resource \(R_2\) to zone \(Z_3\), a cost of \(c_{2,3}=0.745\) is incurred.

Minimum Capacity

Map showing layout of the zones

Map showing layout of the zones 

s <- c(2, 1, 3, 2, 1)

A integer vector \(s\) of length \(J\) will specify the minimum number of resources required for each zone. In this example, \(s_1=2\) means that zone \(Z_1\) requires at least 2 resources.

Problem Formulation

The Integer Programming (IP) problem is formulated as below.

The symbol \(x\) denotes allocation matrix of size \(I \times J\), which is made up of binary values. When \(x_{ij}=1\), it indicates resource \(R_i\) is allocated to zone \(Z_j\).

The constrains also specify that any resources can only be allocated to exactly \(1\) zone, and the total number of resources allocated to zone \(Z_j\) must equal to or exceeding the threshold value indicated by \(s_j\).

The objective function is a function which minimises the total cost \(z\) of the system.

\[ \begin{aligned} & \text{Let} & I & = \text{No. of resources} \\ & & J & = \text{No. of zones} \\ & & x & = \text{Allocation matrix} \\ & & c & = \text{Cost matrix} \\ & & s & = \text{Minimum no. of resources} \\ \end{aligned} \\ \\ \begin{aligned} & \text{Minimise} & z & = \sum_{i=1}^{I}\sum_{j=1}^{J}x_{ij}c_{ij} \\ & \text{Subject to} & x_{ij} & \in \{{0,1\}} & \forall &i &=1,2,3,...,I \\ &&&&& j &=1,2,3,...,J \\ & & \sum_{j=1}^{J}x_{ij} & = 1 & \forall &i &=1,2,3,...,I \\ & & \sum_{i=1}^{I}x_{ij} & \geq s_{j} & \forall &j &=1,2,3,...,J \end{aligned} \]

Problem Solving

To solve the IP problem, load the following packages in R. We will solve this problem using the simplex algorithm.

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

The IP problem is expressed in R in the following code chunk. The computation time is usually significantly longer in a complex, real-world problem.

cost_func <- function(i,j){
  vapply(seq_along(i), function(k) my_costs[i[k], j[k]], numeric(1L))
}

model <- MILPModel() %>%
  add_variable(x[i, j], i = 1:I, j = 1:J, type = "binary") %>%
  add_constraint(sum_expr(x[i, j], j = 1:J) == 1, i = 1:I) %>%
  add_constraint(sum_expr(x[i, j], i = 1:I) >= s[j], j = 1:J) %>%
  set_objective(sum_expr(x[i, j] * colwise(cost_func(i,j)), i = 1:I, j = 1:J), sense = "min")

my_solution <-  model %>% 
  solve_model(with_ROI("glpk", verbose = TRUE))
## <SOLVER MSG>  ----
## GLPK Simplex Optimizer, v4.47
## 15 rows, 50 columns, 100 non-zeros
##       0: obj =  0.000000000e+000  infeas = 1.900e+001 (10)
## *    24: obj =  5.697057522e+000  infeas = 0.000e+000 (1)
## *    45: obj =  2.058630769e+000  infeas = 0.000e+000 (0)
## OPTIMAL SOLUTION FOUND
## GLPK Integer Optimizer, v4.47
## 15 rows, 50 columns, 100 non-zeros
## 50 integer variables, all of which are binary
## Integer optimization begins...
## +    45: mip =     not found yet >=              -inf        (1; 0)
## +    45: >>>>>  2.058630769e+000 >=  2.058630769e+000   0.0% (1; 0)
## +    45: mip =  2.058630769e+000 >=     tree is empty   0.0% (0; 1)
## INTEGER OPTIMAL SOLUTION FOUND
## <!SOLVER MSG> ----

Extracting the variable \(x_{i,j}\) will return the optimal allocation results.

my_allocation_results <- my_solution %>%
  get_solution(x[i,j]) %>% 
  as_tibble() %>%
  filter(value == 1) %>%
  select(i, j) %>%
  arrange(i)

my_allocation_results
## # A tibble: 10 x 2
##        i     j
##    <int> <int>
##  1     1     2
##  2     2     4
##  3     3     3
##  4     4     3
##  5     5     1
##  6     6     3
##  7     7     1
##  8     8     5
##  9     9     4
## 10    10     2

Verifying the Solution

We can count the number of resources allocated to each zone, which is indicated in column n below. Compare that with the target threshold column s.

my_allocation_results %>%
  group_by(j) %>%
  count() %>%
  ungroup() %>%
  arrange(j) %>%
  bind_cols(s = s)
## # A tibble: 5 x 3
##       j     n     s
##   <int> <int> <dbl>
## 1     1     2     2
## 2     2     2     1
## 3     3     3     3
## 4     4     2     2
## 5     5     1     1