Gurobi Cookbook - Column-wise Formation (5/5)

Using a column-wise formulation

Often, is it more convenient to think of the contraints that a variable participates in instead of the variables that are part of a constraint. Also, in a column-generation procedure, this is required. To support this, Gurobi has a Column object that can be used to when defining a model. Let's create the same model as before.

model = grb.Model()
c0 = model.addConstr(0, '>', 10, name='c0')
c1 = model.addConstr(0, '<', 4, name='c1')
model.update()
x = model.addVar(name='x', obj=2, column=grb.Column([1,  1], [c0, c1]))
y = model.addVar(name='y', obj=3, column=grb.Column([1, -1], [c0, c1]))
model.update()
model.optimize()
Optimize a model with 2 rows, 2 columns and 4 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+00]
  Objective range [2e+00, 3e+00]
  Bounds range    [0e+00, 0e+00]
  RHS range       [4e+00, 1e+01]
Presolve time: 0.00s
Presolved: 2 rows, 2 columns, 4 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   1.000000e+01   0.000000e+00      0s
       2    2.3000000e+01   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.01 seconds
Optimal objective  2.300000000e+01

commentary

The column-wise approach often looks awkward. Linear programming formualtions are almost always documented in a row-wise format and it is natural to think in terms of the variables being added to constraints. Because of that, this method of creating models is often restricted to column-generation proceedures, where they are required and are probably under-utilized in other cases.

Obtaining the dual solution

When computing the solution to a linear problem (with no integer constraints), gurobi will also compute dual prices for the contraints. This has the economic interpretation as the change in the value of the objective as the right hand side of the contraint is increased or decreased, and it is used in Column Generation procedures. The attribute is Pi

constr_names = model.getAttr('ConstrName', model.getConstrs())
constr_values = model.getAttr("Pi", model.getConstrs())
zip(constr_names, constr_values)
[('c0', 2.5), ('c1', -0.5)]

commentary

In this example, the dual price of c0 ($x + y \ge 10$) was 2.5. That means, increasing the right hand side from 10, to say 10.1 and resolving would increase the objective from 2.3 to at least 2.55, and reducing the right hand side from 10 to 9.9 would reduce the objective by at most .25 to at least 2.05. We say at least 2.55 because we are minimizing.



Gurobi Cookbook - Linear Model (4/5)

Creating a linear programming model

We will formulate the following linear program

minimize 2 x + 3 y
subject to:
x + y >= 10
x - y <= 4
x, y >= 0
model = grb.Model()
x = model.addVar(name='x', obj=2)
y = model.addVar(name='y', obj=3)
model.update()
c0 = model.addConstr(x + y >= 10, name='c0')
c1 = model.addConstr(x - y <= 4, name='c1')
model.update()
model.optimize()
if model.solCount > 0:
    print("objective = ", model.objVal)
else:
    print("solution status = ", model.Status)
Optimize a model with 2 rows, 2 columns and 4 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+00]
  Objective range [2e+00, 3e+00]
  Bounds range    [0e+00, 0e+00]
  RHS range       [4e+00, 1e+01]
Presolve time: 0.00s
Presolved: 2 rows, 2 columns, 4 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   1.000000e+01   0.000000e+00      0s
       2    2.3000000e+01   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.01 seconds
Optimal objective  2.300000000e+01
objective =  23.0

commentary

We first create a fresh model object with the Model constructor. Then we create the two variables, x and y. By default, Gurobi variables are restrictred to be nonnegative. Since the variables are nonnegative in our mathematical model, we don't need to set any bounds. We do set the objective function. Before adding constraints that reference the variables, we make a call to update. Alternatively, we could have set the parameter UpdateMode to 1.


Gurobi Cookbook - Relaxed Function (3/5)

Obtaining the lp relaxation of a MIP

The function relaxed of the model object creates a copy of the problem with the intgrality constraints relaxed.

relaxed = m.relax()
relaxed.optimize()
Optimize a model with 637 rows, 120 columns and 14280 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+00]
  Objective range [1e+00, 1e+00]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 1e+01]
Presolve time: 0.01s
Presolved: 120 rows, 757 columns, 14400 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0   -0.0000000e+00   0.000000e+00   1.752000e+03      0s
     414    1.7142857e+01   0.000000e+00   0.000000e+00      0s

Solved in 414 iterations and 0.04 seconds
Optimal objective  1.714285714e+01

commentary

The relax function returns a fresh of the original model. The original model will still have the original constraints. You can see this by querying the attribute isMIP.

cplex users

In the cplex interactive shell, this corresponds to the command change problem lp except that this command will modify the global model instead of returning a new model

m.isMIP, relaxed.isMIP
(1, 0)

commentary

The dual price is available only for continous models.

Attempting to query this attribute for MIPs (or an continuous model without a solution) will result in the exception

GurobiError: Unable to retrieve attribute 'Pi'

If you are unsure that the model has integer variables, you can check with m.isMIP


Gurobi Cookbook - Bulk Values (2/5)

Getting values in bulk for a large number of variables or constraints

See previous blog post for the basics of Gurobi.

var_names = m.getAttr('VarName', m.getVars())
constr_names = m.getAttr('ConstrName', m.getConstrs())

commentary

For users of the compute server, this method is almost mantatory. The alternative [v.varName for v in m.getVars()], for example, will be orders of magnitude slower.

Getting nonzero values of the solution

In many models, only a small portion of the variables have nonzero values. In that case, it is usually more convenient to get a list of only the variables that have nonzero values.

epsilon = 1e-6
var_names = m.getAttr('VarName', dvars)
var_values = m.getAttr('X', dvars)
[(name, value) for name, value in zip(var_names, var_values)
               if abs(value) >= epsilon]
[('x0', 1.0),
 ('x7', 1.0),
 ('x18', 1.0),
 ('x27', 1.0),
 ('x33', 1.0),
 ('x40', 1.0),
 ('x46', 1.0),
 ('x50', 1.0),
 ('x60', 1.0),
 ('x66', 1.0),
 ('x73', 1.0),
 ('x81', 1.0),
 ('x88', 1.0),
 ('x94', 1.0),
 ('x101', 1.0),
 ('x109', 1.0),
 ('x112', 1.0),
 ('x116', 1.0),
 ('x117', 1.0),
 ('x119', 1.0)]

Commentary

This solution uses the builk query mechanism m.getAttr to get the values, then test each value for nonzeroness. It is not sufficient to test for 0 exactly, as Gurobi will treat very small numbers as zero. To avoid this, we check that the absolute value of a solution is greater than some small epsilon.

Getting the contributions to the objective function

An interesting subset of variables are the ones that wind up affecting the objective function directly. To contribute to an objective, a variable must be nonzero and have a nonzero coeficient in the objective function.

obj_coeffs = m.getAttr('Obj', dvars)
[(name, obj_coeff * value) 
 for name, value, obj_coeff in zip(var_names, var_values, obj_coeffs)
 if abs(obj_coeff * value) > epsilon]
[('x0', 1.0),
 ('x7', 1.0),
 ('x18', 1.0),
 ('x27', 1.0),
 ('x33', 1.0),
 ('x40', 1.0),
 ('x46', 1.0),
 ('x50', 1.0),
 ('x60', 1.0),
 ('x66', 1.0),
 ('x73', 1.0),
 ('x81', 1.0),
 ('x88', 1.0),
 ('x94', 1.0),
 ('x101', 1.0),
 ('x109', 1.0),
 ('x112', 1.0),
 ('x116', 1.0),
 ('x117', 1.0),
 ('x119', 1.0)]

Gurobi Cookbook - Basics (1/5)

This is a collection of short, idiomatic python code for working with the Gurobi solver using its python interface.

from __future__ import print_function
# get a sample mps file from the benchmark repository
import requests
import os
import toolz
benchmark_url = 'http://miplib.zib.de/download/'
sample_mps = 'cov1075.mps.gz' # "timtab1.mps.gz"
infile = requests.get(benchmark_url + sample_mps)
if not os.path.exists(sample_mps):
    with open(sample_mps, 'w') as outfile:
       outfile.write(infile.content)

Setting up Gurobi

If you will be using the python interface, the best way to install Gurobi and keep it up to date is to use the Anaconda python distribution. After this in installed, you need to import the gurobi library. We usually use the aliases grb for gurobipy, and GRB for the namespace where gurobi keeps its constants.

import gurobipy as grb
GRB = grb.GRB

Checking the gurobi version

under the gurobi library, there is a gurobi class which contains some utilities that aren't part of typical modeling modeling. The gurobi object has a version function which returns a tuple containing the major, minor, and patch level of the gurobi installation.

grb.gurobi.version()
(6, 5, 2)

Reading a saved model

The gurobi python module has a read function that will create an in-memory model from a model saved in an mps, lp, or several other formats. It can also read directly from several standard compressed formats, like, zip, gzip, and bzip.

m = grb.read(sample_mps)

commentary

In this case the variable m is a reference to the in-memory model equivalent to what was stored in the file. Gurobi uses the extension to determine the format of the file. mps files the best way to save a representation of the model. They can be used to compare performance of different solvers, or a single solver with different parameters. In this notebook, we are using a benchmark problem. While the mps file can store an exact copy of the model, it is difficult for a human to read. The lp format is more readable version of the model but it is lossy, so you may get different results reading from an lp file, especially if you have large or small coefficients in your model.

for cplex users

If you are used to the interactive cplex, this is the rough equivalent if read problem. The cplex interactive shell has an implicit model in its global namespace, and all actions operate on the global model. The model

Solving a model

A gurobi model object has an optimize function that will start the solve process.

m.optimize()
Optimize a model with 637 rows, 120 columns and 14280 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+00]
  Objective range [1e+00, 1e+00]
  Bounds range    [1e+00, 1e+00]
  RHS range       [1e+00, 1e+01]
Found heuristic solution: objective 56
Presolve time: 0.01s
Presolved: 637 rows, 120 columns, 14280 nonzeros
Variable types: 0 continuous, 120 integer (120 binary)
Presolved: 120 rows, 757 columns, 14400 nonzeros


Root relaxation: objective 1.714286e+01, 414 iterations, 0.04 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   17.14286    0  120   56.00000   17.14286  69.4%     -    0s
H    0     0                      35.0000000   17.14286  51.0%     -    0s
H    0     0                      27.0000000   17.14286  36.5%     -    0s
H    0     0                      25.0000000   17.14286  31.4%     -    0s
     0     0   17.25000    0  109   25.00000   17.25000  31.0%     -    0s
     0     0   17.28125    0  112   25.00000   17.28125  30.9%     -    0s
     0     0   17.28125    0  111   25.00000   17.28125  30.9%     -    0s
H    0     0                      24.0000000   17.28125  28.0%     -    0s
     0     0   17.29164    0  112   24.00000   17.29164  28.0%     -    0s
H    0     0                      23.0000000   17.29164  24.8%     -    0s
     0     0   17.29941    0  111   23.00000   17.29941  24.8%     -    0s
     0     0   17.30239    0  109   23.00000   17.30239  24.8%     -    0s
     0     0   17.30241    0  111   23.00000   17.30241  24.8%     -    0s
     0     0   17.30959    0  109   23.00000   17.30959  24.7%     -    0s
     0     0   17.31158    0  112   23.00000   17.31158  24.7%     -    0s
     0     0   17.31295    0  108   23.00000   17.31295  24.7%     -    0s
H    0     0                      21.0000000   17.31295  17.6%     -    0s
     0     0   17.31380    0  110   21.00000   17.31380  17.6%     -    0s
     0     0   17.31524    0  111   21.00000   17.31524  17.5%     -    0s
     0     0   17.31562    0  110   21.00000   17.31562  17.5%     -    0s
     0     0   17.31630    0  111   21.00000   17.31630  17.5%     -    0s
     0     0   17.31655    0  114   21.00000   17.31655  17.5%     -    0s
     0     0   17.31739    0  113   21.00000   17.31739  17.5%     -    0s
     0     0   17.31752    0  114   21.00000   17.31752  17.5%     -    0s
     0     0   17.31772    0  113   21.00000   17.31772  17.5%     -    0s
H    0     0                      20.0000000   17.31772  13.4%     -    0s
     0     0   17.31792    0  113   20.00000   17.31792  13.4%     -    0s
     0     0   17.31878    0  113   20.00000   17.31878  13.4%     -    0s
     0     0   17.31933    0  112   20.00000   17.31933  13.4%     -    1s
     0     0   17.31990    0  114   20.00000   17.31990  13.4%     -    1s
     0     0   17.32008    0  112   20.00000   17.32008  13.4%     -    1s
     0     0   17.32011    0  112   20.00000   17.32011  13.4%     -    1s
     0     0   17.32057    0  112   20.00000   17.32057  13.4%     -    1s
     0     0   17.32068    0  113   20.00000   17.32068  13.4%     -    1s
     0     0   17.32111    0  112   20.00000   17.32111  13.4%     -    1s
     0     2   17.32111    0  112   20.00000   17.32111  13.4%     -    1s

Cutting planes:
  Gomory: 5
  Zero half: 16

Explored 301 nodes (20763 simplex iterations) in 3.57 seconds
Thread count was 8 (of 8 available processors)

Optimal solution found (tolerance 1.00e-04)
Best objective 2.000000000000e+01, best bound 2.000000000000e+01, gap 0.0%

commentary

Gurobi will work solving the model until either

  • It finds a provably optimal solution (within tolerance)
  • It proves that the model in infeasible or unbounded
  • It proves no solution exists with an objective better than the cutoff target.
  • It hits a limit
    • Wall clock time
    • Number of solutions found (MIPs only)
    • Number of search nodes explored (MIPs only)
    • Number of Simplex Iterations (rarely used)
  • It receives an interrupt signal
  • It encounters an error
cplex users

This is equivalent to solve in interactive cplex

Checking if Gurobi found a solution

From the log, we can see that an optimal solution was found. In application code, it's important to explicitly check the model status after each solve. The cleaneast way to determine if Gurobi found a feasible solution is to query the SolCount attribute. If the solution count is at least 1, then Gurobi will have a feasible solution available. If not, querying for the values in the solution will through an exception.

print(m.SolCount, m.Status)
8 2

The Status attribute has good information as well, but it is not obvious from the status if there is a feasible solution available. In this case, an optimal solution. A status of optimal or suboptimal imply the existence of a solution and a status of loaded, infeasible, inf_or_unbd, unbounded, cutoff, or numeric definitely mean that no feasible solution was found. However, the statuses iteration_limit, node_limit, time_limit, solution_limit or interrupted refer to why the solve terminated. For these cases, all you know is that the solution, if any, isn't optimal according to the optimality criteria. That's why checking the solution count is recommended.

Its always a good idea to explitly check for feasibility in any production application, even if you reliably get feasible solutions from Gurobi. Unless you have complete control over the data you receive, inevitibly your application will receive data that results in an infeasible model. If you don't check for fesibility, your application will throw the unintuive exception.

GurobiError: Unable to retrieve attribute 'X'

Getting the objective value of a solution

If a solution does exist, you will usually be interested in the objective value. Most often, the objective value, as this is usually tied to some business metric. This can be done with the ObjVal attribute. Gurobi also computes a bound on the best possible objective in a feasible solution, which can be queried with the ObjBound attribute.

m.ObjVal, m.ObjBound
(20.0, 20.0)

Getting solution quality information

Gurobi works with floating point arithmetic, which has inevitible round-off errors. Even if you have a pure integer problem with all integer coeficients, Gurobi still uses floating point numbers to solve the LP subproblems. Gurobi can produce solutions that mathematically violate constraints in what are usually very small amounts. For each variable bound, and each constraint, there is potentially a small violation. By default, Gurobi allows variables and constraints to be violated by 0.000001, but for most problems, the violation is much smaller.

m.BoundVioSum, m.ConstrVioSum, m.BoundVio, m.ConstrVio
(0.0, 0.0, 0.0, 0.0)

Getting solution values

Of course, you will want to get the actual values of the solution. In any non-trivial model, there will be a set of decision variables. The complete list of deicision variables can be obtained using the function getVars. If the variables have all been given unique names, you can get an individual variable with the getVarByName function. Values can be obtained by querying the X attribute of a variable. The name of the variable.

dvars = m.getVars()
dvars[0].VarName, dvars[0].X
('x0', 1.0)