Gurobi Cookbook - Column-wise Formation (5/5)
written by Abremod, LLC on 2019-07-01
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)
written by Abremod, LLC on 2019-06-30
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)
written by Abremod, LLC on 2019-06-29
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)
written by Abremod, LLC on 2019-06-28
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)
written by Abremod, LLC on 2019-06-27
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
cutofftarget. - 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)