Why is Gurobi returning MIP solution with gap of 0 which is not close to the actual optimal solution?
Ongoing
Using the below code and the following two .csv files, Gurobi returns an optimal solution of 62.42 when the parameter SET_SOME_VARS is set to 0, but returns a much better solution (of 12.32) when it is set to 1. But setting it to 1 is just forcing one feasible solution for the problem, so why is the version without these constraints finding a much worse "optimal" solution? The MIP solution ends with a gap of 0 in both cases.
Thanks in advance for any insights on this!
(Census_Population_1970.csv can be obtained via this google sheet):
https://docs.google.com/spreadsheets/d/1yS7FWRrknAEUnQzDXu1Fv0eiY2UplY_ElQleXcyMqOQ/edit?usp=sharing
(moo93.csv can be obtained via this google sheet):
https://docs.google.com/spreadsheets/d/1NkuxVslk9PhpFv0eujiNvnOaBK7nWf6ekexn3dkcsjk/edit?usp=sharing
library(readxl)
library(writexl)
library(Matrix)
library('gurobi')
setwd("C:/Users/Steven Shechter/Dropbox/R/Gurobi/Apportionment_Gurobi")
SET_SOME_VARS <- 0 # setting to 1 forces some of the decision variables to take on certain pre-set value
# setting to 0 lets Gurobi optimize these and the other variables
total_reps <- 435
num_states <- 50
obj_weight_pairwise <- 94
obj_weight_range <- 100-obj_weight_pairwise
data <- read.csv("Census_Population_1970.csv")
total_pop <- sum(data$Population)
data$quota <- data$Population/total_pop*total_reps # the # of seats the state should get under continuous reps by
###########################
# SET UP DECISION VARIABLES
###########################
num_x_vars <- num_states
num_v_vars <- (num_states-1)*num_states/2
num_y_vars <- 1
num_vars <- num_x_vars + num_v_vars + num_y_vars
model <- list()
#setting up the x vars
for (i in 1:num_x_vars)
{
model$vtype[i] = 'I'
model$varnames[i] = paste0('X', i)
model$obj[i] = 0
model$lb[i] = 1
}
#setting up the v vars
for (i in 1:num_v_vars)
{
var_index <- num_x_vars + i
model$vtype[var_index] = 'C'
model$varnames[var_index] = paste0('V', i)
model$obj[var_index] = obj_weight_pairwise #the weight given to the v vars
model$lb[var_index] = 0
}
# setting up the y vars
for (i in 1:num_y_vars)
{
var_index <- num_x_vars + num_v_vars + i
model$vtype[var_index] = 'C'
model$varnames[var_index] = paste0('Y', i)
model$obj[var_index] = obj_weight_range #the weight given to the y var
model$lb[var_index] = 0
}
###########################
# SET UP CONSTRAINTS
###########################
# Constraint 1: sum of x vars = 435
model$A <- matrix(c(rep(1, num_x_vars), rep(0, num_v_vars), rep(0, num_y_vars)), nrow = 1, byrow = FALSE) # constraints of type labeled constraint 1 above (which is just 1 constraint)
model$sense <- '='
model$rhs <- total_reps
#TESTING: activating the following code chunk returns a feasible solution with a better objective value (12.32),
# compared to the objective value of 62.42 the optimization returns when it is de-activated
if(SET_SOME_VARS) {
# Constraints to fix x vars at certain values given in the following .csv file (obtained by a weighting of 93 and 7 instead of 94 and 6)
moo93 <- read.csv("moo93.csv")
B <- spMatrix(num_x_vars, num_vars,
i = (1:num_x_vars),
j = (1:num_x_vars),
x = rep(1, num_x_vars))
model$A <- rbind(model$A,B)
model$sense <- c(model$sense, rep('=', num_x_vars))
model$rhs <- c(model$rhs, moo93$MOO93)
}
# Constraints to set up the v and y variables to take on the intended objective values
v_index <- 1
qc_index <- 1
cumulative_quadcon <- list()
temp_list <- list()
for (i in 1:(num_states-1))
{
for (j in (i+1):num_states)
{
temp_list$Qc <- spMatrix(num_vars, num_vars, c(i, j), c(i, j), c((data$quota[i])^(-2), -(data$quota[j])^(-2)))
temp_list$q <- sparseVector(c(-(data$quota[i])^(-2), -(data$quota[j])^(-2), -1), c(i,j,num_x_vars+v_index), num_vars)
temp_list$rhs <- 0
cumulative_quadcon[[qc_index]] <- temp_list
qc_index <- qc_index + 1
temp_list$Qc <- spMatrix(num_vars, num_vars, num_vars, i, -1/data$quota[i])
temp_list$q <- sparseVector(1/data$quota[j], j, num_vars)
temp_list$rhs <- 0
cumulative_quadcon[[qc_index]] <- temp_list
qc_index <- qc_index + 1
temp_list$Qc <- spMatrix(num_vars, num_vars, c(i, j), c(i, j), c(-(data$quota[i])^(-2), (data$quota[j])^(-2)))
temp_list$q <- sparseVector(c(-(data$quota[i])^(-2), -(data$quota[j])^(-2), -1), c(i,j,num_x_vars+v_index), num_vars)
temp_list$rhs <- 0
cumulative_quadcon[[qc_index]] <- temp_list
qc_index <- qc_index + 1
temp_list$Qc <- spMatrix(num_vars, num_vars, num_vars, j, -1/data$quota[j])
temp_list$q <- sparseVector(1/data$quota[i], i, num_vars)
temp_list$rhs <- 0
cumulative_quadcon[[qc_index]] <- temp_list
qc_index <- qc_index + 1
v_index <- v_index + 1
}
}
model$quadcon <- cumulative_quadcon
model$modelsense <- 'min'
params<-list()
params$NonConvex <- 2
# solve model
result <- gurobi(model, params)
if (result$status == 'OPTIMAL') {
data$reps <- result$x[1:num_states]
} else {
cat('Optimization finished with status',result$status)
}
-
Hi Steven,
Could you use the gurobi_write() function to write 2 LP files (one with SET_SOME_VARS=0 and the other with =1) and upload those? This would make helping you easier.
Best regards,
Jaromił0 -
Thank you, Jaromil. I don't see how to upload a file, but here are links to google docs with the .lp outputs. These are long files, but I ran a "compare" of these in Word and confirmed that the only differences between the files are that the first one explicitly forces the integer variables to take on specific pre-set values (shown in the "differences" doc.
when setting SET_SOME_VARS <- 1:
https://docs.google.com/document/d/1mqcMFi2QFqCijP7YY6mSOSGvVhA6hzBH3WZ3U3AmTlw/edit?usp=sharing
when setting SET_SOME_VARS <- 0:
https://docs.google.com/document/d/1hjtd7vec5-Tka0uu10cCcPphxpy8Z2S2ZLBIxP-k5pI/edit?usp=sharing
The differences between these two files results in:
https://docs.google.com/document/d/1hjtd7vec5-Tka0uu10cCcPphxpy8Z2S2ZLBIxP-k5pI/edit
0 -
sorry, the differences link should instead be:
https://docs.google.com/document/d/1NePZgZ0PcSFUG8Mwh0byCEQAs0KNJIKLEdA1iR9g5GA/edit?usp=sharing
0 -
I observed a few cases of incorrect solutions with Gurobi 9.1.0-9.1.2, similar to what you described. However, I couldn't replicate this behavior with Gurobi 9.5.0. Which version of Gurobi are you using? Could you please post the log output printed by Gurobi when \( \texttt{SET_SOME_VARS} \) equals \( 0 \)?
0 -
Thank you, Eli! That did the trick. I was using version 9.1.2. The problem no longer appears after I updated to version 9.5.0. Was the problem in 9.1.2 related to my model's use of non-convex quadratic constraints?
thanks again!
0 -
That's great to hear! Yes, your intuition is correct - the problem in Gurobi 9.1.2 was related to the non-convex quadratic constraints.
0 -
Hi Eli,
I'm actually running into similar problems again as described above, even for the 9.5.0 version of Gurobi.
Same situation, where the "optimal" solution returned by Gurobi is worse (by two orders of magnitude) than another feasible solution obtained after fixing several decision variables. I set MIPGap to 0 in both cases (and the log files show the problems terminating with gap 0.0000%). These are also again occurring in a non-convex quadratically constrained problem.
The log files are here:
Solution without pre-fixing decision vars:
https://docs.google.com/document/d/1kvomcFwnFLMfersRXSxGOytictNTzAi6yP_KobUh-Fk/edit?usp=sharing
Objective value (for a minimization problem): 85.7
Solution with pre-fixing several decision vars:
https://docs.google.com/document/d/1s0y34J8LFw9OtuL2-oPOXKsCFbUX1NAw9uG2Dx5euOc/edit?usp=sharing
Objective value: 0.408
thanks in advance for any insights you may have on this!
best,
Steven
0 -
Could you post a link to MPS model files for the two models that produce these different results? You can create these model files using the gurobi_write() function, as you did before. Thanks!
0 -
Sure.
The MPS file for the version that fixes the 48 integer variables at pre-determined values is here:
https://drive.google.com/file/d/1DJKiieKbyhwBoda0LBkumJAXZiEQC5su/view?usp=sharing
and the MPS file that does not fix the 48 integer variables is here:
https://drive.google.com/file/d/1sAeDBgz--LawVRphs267yba7pj4izQ3H/view?usp=sharing
thanks again for having a look!
Steven
0 -
Hi Eli,
I just wanted to follow up on this and see if you were able to find out what the issue might be?
Thank you!
Steven
0 -
Sorry, Steven. I completely missed your previous post.
I could reproduce the issue with Gurobi 9.5.1. I will open a ticket in our support portal so we can investigate this further. Thanks!
0 -
I have the same issue with Gurobi 10.0.1.
The solver shows a gap = 0%, but the solution is not optimal.
I am sure the solution is not optimal by fixing the following variablesf1.irr[0,2,0] 50
f1.irr[0,2,1] 60f1.i_area[0,0,0] 0
f1.i_area[0,1,0] 0
f1.i_area[0,2,0] 1
f1.i_area[0,3,0] 0and the model returns a better objective value.
I provide the corresponding .lp, .mps, and .sol here.
https://1drv.ms/f/s!AvYlq7fSX1V1oZBXHyrg58LqTyNcVQ?e=Rdxl22
Example_org is the model that produces a non-optimal solution while claiming gap = 0.
Example_fix_decision_variables_irr is the one I fixed the above variables and got a better (higher) objective with also gap = 0.Please let me know how I can move forward to solve the issue.
0 -
When running model Example_fix_decision_variables_irr.lp, you should see that Gurobi prints a warning at the end of the solution process
Optimal solution found (tolerance 1.00e-04)
Warning: max constraint violation (1.2500e-01) exceeds tolerance
(model may be infeasible or unbounded - try turning presolve off)
Best objective 3.255595416759e-02, best bound 3.255595416759e-02, gap 0.0000%When you run model Example_fix_decision_variables_irr.lp with Presolve=0, Gurobi declares the model to be infeasible. From here, you can proceed as described in the Knowledge Base article How do I determine why my model is infeasible?
0 -
Hi, I am facing a similar issue with Gurobi 10.0.1.
By considering a minimization problem, I observed the following solutions:
- 1) with MIPGap=1%, in 21 seconds I found a solution with gap=0.9898% and objective= 70,390,180.4
- 2) with MIPGap=0, in 352 seconds I found a solution, with gap=0.0000% and objective = 70,605,324.2
(Note that in the 2nd case there is no warning about constraint violation).
How is this possible? Could the use of Heuristics=1 cause this result?
By removing Heuristics=1,
- 3) with MIPGap=0 and TimeLimit=353, I found a solution with gap =0.7871% and objective= 70,251,694.8
Here are the logs:
https://drive.google.com/drive/folders/1zHyxZCaYJ1DEWSFZCrjejB3V6g5Vb8yR?usp=sharing
Thank you in advance for the help,
Martina
0 -
Hi Martina,
How is this possible? Could the use of Heuristics=1 cause this result?
The Heuristics parameter should not have any effect on the final best bound.
Just from the logs I cannot tell what is going wrong. This may be a numerical issue and you might want to experiment with the NumericFocus parameter. It looks like you do not provide bounds for your variables. Could you try providing finite (and in best case tight) bounds for all variables participating in bilinear terms? It is possible that the unboundedness of bilinear terms is causing some numerical trouble. However, I am only guessing here.
Could you please share this particular model so I can have a closer look?
Best regards,
Jaromił0 -
Hi Jaromił,
thank you for your answer.
I set the bounds to the variables participating in bilinear terms by adding additional constraints to the model. Can you suggest a better way to provide these bounds?
The only way I can share the code with you is to send it under confidentiality to your private e-mail because it is part of my doctoral work and I have yet to publish it. Is it possible to proceed in this way?
Best regards,
Martina
0
Please sign in to leave a comment.
Comments
16 comments