I am trying to write an MILP to do a grouping (clustering) given one difference condition.

Suppose, I have the following elements in a list: 1,2,3,4,6,8,9,10.

My goal here is to group them so that the difference between two elements won't cross a given threshold for each group and the number of groups should be minimised. So, if I give a difference threshold of 6, my two groups should be: 1,2,3,4,6 and 8,9,10. (1 & 8 can't be in the same group because the difference is 7)

To formulate the MILP, I am considering the following:

Objective function: Minimise the number of groups


  1. Each element must be in exactly one group
  2. If a group is used then it should be considered
  3. Two elements of the same group shouldn't have a difference greater than the given threshold.

I have modelled this in R like this but it is converging very slowly and often the output is not correct. Is there any more efficient way to model this?


# Define your element ids and the difference
element_ids <- c(1, 2, 3, 4, 5, 6, 7, 9, 10, 12)
diff <- 5

model <- MIPModel() %>%
  # Binary variable for each element id and group
  add_variable(x[i, g], i = 1:length(element_ids), g = 1:length(element_ids), type = "binary") %>%
  # Binary variable for each group
  add_variable(y[g], g = 1:length(element_ids), type = "binary") %>%
  # Objective: minimize the number of groups
  set_objective(sum_expr(y[g], g = 1:length(element_ids)), "min") %>%
  # Constraint: each element id must be in exactly one group
  add_constraint(sum_expr(x[i, g], g = 1:length(element_ids)) == 1, i = 1:length(element_ids)) %>%
  # Constraint: if a group is used, then y[g] = 1
  add_constraint(sum_expr(x[i, g], i = 1:length(element_ids)) <= length(element_ids) * y[g], g = 1:length(element_ids))

# Constraint: no two element ids in the same group can have a difference equal to 'diff'
for (i in 1:(length(element_ids) - 1)) {
  for (j in (i + 1):length(element_ids)) {
    if (element_ids[j] - element_ids[i] > diff) {
      model <- model %>%
        add_constraint(x[i, g] + x[j, g] <= 1, g = 1:length(element_ids))

# Solve the model
result <- solve_model(model, with_ROI(solver = "glpk", verbose = TRUE))

# Print the solution
solution <- result %>% get_solution(x[i, g]) %>% filter(value > 0) %>% select(i, g, value)

i g value
1   1   1       
2   1   1       
3   1   1       
4   1   1       
5   1   1       
6   1   1       
7   10  1       
8   10  1       
9   10  1       
10  10  1   

  • 1
    $\begingroup$ It seems to me that your model doesn't make use of the 1-dimensional nature of the problem? For instance if 1 and 3 are in the same group, then of course 2 is in that group too. In addition, you already have a good upper-bound on the number of groups. So the only question is "where to make the cuts", which seems a lot easier than a general clustering problem. $\endgroup$
    – Stef
    Commented Jan 18 at 17:34
  • 1
    $\begingroup$ For instance if the smallest element is $x_0$ and the threshold is $t$, and all you care about is minimising the number of clusters, then there is no reason not to put every element between $x_0$ and $x_0+t$ in the same cluster. And then if the smallest element $> x_0+t$ is $x_1$, then there is no reason not to put every element between $x_1$ and $x_1+t$ in the same cluster. Etc. So here you have a linear-time algorithm that solves your problem, much faster than integer programming can hope to be. $\endgroup$
    – Stef
    Commented Jan 18 at 17:50

3 Answers 3


Some of this may be solver dependent. I just ran your model with only one change -- switching the solver from glpk to CPLEX -- and got a valid solution (the same as yours, changing group 10 to group 2) in 0.06 seconds.

There are tweaks you can make to your model that may or may not improve execution time. It contains a fair bit of symmetry, since group indexing is arbitrary. (Permuting group indices in a solution produces another valid solution.) You can reduce symmetry either by using a parameter setting (if one exists) telling your solver to work on symmetry reduction or by adding the following constraint.

add_constraint(y[g] <= y[g - 1], g = 2:length(element_ids))

Also, your constraint "if a group is used, then y[g] = 1" can be disaggregated, leading to more but less dense constraints, by changing it to the following.

add_constraint(x[i, g] <= y[g], i = 1:length(element_ids), g = 1:length(element_ids))

Whether either of those tweaks helps is an empirical question.

As to the output being incorrect, if you are referring to having group numbers 1 and 10 rather than 1 and 2, that is addressed above. If you mean that some instances produce infeasible or suboptimal solutions, I suggest you start a new question and post the source code for such an instance. Your model is mathematically correct, so infeasible/suboptimal solutions should not occur.

  • $\begingroup$ Thank you! These two tweaks did help! But still, when the length of element_ids is increasing, it is converging very slowly. Took 11 secs for 25 elements, and it has been more than 2 mins, and haven't found the optimal solution for 30 elements. But I guess this might just be a solver problem, not a math problem. $\endgroup$ Commented Jan 19 at 4:24

Since you've a threshold (diff) you may pre-define max of diff number of groups, not the full length of the list. As for the if condition try this
$L(x_{i,g}+x_{j,g}-1) \le d - (L_j - L_i)$
where $L$ is the largest number is the sorted list.
Try sorting the list first (outside of the model)

  • 3
    $\begingroup$ In fact, the whole algorithm can be run outside of an integer programming model: (1) Sort the list. (2) Assign the smallest element to the first group. Assign all elements who are $\leq$ (smallest element + threshold) to the first group. (3) Assign next unassigned element to second group. Assign all next elements who are smaller than this element + threshold to the second group. (4) Etc etc. $\endgroup$
    – Stef
    Commented Jan 18 at 17:42
  • $\begingroup$ Thanks @Stef , I might be over-engineering this. Do you think your proposed method will give optimal solution every time? Just thinking out loud, if it will fail (sub-optimal) anytime or not. $\endgroup$ Commented Jan 19 at 4:54

Used the approach suggested in the comment to solve this:

element_ids <- c(1, 2, 3, 4, 5, 6, 7, 9, 10, 12)
diff <- 5
element_ids <- sort(element_ids)

# Initialize an empty list to hold the groups
groups <- list()

# While there are still elements left to assign
while (length(element_ids) > 0) {
  # The first element is the smallest one
  first_element <- element_ids[1]
  # Find all elements that are less than or equal to first_element + diff
  group <-
    element_ids[element_ids <= first_element + diff]
  # Add this group to the list of groups
  groups[[length(groups) + 1]] <- group
  # Remove the elements in this group from element_ids
  element_ids <- setdiff(element_ids, group)


[1] 1 2 3 4 5 6

[1]  7  9 10 12

It is quite fast and gives "good enough" answer. Though the approach suggested by prubin also helped.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.