Get presolved model from Gurobi

Hi, in Julia is it possible to get the model after Gurobi finishes “presolve”? That is the model after being simplified by Gurobi (less variables, constraints,…).

The method Model.presolve() in the Gurobi Python interface seems to do this, but can we do it in Julia?


This is not supported by Gurobi.jl at the JuMP or MOI level.

But if you are familiar with the C API, you could use that.

using JuMP, Gurobi
model = direct_model(Gurobi.Optimizer())
# ... stuff here ...
presolvedP = Ref{Ptr{Cvoid}}(C_NULL)
GRBpresolvemodel(backend(model), presolvedP)
GRBwrite(presolvedP[], "presolve.lp")
1 Like

I tried to look into this a while ago. I came up with the following which might give some examples:

using JuMP, Gurobi
model = direct_model(Gurobi.Optimizer(GRB_ENV))
# ...

ps_model_ref = Ref{Ptr{Gurobi.GRBmodel}}()
ps_ret = Gurobi.GRBpresolvemodel(backend(model).inner, ps_model_ref)
ps_model = ps_model_ref[] # Get the GRBmodel pointer

# Direct optimization call:
# Query attributes directly:
nvars = Ref{Cint}() # Number of variables
Gurobi.GRBgetintattr(ps_model, Gurobi.GRB_INT_ATTR_NUMVARS,nvars)
nintvars = Ref{Cint}() # Number of integer variables
Gurobi.GRBgetintattr(ps_model, Gurobi.GRB_INT_ATTR_NUMINTVARS,nintvars)
nbinvars = Ref{Cint}() # Number of binary variables
Gurobi.GRBgetintattr(ps_model, Gurobi.GRB_INT_ATTR_NUMBINVARS,nbinvars)
nlincons = Ref{Cint}() # Number of linear constraints
Gurobi.GRBgetintattr(ps_model, Gurobi.GRB_INT_ATTR_NUMCONSTRS,nlincons)
ngencons = Ref{Cint}() # Number of general constraints
Gurobi.GRBgetintattr(ps_model, Gurobi.GRB_INT_ATTR_NUMGENCONSTRS,ngencons)

nnodes = Ref{Cdouble}() # Number of nodes
Gurobi.GRBgetdblattr(ps_model, Gurobi.GRB_DBL_ATTR_NODECOUNT,nnodes)

## Saving to file:
postps_filename = "my_presolved_model.mps"
postps_filepath = joinpath(@__DIR__,postps_filename)

## Reading it back:
postps_model = MOI.FileFormats.Model(format = MOI.FileFormats.FORMAT_MPS)
MOI.read_from_file(postps_model, postps_filepath)

rmodel = direct_model(Gurobi.Optimizer(GRB_ENV))
MOI.copy_to(rmodel, postps_model);

Thank you very much @odow and @jd-foster ! These are exactly what I was looking for!
Additionally, is there a way to do the following?
Say the model after presolve is:

variables: x[1:100], y[1:100, 1:50]
min x[1]+x[2]-6*y[1,2]
s.t. x[1]-2*x[3]+y[1,2] <= 1
      ...other constraints...

Can I get the coefficients and the corresponding variable indices in the constraints and objective? Also can I know if that constraint is <=, >= or ==? I would like get these in Julia vectors, matrices, dictionaries, or anything, but not as JuMP objects, so that I can play with them without JuMP.

In my example, for the constraint x[1]-2*x[3]+y[1,2] <= 1
this would be something like

LHS = [1;-2;1], 
RHS = [1], 
var = [x[1];x[3];y[1,2]], 
constraint_sense = ["<="]

and similar for the objective function.

And if I have nonlinear constraints (or objective), like x[1]-5*x[3]*y[1,2] +x[4]<= 1, can I have

LHS = [1;-5;1], 
RHS = [1], 
var = [x[1]; x[3]*y[1,2]; x[4]], 
constraint_sense = ["<="]

Anything that can achieve something like this would be great!


1 Like

Given the presolved Gurobi model, you can use the C API:

The C API can be difficult to work with if you aren’t familiar with it, so the best place to look for examples is inside the Gurobi.jl source code:

The easier option would be to read the model back into JuMP from the file:

using JuMP
presolved = read_from_file("my_presolved_model.mps")

and then use the JuMP API to get various bits:

but this will leave you with JuMP variables and constraints, not with matrices.

There is an internal function to create a “standard form” version for inspecting the problem

## JuMP.jl/src/lp_sensitivity2.jl:
# Docstring:
# Given a problem:
#     r_l <= Ax <= r_u
#     c_l <=  x <= c_u
# Return the standard form:
#            [A -I] [x, y] = 0
#     [c_l, r_l] <= [x, y] <= [c_u, r_u]

If you have an integral problem first do:

undo_relax = relax_integrality(rmodel);

then you can do

s = JuMP._standard_form_matrix(rmodel)
julia> keys(s)
(:columns, :lower, :upper, :A, :bounds, :constraints)


A = s.A

#  If n is the number of variables in the original problem
# and m is the number of affine constraints in the problem
# then
m, p = size(s.A)
# where p == n + m.

# Hence, we can recover the "original" constraint matrix as
rA = A[:,1:(p-m)]

## Restore the original integral problem if necessary with:

In the case of an existing equality constraint at row i, we have that
0 = r_l[i] <= y[i] <= r_u[i] = 0; i.e. y[i] == 0.