sum_expr with index dependent on first index (object not found) & MIP to MILP constraint transformation

[tjseydel]

First off all, thank you for creating a terrific intuitive package for modelling MILP's.
I'm modelling a scheduling problem and I have a question concerning the use of sum_expr.

In my model I'm trying to model the following constraint:

I tried to model the constraint as follows:
(this is a simplified example that does not account for negative values in t_t)

timeslots <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17)
r <- c(1,2,3,4,5,6,7,8,9,10)
add_constraint(sum_expr(x[r,t_t], 
            r = requests, 
            t_t = t-l[r,"Duration"]+1), 
            <= chairs, t = timeslots)

l is a dataframe containing the requests (r) and duration of the request, respectively.
The variable X[r,t_t] is 1 when a request starts at time t_t and 0 otherwise.
The idea behind the constraint is that it ensures that the number of chairs taken by the requests that start on any time t and the requests that started before t (but have not ended yet) do not exceed the number of available chairs.

For this to work t_t is dependent on the current time (t) and the duration of each individual request (r).
Though when running the code I get the error message of object 'r' not found. I couldn't find a way to solve my problem in previously posted issues or documentation.
I hope someone is able to help me out.

[hugolarzabal]

I think you can do it with a filter:
add_constraint(sum_expr(x[r,t_t], r = requests, t_t = timeslots, filter = f(t, request))
<= chairs, t = timeslots)

The function f should return a vector with True and False for every combination of t_t and r for a given t.

[tjseydel]

Thank you for your quick response hugolarzabal!
I must say I'm not familiar with using a function as a filter for sum_expr.
Just to be clear, should the function f return a vector with multiple elements (like TRUE TRUE FALSE FALSE) or should it just return a TRUE or FALSE for every combination of t_t and r for a given t whenever it is called?
Thanks again for trying to help me out

[hugolarzabal]

The function should return a vector with a TRUE or FALSE for every combination of t_t and r for a given t.

sum_expr doesn't always work very intuitively so it is not an easy function to write.

Maybe it would be simpler for you to use an auxiliary variable telling if a chair is occupied or not because of a request. Then you will be able to separate your constraint into 2 simpler constraints.

(1) y[r, t] >= x[r, t_t] for every t_t in [ t, t + L[r] ] for every r

(2) sum(y[r,t], r = request) <= chair for every t 

You will still have to implement the filter in (1) but I think it will be simpler in a regular constraint than in a sum_expr. You may not even need a function.

[tjseydel]

Thanks again for the fast response and the idea of using auxiliary variables.

I've tried to construct a function that would return a TRUE or FALSE for every combination of r and t_t for a given t, though I can't seem to figure out why it won't work.

If I can't make the sum_expr work, I will try the method with the auxiliary variables. Though I hope I can make the sum_expr work since it would probably take less computational effort.

In my formula I want to pass the t_t, r and t [so f(t_t,r,t)]. The model doesn't seem to have any trouble to pass the r and t to the function but doesn't want to pass the t_t resulting in an error when trying to construct the model.

Do you have any idea why the model can pass the r and t but not the t_t?

[hugolarzabal]

Could you provide a reproductible example for that problem ?

It is quite hard to find a solution with only a description of the error.

[tjseydel]

In essence, the following code describes what I want to accomplish. I left out all other constraints. When the function is loaded into the model it returns an error: R Error in f(t_t,r,t): object 't_t' not found.
When I pass a function with only t and r [f(t,r)] it doesn't give an error.
Thanks again for trying to help me!

require(ompr)
require(ompr.roi)
require(ROI.plugin.glpk)
require(knitr)
require(dplyr)

timeslots <- 1:17
chairs <- 14 
requests <- seq(from= 1, to = 25, by= 1)
l<- data.frame(Treatmentcode= sample(1:47,length(requests)),Duration= replicate(length(requests),sample(1:12,1)))

f <- function(t_t,r,t){
  if (t_t >= t-l[r,"Duration"]+1){
    if (t_t <= t){
      return(TRUE)
    }
  }
  return(FALSE)
}

model <- MIPModel() %>%
  #Create a variable that is 1 if request r is planned to start on time t
  add_variable(x[r,t], r= requests, t= timeslots,
               type = "binary") %>%
  
  #Create S_End variable which correspsonds to the start time of the dummy variable End, can not be negative
  add_variable(S_End, type= "integer", lb=0, ub = 19) %>%
  
  #Minimize start time of End (makespan)
  set_objective(S_End,"min") %>%

  #Constraint in which the t_t variable gives an error when passed to the f function. 
  add_constraint(sum_expr(x[r,t_t], r = requests, t_t = timeslots, filter = f(t_t,r,t)) <= chairs, t = timeslots)

[hugolarzabal]

Putting the filter that way instead of using a function seems to work:

add_constraint(sum_expr(x[r,t_t], r = requests, t_t = timeslots, (t_t >= (t-l[r,"Duration"]+1)) | (t_t <= t) ) <= chairs, t = timeslots)

[tjseydel]

This is exactly what I meant. The constraint makes so much sense now I see it. I only replaced the OR (|) operator for an AND (&) operator. Thanks again for helping me out!

[tjseydel]

I updated my ompr package to the master version so I could better update the MIP model to a MIPL model to increase computational speed. So far I was able to alter almost all of the constraints so they would work with a MILP model. You can probably guess that constraint discussed in this issue is causing me some trouble again.

First of all, since the update to the master version the constraint in the MIP model did not work anymore. I figured that it was because the model only took one value for r in l[r,"Duration"].
I solved this by:

add_constraint(sum_expr(x[r,t_t], r = requests, t_t = timeslots, (t_t >= (t-l[requests,"Duration"]+1)) & (t_t <= t) ) <= chairs, t = timeslots)

After adding the constraint I added an additional variable (C_Max) instead of the chairs parameter such that I could minimize the number of total chairs used.

Though, when directly using this constraint in a MILP it gives the warning message:
In t - template_requests[requests, "durations"] : longer object length is not a multiple of shorter object length

I'm not able to produce a simple reproducible example. I used the following code (which is a simplified version of my whole model and added a small check at the end which produces an overview of when which request is scheduled):

require(ompr)
require(ompr.roi)
require(ROI.plugin.glpk)
require(knitr)
require(dplyr)

#Set allowed number of patients that can start therapy at respectively: 
n <- c(2     ,2    ,2    ,2    ,2    ,1    ,2    ,2    ,1    ,2    ,1    ,2    ,2    ,2    ,2    ,2    ,2   )
#08:00,08:30,09:00,09:30,10:00,10:30,11:00,11:30,12:00,12:30,13:00,13:30,14:00,14:30,15:00,15:30,16:00)

#Create data set and set parameters
template_requests <- data.frame(request = 1:21, durations = c(2,4,6,2,2,6,10,16,6,10,2,2,2,2,2,4,2,2,2,2,2))
requests <- template_requests$request
timeslots <- seq(1,17,1)


#Create a makespan matrix and funtion that gives a vector of the makespan for each value of X given request r and timeslot t
Makespan <- matrix(0, nrow = length(requests), ncol = length(timeslots))
for (r in 1:length(requests)){
  for (t in 1:length(timeslots)){
    Makespan[r,t] <- t + template_requests[r,"durations"]-1
  }
}

Makespan_fun <- function(r,t){
  Makespan_sub <- t(Makespan[r, t]) 
  Makespan_vec <- as.vector(Makespan_sub) 
  return(Makespan_vec)
}

model <- MILPModel() %>%
  #Create a variable that is 1 if request r is planned to start on time t
  add_variable(x[r,t], r= requests, t= timeslots,
               type = "binary") %>%
  
  #Create S_End variable which correspsonds to the start time of the dummy variable End, can not be negative
  add_variable(S_End, type= "integer", lb=0, ub = 19) %>%
  
  #Create C_Max variable which corresponds to the maximum numbers of chairs used on a day. The variable is constrained between 0 and 14 (14 is the chair capacity)
  add_variable(C_Max, type= "integer", lb=0, ub = 14) %>%
  
  #Minimize start time of End (makespan) and within that solution minimize the number of chairs used
  set_objective(14*S_End + C_Max,"min") %>%
  
  #Constraint 1: Each activity starts exactly once
  add_constraint(sum_expr(x[r,t],t = timeslots) == 1, r = requests) %>%
  
  #Constraint 2: The number of patients that can start treatment at a given time is respected
  add_constraint(sum_expr(x[r,t], r = requests) <= n[t], t = timeslots) %>%
  
  #Constraint 3: The start time of End is bigger or equal to the end of the last patient on a scheduled day
  add_constraint(sum_expr(x[r,t]*colwise(Makespan_fun(r,t)), t = timeslots) <= S_End, r = requests) %>%
  
  #Constraint 4: The maximum number of chairs used is bigger or equal to the sum of the number of chairs at any given time t
  add_constraint(sum_expr(x[r,t_t], r = requests, t_t = timeslots, (t_t >= (t-template_requests[requests,"durations"]+1)) & (t_t <= t) ) <= C_Max, t = timeslots)

#Solve model
result <- solve_model(model,with_ROI(solver = "symphony")) 
result$solution[1]
result$solution[2]

#Collect solution
matching <- result %>% 
  get_solution(x[r,t]) %>%
  filter(value > 0) %>%  
  dplyr::select(r, t)
for (i in 1:length(requests)){
  matching[i,"Duration"] <- template_requests$durations[matching[i,"r"]]
}

#Check solution
Roster <- array(0,dim=c(1,14,19))

for (i in 1:dim(matching)[1]){
  spot <- FALSE
  c <- 1
  while (spot == FALSE & c<= 14) {
    if(Roster[1,c,matching[i,"t"]]==0){
      start <- matching[i,"t"] 
      end <- start + matching[i,"Duration"] - 1
      for (j in start:end){ #-1 because it starts at 0 instead of 1
        Roster[1,c,j] <- matching[i,"r"]
      }
      spot <- TRUE
      break
    }
    c <- c + 1
  }
  if(c > 14){print(paste0("too many patients on chairs at time ",  matching[i,2]))}
}

Roster[1,,]

The constraint doesn't work anymore since the result shows C_Max = 1 even though MIP and the roster of the #check solution part clearly show more chairs are being used.

@hugolarzabal Do you have an idea of how I could transform constraint 4 to work with a MILP model? I tried looking into the model structure itself, but the MILP is a lot harder to understand than a MIP model.

[hugolarzabal]

The function should return a vector with a TRUE or FALSE for every combination of t_t and r for a given t.

sum_expr doesn't always work very intuitively so it is not an easy function to write.

Maybe it would be simpler for you to use an auxiliary variable telling if a chair is occupied or not because of a request. Then you will be able to separate your constraint into 2 simpler constraints.

(1) y[r, t] >= x[r, t_t] for every t_t in [ t, t + L[r] ] for every r

(2) sum(y[r,t], r = request) <= chair for every t 

You will still have to implement the filter in (1) but I think it will be simpler in a regular constraint than in a sum_expr. You may not even need a function.

I think that would be my best answer with a MILPModel.

The sum_expr function is quite unintuitive and creating a filter which use indexes from inside and outside the sum_expr doesn't work well.

Still, using an auxiliary variable can improve the readability of you model and may not necessarily increase solver time. Sometime, it may even be easier for the solver to use the auxiliary variable than to work with a complicated expression.