MILPModel and MIPModel give different objective values

[moaarab]

Hello everyone,

I am trying to get to know the OMPR package by setting up and solving some small examples. My first attempt was using the MIPModel backend which gave me the correct objective value as well as decision variable values. However when I try to replicate the model using the MILPModel backend I get different/wrong results. The code for both are the same with the only difference of using different backends. Any explanation for this?

Reproducible_Example.txt

[dirkschumacher]

Thanks! Happy to look into it.

Could you post a smaller, reproducible example? Currently your example can only run by you. A good strategy is to remove constraints and components from the objective function as long as the problem persists. Also try using artificial data not stored in external files.

[moaarab]

Hi Dirk, thanks for the quick reply. I edited my main post with a reproducible example. If you change the backend from MIP to MILP the objective value changes from 198500 to 215000 where I know the correct solution is the former.

[hugolarzabal]

Hello, Your problem is in the objective function with the 3 functions with parameters.

With MILP models, everything is calculated with vectors at the same time.

In other words, fun_Cost_Factory_To_Depot(f,d) * x[f,d] is exactly the same than fun_Cost_Factory_To_Depot(1:nbr_fac,1:nbr_dep) * x[1:nbr_fac,1:nbr_dep]

In your script fun_Cost_Factory_To_Depot(1:nbr_fac,1:nbr_dep) return [1] 0.5 0.3 1.0 0.2 which is only part of the matrix. That is why the solver doesn't find the right answer.

To solve this problem there is an ompr function named colwise. You can either use it in the function fun_Cost_Factory_To_Depot:

fun_Cost_Factory_To_Depot  <- function(f, d) {
  mapply(function(x, y) {
    Cost_Factory_To_Depot[x, y]
  }, f,d)%>%colwise()
}

or in the objective function:

set_objective(
    sum_expr(colwise(fun_Cost_Factory_To_Depot(f,d)) * x[f,d], f = 1:nbr_fac, d = 1:nbr_dep) + 
      sum_expr(colwise(fun_Cost_Factory_To_Customer(f,c)) * y[f,c], f= 1:nbr_fac,  c = 1:nbr_cust) + 
      sum_expr(colwise(fun_Cost_Depot_To_Customer(d,c)) * z[d,c], d= 1:nbr_dep,  c = 1:nbr_cust), 
    sense = "min")

[dirkschumacher]

I need to find a solution for colwise - it seems to be a source of potential errors. @hugolarzabal thanks so much for answering 👍

[moaarab]

Thank you for the answer @hugolarzabal ! Your hint was indeed the key to solving my problem. Does this mean that every time I define a function with parameters I have to the use the colwise() function? What about when there is only one parameter in the function e.g. f(x) instead of f(x,y)?

[hugolarzabal]

If it is just f(x), you won't need colwise. You won't need a function either since the vector will work). For example, CustomerDemand[c] works fine in your script.

For f(x,y), I think you will almost always need colwise.

[moaarab]

Thanks again @hugolarzabal for the swift and clear reply. I am already trying it out on a new model I am building. When I add following constraint:

 add_constraint(
    (dummy[r,r_help]*fun_distance_ret_ret(r,r_help)) >=  5*dummy[r,r_help], r=1:(nbr_pos_ret-1), r_help=2:nbr_pos_ret, r_help>r
  )

I get the following error:

Error in `[[<-.data.frame`(`*tmp*`, "coef", value = c(14, 10, 15, 4, 19,  : 
  replacement has 1500625 rows, data has 1225

The 1500625 is the 1225 squared. The colwise() function was applied when defining the fun_distance_ret_ret function. Any idea what the reason for this error could be? I use another function in the model elsewhere without any issues, however that one is within a sum_expr().

[hugolarzabal]

I think that the problem may be in the function. Could you post the function fun_distance_ret_ret here too ?

Is nbr_pos_ret equal to 50 ?

[moaarab]

The function in which distance_ret_ret is a matrix:

fun_distance_ret_ret  <- function(r,r_help) {
  mapply(function(x, y) {
    distance_ret_ret[x, y]
  }, r,r_help) 
} %>%colwise()

Yes nbr_pos_ret is equal to 50. Since I am only interested in the constraints for indices where r_help > r, I should have (50*49)/2 = 1225 constraints.

[hugolarzabal]

It just occured to me that your constraint is pointless since you have the same variable in both ends of the equation. It is the same as

fun_distance_ret_ret(r,r_help)) >= 5

It will either be always true or always infeasible.

I don't know if the problem comes from there but that may be part of the problem.

[moaarab]

Actually I forgot to mention that the binary variable dummy[r,r_help] is equal to the product of two other binary variables namely x[r] and x[r_help] for r_help > r. Since the product of two binary variables cannot enter a linear optimization problem I had to use a "trick" to transform that product into a new binary variable (for more info: https://www.leandro-coelho.com/linearization-product-variables/). Which means that that particular constraint will only be "activated" when both x[r]and x[r_help] are both equal to one, exactly as I want it.

When I try to run the model using the MIP backend it works correctly. So I also think it has something to do with the function.

[hugolarzabal]

You may try writing your constraint that way:

dummy[r,r_help] <= (fun_distance_ret_ret(r,r_help)) >= 5)

(fun_distance_ret_ret(r,r_help)) >= 5) is a condition that may be true(1) or false(0) and will force the variable dummy to be equal to 1 or 0.

It may be more elegant to put the condition directly in the function in order to have the function return 1 or 0.

Since the numbers 1500625 and 1225 suggest that the model try to combine the combinations of indexes too many times, removing one ocurence of the variable may solve the problem.

[moaarab]

Your reformulation of the constraint does indeed work as a charm! I have also found out that I can get the original constraint to work when I do not apply the%>%colwise()function on fun_distance_ret_ret(). I think it is best for me personally to first build a small example of the problem using the MIP backend and MILP backend and only continue using the MILP backend when I manage to reproduce the same results with both backends. Thank you again both for the help.

[dirkschumacher]

Thanks from me as well. I will need to work more on the ´colwise´thing to not cause so much confusion. ´MILPModel´still is WIP.

[rickyars]

@hugolarzabal can you give an example using the function notation and colwise with a variable that has 3 or more indices? I often have code with i's, j's and k's but I'm stuck using the slower MIP back-end. Being able to switch would save me probably half my run-time. Which on my 8 core machine is already several days. Thanks!

[hugolarzabal]

This function will allow you to transform a dataframe with your values into a function that return the right values for i, j and k

# .df is a dataframe with the following columns : values of i, values of j, values of k, result in that order
param <- function(.df,i,j,k,result){

  .df[(.df[[1]] %in% i) & (.df[[2]] %in% j) & (.df[[3]] %in% k),][[result]] %>% colwise()

}

If you have the values in a dataframe called my_df, you will use the function like this: add_constraint( v[i,j,k] * param(my_df, i,j,k,"result") < 1, i = 1:n, j = 1:m, k = 1:p)

[rickyars]

@hugolarzabal i'll give it a shot, thank you!

[Samidha30]

Hi everyone, re-opening this chain as I am working on something similar and need some help with framing nested constraints for 3D variables in MILP. I am trying to optimize a supply network with m customers, n DCs, k plants and q SKUs. The constraint to equate plant outflow with customer demand by each DC and SKU is not working properly.

add_constraint(sum_expr(aa[l,j,q], l = 1:k) == sum_expr(VCustDemand[i+(q-1)m] x[i,j,q], i = 1:m), j = 1:n, q = 1:r)

Also attaching a reproducible version of the complete R code for details.

sample code.txt

[dirkschumacher]

@Samidha30 can you try to reduce your problem further? Ideally a model with as few variables and constraints as possible but still reproducing the problem. See also the reprex package that helps generated reproducible code chunks.