Example for solving a QP (or LP) in Gurobi directly (without JuMP)

Similar to this question, I am wondering if anyone has example Julia code that they would be willing to share for how to pass vectors and matrices directly into Gurobi.jl to construct/solve a model. The Gurobi C interface qp example looked a bit limited.

The C example is almost exactly what you’re after. Here’s the equivalent of the first few lines:

using Gurobi
env_p = Ref{Ptr{Cvoid}}()
error = GRBloadenv(env_p, "qp.log")
env = env_p[]
model_p = Ref{Ptr{Cvoid}}()
error = GRBnewmodel(env, model_p, "qp", 0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
model = model_p[]
error = GRBaddvars(model, 3, 0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
qrow = Cint[0, 0, 1, 1, 2]
qcol = Cint[0, 1, 1, 2, 2]
qval = Cdouble[1, 1, 1, 1, 1]
error = GRBaddqpterms(model, 5, qrow, qcol, qval)

The things to pay attention to are:

  • Indices are Cint, not Int
  • Indices are 0-indexed, not 1-indexed

Both points together mean you need qrow = Cint[0, 0, 1, 1, 2], not qrow = Int[1, 1, 2, 2, 3].

1 Like

Working toward completeness, here is a script that solves the problem with constraints. To play with setting variable bounds, I removed the implicit [0,+\infty) bounds and added them explicitly later. One remaining question I have is how to access the variables by name with GRBgetvarbyname (see second code block). How should I name the constraints / variables so that I can get their values by name after optimizing?

#=
minimize    x^2 + x*y + y^2 + y*z + z^2 + 2 x
subject to  x + 2 y + 3 z >= 4
            x +   y       >= 1
            x, y, z non-negative
=#

using Gurobi

# initialize model
env_p = Ref{Ptr{Cvoid}}()
error = GRBloadenv(env_p, "qp.log")
env = env_p[]
model_p = Ref{Ptr{Cvoid}}()
error = GRBnewmodel(env, model_p, "qp", 0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
model = model_p[]

# variables bounds
vlb = Cdouble[-Inf, -Inf, -Inf]
vub = Cdouble[+Inf, +Inf, +Inf]

# linear objective coeffs
cobj = Cdouble[2.0, 0.0, 0.0]

# quadratic objective coeffs
qrow = Cint[0, 0, 1, 1, 2]
qcol = Cint[0, 1, 1, 2, 2]
qval = Cdouble[1, 1, 1, 1, 1]

# linear constraint1 coeffs
c1ind = Cint[0, 1, 2]
c1val = Cdouble[1.0, 2.0, 3.0]
c1sense = GRB_GREATER_EQUAL
c1rhs = 4.0

# linear constraint2 coeffs
c2ind = Cint[0, 1]
c2val = Cdouble[1.0, 1.0]
c2sense = GRB_GREATER_EQUAL
c2rhs = 1.0

# linear constraint3 coeffs: nonnegativity
c3xind = Cint[0]
c3xval = Cdouble[1.0]
c3xsense = GRB_GREATER_EQUAL
c3xrhs = 0.0
c3yind = Cint[1]
c3yval = Cdouble[1.0]
c3ysense = GRB_GREATER_EQUAL
c3yrhs = 0.0
c3zind = Cint[2]
c3zval = Cdouble[1.0]
c3zsense = GRB_GREATER_EQUAL
c3zrhs = 0.0
c3ind = [c3xind, c3yind, c3zind]
c3val = [c3xval, c3yval, c3zval]
c3sense = [c3xsense, c3ysense, c3zsense]
c3rhs = [c3xrhs, c3yrhs, c3zrhs]

# variables and objective
vname = Cstring["xyz"]
error = GRBaddvars(
  model, # model
  3,      # : numvars
  0,      # : numnz
  C_NULL, # : *vbeg
  C_NULL, # : *vind
  C_NULL, # : *vval
  cobj,   # : *obj
  vlb,    # : *lb
  vub,    # : *ub
  C_NULL, # : *vtype
  # "xyz" # : **varnames
  ["xyz"]   # : **varnames
  # Cuchar["xyz"]   # : **varnames
  # Base.unsafe_convert(Cstring,"xyz")
  # Base.cconvert(Cstring,"xyz")
)
error = GRBaddqpterms(model, 5, qrow, qcol, qval)

# constraint1
error = GRBaddconstr(
  model,   # : *model
  3,       # : numnz
  c1ind,   # : *cind
  c1val,   # : *cval
  c1sense, # : sense
  c1rhs,   # : rhs
  "c1",    # : *constrname
)

# constraint2
error = GRBaddconstr(
  model,   # : *model
  2,       # : numnz
  c2ind,   # : *cind
  c2val,   # : *cval
  c2sense, # : sense
  c2rhs,   # : rhs
  "c2",    # : *constrname
)

# constraint3: nonnegativity
for i in 1:3
  error = GRBaddconstr(
    model,   # : *model
    1,       # : numnz
    c3ind[i],   # : *cind
    c3val[i],   # : *cval
    c3sense[i], # : sense
    c3rhs[i],   # : rhs
    "c3_$i",    # : *constrname
  )
end
error = GRBoptimize(model)
error = GRBwrite(model, "qp.lp");

pinfeas = Ref{Cdouble}()
dinfeas = Ref{Cdouble}()
relgap = Ref{Cdouble}()
NumVars = Ref{Cint}()
NumConstrs = Ref{Cint}()
IterCount = Ref{Cint}() # simplex iters
BarIterCount = Ref{Cint}() # barrier iters
optimstatus = Ref{Cint}() # barrier iters
objval = Ref{Cdouble}() # barrier iters
sol = ones(3)

GRBgetdblattr(model, "ConstrVio", pinfeas) # maximum (primal) constraint violation
GRBgetdblattr(model, "MaxVio", dinfeas) # sum of (dual) constraint violations
GRBgetdblattr(model, "ComplVio", relgap) # complementarity violation
GRBgetintattr(model, "NumVars", NumVars) # sum of (dual) constraint violations
GRBgetintattr(model, "NumConstrs", NumConstrs) # sum of (dual) constraint violations
GRBgetintattr(model, "BarIterCount", BarIterCount) # sum of (dual) constraint violations

error = GRBgetintattr(model, GRB_INT_ATTR_STATUS, optimstatus);
error = GRBgetdblattr(model, GRB_DBL_ATTR_OBJVAL, objval);
error = GRBgetdblattrarray(model, GRB_DBL_ATTR_X, 0, 3, sol);

Now trying to access the variable values by name "xyz":

julia> GRBgetvarbyname(model, "xyz", Cint[0,1,2])
0

julia> GRBgetvarbyname(model, ["xyz"], Cint[0,1,2])
0

julia> xyz = Cdouble[1.23, 1.23, 1.23]
3-element Vector{Float64}:
 1.23
 1.23
 1.23

julia> GRBgetvarbyname(model, ["xyz"], xyz)
0

julia> xyz
3-element Vector{Float64}:
 1.2300004959106443
 1.23
 1.23

julia> GRBgetvarbyname(model, "xyz", xyz)
0

julia> xyz
3-element Vector{Float64}:
 1.2299995422363281
 1.23
 1.23

See the documentation for usage: https://www.gurobi.com/documentation/9.5/refman/c_getvarbyname.html

varnumP = Ref{Cint}(-1)
ret = GRBgetvarbyname(model, "xyz", varnumP)
@assert ret == 0
varnumP[]

Thanks, but I am having trouble understanding what GRBgetvarbyname returns. Based on the documentation, I thought it would return the values in the array called "xyz", but this call doesn’t seem to return anything. Also, I was unable to name my variable "xyz" and instead had to call it ["xyz"], e.g., in:

error = GRBaddvars(
  model, # model
  3,      # : numvars
  0,      # : numnz
  C_NULL, # : *vbeg
  C_NULL, # : *vind
  C_NULL, # : *vval
  cobj,   # : *obj
  vlb,    # : *lb
  vub,    # : *ub
  C_NULL, # : *vtype
  ["xyz"]   # : **varnames #! worked
  # "xyz" # : **varnames #! didn't work
  # Cuchar["xyz"]   # : **varnames  #! didn't work
  # Base.unsafe_convert(Cstring,"xyz")  #! didn't work
  # Base.cconvert(Cstring,"xyz")  #! didn't work
)

Am I interpreting GRBgetvarbyname correctly? The docs only say:

Retrieves a variable from its name. If multiple variables have the same name, this routine chooses one arbitrarily.

I am having trouble understanding what GRBgetvarbyname returns

It returns an error code, which is 0 if no error occurred. varnumP[] is the 0-based Cint index of the variable in Gurobi.

You’re also adding three variables. Did you mean ["x", "y", "z"]?

What are you trying to achieve with GRBgetvarbyname?

I am trying to obtain the value of variables in the solved model, something like equivalent of value.(x) in JuMP. I see that I need to also use GRBgetvars and am trying to figure that out.

Ah yes, I think I did mean ["x","y","z"] for the above code, but in general can you add/access vectors of variables? What I am trying to do is provide a name vector with a single element "xyz" that corresponds to the vector of variables [x,y,z]. I tried the following unsuccessfully (though I also cannot get it to work by accessing the variables individually with names ["x","y","z"] either):

varnumPxyz = Ref{Cint}(-1)
ret = GRBgetvarbyname(model, "xyz", varnumPxyz)
@assert ret == 0
varnumPxyz[] # = 0
xyz = Cdouble[0.1, 0.1, 0.1] # want to overwrite this
GRBgetvars(
  model,
  3, # : *numnzP
  Cint[varnumPxyz[]], # : *vbeg
  Cint[varnumPxyz[],varnumPxyz[]+1,varnumPxyz[]+2], # : *vind
  xyz, # : *vval
  0, # : start
  3, # : len
)

but in general can you add/access vectors of variables?

Yes. That’s what you did with GRBaddvars?

GRBgetvars

GRBgetvars is almost certainly not what you want to do, nor is GRBgetvarbyname.

I am trying to obtain the value of variables in the solved model

You need to query the X variable attribute:

See:

solution = zeros(Cdouble, num_variables)
ret = GRBgetdblattrarray(model, "X", 0, num_variables, solution)
@assert ret == 0

As you’re probably starting to realize, the C API is very bare-bones. It provides no helper functionality to make things easy.

If you need help for which functions to call, I’d look through the Gurobi.jl source code. For example, here is where we query the primal solution value of a single decision variable:

The key points are:

  • Every variable and constraint is a 0-based Cint index in Gurobi. You’re responsible for remembering which variables and constraints correspond to which indices. In general, you don’t need to provide names, and using Gurobi to lookup string names is usually the wrong thing to do.
  • Building the model is mostly function-based
  • Querying the model is mostly via attributes with GRBgetdblattr...

You could also follow the Gurobi C examples:

You should be able to see how they map onto the Julia functions.

Ah ok I see how that is one way that can work. It will get all the variables in the model, then I just have to remember the indexing and split it up.

Why is that? Am I wrong in thinking that there is an alternative using GRBgetvarbyname and then GRBgetvars? Namely:

  1. use GRBgetvarbyname to get the start index of a named variable
  2. pass this start index + length + output array to GRBgetvars
  3. GRBgetvars fills the output array (this array is what I try to pass in the argument location corresponding to vval in the Gurobi docs)

Ah ok I see how that is one way that can work

There is also GRBgetdblattrelement and GRBgetdblattrlist.

Why is that?

It’s easy to introduce a bug if two variables have the same name. (Gurobi picks randomly, so it may work initially, and then fail with other changes.)

Instead, create your own data structure in Julia that maps variables to their Gurobi indices.

If you’re going to be adding more complexity and data structures, then just use Gurobi.Optimizer() and the MOI API.

If you really want to though, you could use something like:

function value(model, name::String)
    colP = Ref{Cint}(-1)
    ret = GRBgetvarbyname(model, name, colP)
    @assert ret == 0
    valP = Ref{Cdouble}(NaN)
    ret = GRBgetdblattrelement(model, "X", colP[], valP)
    @assert ret == 0
    return valP[]
end

GRBgetvars

GRBgetvars doesn’t return the solution. It’s used to query bounds and constraint coefficients. See the documentation: https://www.gurobi.com/documentation/10.0/refman/c_getvars.html

1 Like

Oh I see, I kept looking through the docs but I didn’t understand. Thanks for your patience!!

1 Like

This is the main reason we don’t actively advertise the Gurobi C API. It’s powerful, but it’s not very friendly to users :laughing:

1 Like