Hidden Gems

Useful Gurobi features you may not know

Elisabeth Rodríguez Heck

Gurobi Live, Barcelona, 2023


  • You probably have used Gurobi to solve LPs, MIPs or QPs
  • We’ll use our most popular API through the talk: gurobipy
import gurobipy as gp
from gurobipy import GRB

with gp.Env() as env, gp.Model(env=env) as model:
    model = gp.Model("small")

    x = model.addVar(vtype=GRB.BINARY, name="x")
    y = model.addVar(vtype=GRB.BINARY, name="y")
    z = model.addVar(vtype=GRB.BINARY, name="z")

    model.setObjective(x + y + 2 * z, GRB.MAXIMIZE)

    # Linear constraints
    model.addConstr(x + 2 * y + 3 * z <= 4, "c0")
    model.addConstr(x + y >= 1, "c1")
    # Quadratic constraints
    model.addConstr(x*x + y*y <= 4.0, "qc0")

    model.optimize()

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

General Constraints




  • There are different strategies for the algorithmic implementation…


  • … but the API is consistent
model.addGenConstrXXX()

# XXX = Max, Min, Exp, ...

Gurobi API for Function Constraints - Example


Example: Natural Exponential

\[ \begin{align*} \min ~& -2x + e^x \\ \text{s. t. } ~& 0 \leq x \leq 1 \end{align*} \]

model = gp.Model("gen")

x = model.addVar(lb=0, ub=1, name="x")
y = model.addVar(name="y")

# y = exp(x)
gc = model.addGenConstrExp(x, y)

model.setObjective(-2 * x + y)

model.optimize()

Simple Constraints vs. Function Constraints


Simple Constraints

  • Modeling convenience feature
  • Represented exactly in terms of binary variables, linear and SOS constraints
  • Takes the burden of knowing how to model these off the user
  • Good for code readability and maintenance

Function Constraints

  • Much more complex relationships between variables: \(y = f(x)\) where \(x\) and \(y\) are Gurobi variables, \(f(.)\) can be selected out of:
    • Polynomial, Natural exponential, Exponential, Natural Logarithm, Logarithm, Logistic, Power, Sine, Cosine, Tangent
  • No simple way to re-state these constraints using integer variables and linear constraints
  • Alternative algorithms to standard MIP are needed
  • Models much more difficult to solve
  • Only univariate functions supported

Function Constraints - Current State


Static Piecewise-Linear (PWL) approximation: nonlinear function is approximated once before MIP solving starts


Cost vs. Accuracy

  • More accurate PWL → more segments
  • More segments → problem harder to solve

Accuracy controlled with constraint attributes or globally with model parameters

Underestimate, overestimate or somewhere in between

Function Constraints - Starting with Gurobi 11


  • Dynamic PWL approximation using Outer Approximation and Spatial Branch-and-Bound:
    • At each node of the Branch-and-Bound search the solver computes an approximation in the neighborhood of the current relaxation solution
  • API details: set globally as a parameter or for specific constraints as an attribute:
    • FuncNonlinear = 1 (dynamic outer approximation)
    • FuncNonlinear = 0 (default, static PWL approximation)
  • Usually produces more accurate solutions (with less violations in the original nonlinear constraint), but it can be expensive

What is better, static or dynamic?

  • Many pieces with static PWL may be even harder to solve than using Outer Approximation
  • The most effective approach will depend on your model and accuracy goals
  • Test before deciding!

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Multiple Objectives


Real-world optimization problems often have multiple, competing objectives


Maximize Profit

&

Minimize Late Orders

Minimize Shift Count

&

Maximize Worker Satisfaction

Minimize Cost

&

Maximize Product Durability

Maximize Profit

&

Minimize Risk

How does Gurobi handle the trade-offs?


  • Weighted or Blended: Optimize a weighted combination of the individual objectives
flowchart LR
    id1[Obj 1] --- id2[Obj 2] --- id3[Obj 3]

\[ \begin{align*} \min & ~w_1 f_1(x) + w_2 f_2(x) + w_3 f_3(x) \\ \text{s. t. } & ~x \in C \end{align*} \]

  • Hierarchical or Lexicographical: Optimize each objective in a given priority order while limiting the degradation of the higher-priority objectives
flowchart TD
    id1[Obj 1] --> id2[Obj 2] --> id3[Obj 3]

\[ \begin{align*} \min & ~f_1(x) \\ \text{s. t. } & ~x \in C \end{align*} \]

\[ \begin{align*} \min & ~f_2(x) \\ \text{s. t. } & ~x \in C \\ & ~f_1(x) \leq \epsilon_1 \end{align*} \]

\[ \begin{align*} \min & ~f_3(x) \\ \text{s. t. } & ~x \in C \\ & ~f_1(x) \leq \epsilon_1 \\ & ~f_2(x) \leq \epsilon_2 \end{align*} \]

  • Weighted + Hierarchical

Multiple Objectives API


model.setObjectiveN(expr, index, priority=0, weight=1,
                    abstol=1e-6, reltol=0, name="")

# expr (LinExpr): New alternative objective
# index (int): Index for new objective (to set parameters or
#              query solution per objective)
# priority (int, optional): Objective's priority (ObjNPriority attribute)
# weight (float, optional): Objective's weight (ObjNWeight attribute)
# abstol (float, optional): Absolute tolerance used in calculating
#                           allowable degradation (ObjNAbsTol attribute)
# reltol (float, optional): Relative tolerance used in calculating
#                           allowable degradation (ObjNAbsTol attribute)
# name (string, optional): Objective's name

Multiple Objectives API - Example


Workforce Scheduling:

  • First, minimize cost
  • Second, make the model “fair”


# Set global sense for ALL objectives
model.ModelSense = GRB.MINIMIZE

# Set up primary objective: minimize total pay cost
model.setObjectiveN(gp.quicksum(pay[w] * x[w, s] for w, s in availability),
                index=0, priority=2, abstol=20,
                reltol=0, name='Cost')

# Set up secondary objective: minimize difference in shift length
model.setObjectiveN(maxShift - minShift, index=1, priority=1, name='Fairness')

How is the degradation value calculated?


flowchart TD
    id1[Obj 1] --> id2[Obj 2]

Additional constraint:

Obj 1 \(\leq\) rhs

\[ \begin{align*} \min & ~f_2(x) \\ \text{s. t. } & ~x \in C \\ & ~f_1(x) \leq \epsilon_1 \end{align*} \]

base_value = max(objbnd +|objval|*MIPGap, objbnd + MIPGapAbs, objval)
relaxed    = max(ObjNRelTol*|base_value|, ObjNAbsTol)
rhs (𝝐_𝟏)    = base_value + relaxed

objbnd    : best bound of objective OBJ 1
objval    : best solution value for objective OBJ 1
MIPGap    : relative MIP gap
MIPGapAbs : absolute MIP gap
ObjNRelTol: further allowable relative degradation for OBJ 1
ObjNAbsTol: further allowable absolute degradation for OBJ 1

Why use the Multiple Objectives API?


  • Make the objective functions easy to understand and maintain
  • Get faster performance with warm starts for hierarchical objectives
  • Avoid numerical issues with large objective coefficients
    • “There are 45 coefficients in 2 distinct groups. Is this a multi-objective case in hiding?”

Multiple Objectives API - Details


  • Single objective sense for all objectives controlled via the ModelSense attribute
    • Can always change the sense of the objective multiplying by -1
  • Objective expressions must be linear
  • Parameters can be set on each objective using multi-objective environments
  • Callbacks are available

Too many objectives?

  • If your model has too many objectives (e.g. more than 10) consider whether you really need them
    • Hierarchical: too many objectives can result in a too small search space
    • Weighted: too many objectives need many groups of distinct, widespread weights, which can result in numerical issues
  • Possible alternative: use precedence constraints to model customer priorities

Multiple Objectives - Logging


  • Weighted: log will be the same as for a single-objective model

  • Hierarchical:

--------------------------------------------------------------------
Multi-objectives: starting optimization with 3 objectives ...
--------------------------------------------------------------------
[...]
---------------------------------------------------------------------------
Multi-objectives: optimize objective 1 (Name) ...
---------------------------------------------------------------------------
  • Mixed:
---------------------------------------------------------------------------
Multi-objectives: starting optimization with 5 objectives (3 combined) ...
---------------------------------------------------------------------------
[...]
---------------------------------------------------------------------------
Multi-objectives: optimize objective 1 (weighted) ...
---------------------------------------------------------------------------

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Scenario Analysis: Modeling What-If?


Idea: One model → solutions to multiple scenarios

Examples:

  • “What happens to the optimized production if our forecast is 10% lower than anticipated?”
  • “What happens to the building costs if phase 1 is 1 day late? 2 days late? 5 days late?”

Motivation to use the Gurobi Multi-Scenario API:

  • Solves faster than separate models
  • Easier to maintain and understand

Changes you can make:

  • Linear objective function coefficients
  • Variable lower and upper bounds
  • Linear constraint right-hand-side values

Multiple Scenarios API


Multiple Scenarios - Tips & Tricks

Multiple scenarios API restrictions

It is not possible to explicitly:

  • Add/remove variables or constraints
  • Change variable types
  • Change the sense of constraints

Tricks to circumvent some restrictions

  • To remove a variable, set its bounds to zero
  • To add a variable to a scenario, add it to the base model with bounds set to zero and then change the bounds accordingly
  • To remove a constraint, change its RHS values to GRB.INFINITY / -GRB.INFINITY
  • To add a constraint to a scenario or change its sense, add it as a pair of inequalities to the base model and change its RHS values accordingly

Multiple Scenarios - Logging

Phase 1: Find the best solution over all scenarios

[...]
Solving a multi-scenario model with 2 scenarios...
[...]
    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

H    0     0                    1.850017e+09 8.0000e+08  56.8%     -    0s

Optimal solution found at node 131289 - now completing multiple scenarios...

Phase 2: Find the best solution for each scenario

  • Worst Incumbent: The worst solution found over all scenarios
  • BestBd: Best possible objective value for any solution that has not yet been found
    Nodes    |    Current Node    |   Scenario Obj. Bounds    |     Work
             |                    |   Worst                   |
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

 1312887 143640 1.5000e+09   58   14 1.5750e+09 1.5000e+09  4.76%   8.5  218s
[...]
Best objective 1.500012833333e+09, best bound 1.500012833333e+09, gap 0.0000%

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Multiple Solutions


  • Idea: present alternative solutions
    • The model may lack implicit elements like preferences, or some aspects of the objective may be difficult to quantify
    • Demonstrate value by comparing alternatives to the optimal solution
    • Gives a greater feeling of control
    • Get feedback, alternate solution may show missing modeling elements

  • How can you quickly report several feasible solutions?
    • Re-run the model with a manually added constraint cutting off the previously reported solution
    • Define a Solution Pool and report multiple solutions automatically after a single run

Subtleties and Limitations

  • Continuous variables: multiple equivalent solutions will not be reported as per our definitions.

Solution Pool - Setup

Parameter settings Behavior
PoolSearchMode=0 - Stores all solutions found in the regular optimization
- No additional tree search performed
PoolSearchMode=1
PoolSolutions=n		 - Stores n-1 additional solutions to the optimal solution
- Controls how many solutions to save
PoolSearchMode=2
PoolSolutions=n
PoolGap=x		 - Stores n-1 best solutions with a MIPGap less than x% in addition to the optimal solution
- Requires further tree search than PoolSearchMode=1		 


# Limit how many solutions to collect
model.setParam(GRB.Param.PoolSolutions, 100)

# Limit the search space by setting a gap for the
# worst possible solution that will be accepted
model.setParam(GRB.Param.PoolGap, 0.10)

# Do a systematic search for the k best solutions
model.setParam(GRB.Param.PoolSearchMode, 2)

Multiple Solutions - Logging


Phase 1: Find one provably optimal solution (identical to MIP log except for last line)

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

     0     0   65.87500    0    1   65.00000   65.87500  1.35%     -    0s

Optimal solution found at node 0 - now completing solution pool...


Phase 2: Populate solution pool

  • Worst Incumbent: objective value of the worst solution among the ones in the pool
    Nodes    |    Current Node    |      Pool Obj. Bounds     |     Work
             |                    |   Worst                   |
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0   65.00000    0    1          -   65.00000      -     -    0s

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Python Matrix-Friendly API


Idea: build optimization models directly in terms of matrices and vectors

\[ \begin{align*} \min &~x^T Q x \\ \text{s. t. } &~Ax \geq b \\ &~x \geq 0 \end{align*} \]

Matrix variables: MVar objects:

  • An object representing a matrix (or vector) of optimization variables
model.addMVar(shape, lb=0.0, ub=float('inf'), obj=0.0,
              vtype=GRB.CONTINUOUS, name="")

# shape: tuple with dimensions of the variable (like Numpy's ndarray)

# Example: Add a 4-by-2 matrix binary variable
x = model.addMVar((4,2), vtype=GRB.BINARY)

Python Matrix-Friendly API - How To


Building linear constraints with MVars:

  • use python3’s matrix multiplication operator @ to build linear expressions and constraints
  • typical usage pattern: model.addConstr(A @ x == b)
    • A is a Numpy ndarray, or a Scipy.sparse matrix
    • b is a Numpy ndarray

Python Matrix-Friendly API - Example


import gurobipy as gp

m = gp.Model()

x = m.addVar(ub=1.0)
y = m.addVar(ub=1.0)
z = m.addVar(ub=1.0)

m.setObjective(x*x+2*y*y+3*z*z)
m.addConstr(x + 2*y + 3*z >= 4)
m.addConstr(x + y >= 1)


import gurobipy as gp
import numpy as np

m = gp.Model()

Q = np.diag([1, 2, 3])
A = np.array([ [1, 2, 3],
             [1, 1, 0] ])
b = np.array([4, 1])

x = m.addMVar(3, ub = 1.0)
m.setObjective( x @ Q @ x)
m.addConstr( A @ x >= b)

Resources

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Callbacks


Idea:

  • Monitor the optimization progress
  • Modify the behavior of the solver

Examples:

  • Add lazy constraints or cuts
  • Define custom termination criteria
  • Add solutions to the solver

Arguments:

  • where - from where it is called (Presolve, Simplex, MIP, …)
  • what - more detailed information about the state of the optimization

Callbacks - Example


import gurobipy as gp

def mycallback(model, where):
  if where == gp.GRB.Callback.MIP:
    solcnt = model.cbGet(gp.GRB.Callback.MIP_SOLCNT)
    nodecnt = model.cbGet(gp.GRB.Callback.MIP_NODCNT)
    if nodecnt >= 10 and solcnt >= 4:
      print('Stop early')
      model.terminate()

model = gp.read("data/glass4.mps")
model.optimize(mycallback)

Callbacks - Example


    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
[...]
     0     2 8.0000e+08    0   77 2.0000e+09 8.0000e+08  60.0%     -    0s
H   64    79                    1.950016e+09 8.0000e+08  59.0%  22.5    0s
Stop early
H   74    79                    1.900016e+09 8.0000e+08  57.9%  20.7    0s
[...]
Explored 79 nodes (2138 simplex iterations) in 0.11 seconds (0.08 work units)
[...]
Solution count 5: 1.90002e+09 1.95002e+09 2.00002e+09 ... 3.13336e+09

Solve interrupted
Best objective 1.900016200000e+09, best bound 8.000042222222e+08, gap 57.8949%
User-callback calls 527, time in user-callback 0.00 sec

Callbacks

Remarks

  • Passing data to callback can be done by defining custom data attributes before optimization begins
def mycallback(model, where):
      ...
      if model._value == 1:
        ...

model._value = 1
model.optimize(mycallback)

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Infeasibility Analysis


Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (mac64[arm])
[...]
Optimize a model with 14 rows, 72 columns and 72 nonzeros
[...]

Solved in 0 iterations and 0.00 seconds (0.00 work units)
Infeasible model


Why is the model infeasible?

→ Compute an Irreducible Inconsistent Subsystem (IIS)


What changes do I need to make to recover feasibility?

→ Compute a Feasibility Relaxation: the smallest perturbation needed to recover feasibility

Irreducible Inconsistent Subsystem (IIS)


Given an infeasible system of constraints

  • Find a subset of constraints / variable bounds that:
    • Is infeasible
    • Removing a single constraint / bound makes it feasible
  • IIS is minimal (not minimum)

Meant to be read and analyzed by a human

  • The smaller, the better

Computational complexity

  • Cheap for LPs and expensive for MIPs

Irreducible Inconsistent Subsystem (IIS)

API

if model.Status == GRB.INFEASIBLE:
    model.computeIIS()
    model.write("iis.ilp")

ILP format

\ Model assignment_copy
\ LP format - for model browsing.
\ Use MPS format to capture full model detail.
Minimize

Subject To
 Thu4: x[Cathy,Thu4] + x[Ed,Thu4] = 4
Bounds
 -infinity <= x[Cathy,Thu4] <= 1
 -infinity <= x[Ed,Thu4] <= 1
End


User control to guide IIS computation: Attributes to include or exclude constraints / bounds from the IIS (IISConstrForce, IISLBForce, IISUBForce, IISSOSForce…)

Feasibility Relaxation


  • Modified model minimizing the amount by which the bounds and linear constraints of the original model are violated

  • Violation can be measured as:

    • Number of violations (0-norm)
    • Sum of the violations (1-norm)
    • Sum of squares of violations (2-norm)

Infeasible model

\[ \begin{align*} \min &~c^T x \\ \text{s. t. } &~Ax \leq b \\ &~x \geq 0 \end{align*} \]

Feasibility relaxation

\[ \begin{align*} \min &||(s, u)||_p \\ \text{s. t. } &~Ax -s \leq b \\ &~x + u \geq 0 \\ &~s, u \geq 0 \end{align*} \]

APIs

# Simplified version
model.feasRelaxS(relaxobjtype, minrelax, vrelax, crelax)

# relaxobjtype (integer): 0-norm, 1-norm, 2-norm
# minrelax (boolean): controls type of feasibility relaxation
# vrelax (boolean): indicates whether variable bounds can be relaxed
# crelax (boolean): indicates whether constraints can be relaxed

# Full version
model.feasRelax(relaxobjtype, minrelax, vars, lbpen, ubpen, constrs, rhspen)

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Start Features


  • Sometimes Gurobi has trouble finding an initial feasible solution
  • This can affect overall performance


Start Values and Variable Hints

→ Complete or partial solutions provided by the user


Start Heuristics

→ Specialized heuristics in Gurobi to produce initial solutions

Start Values and Variable Hints


Situation:

  • some variables are expected to take a certain value
  • (parts) of a feasible solution are known or can be guessed (heuristics, similar model…)

Idea: Reduce solve times by specifying these values

Example: Planning with rolling horizon: every month plan next 12 months

→ values from prior solution of first 11 months can be used as start or hints

Start Values and Variable Hints - Comparison


Start Values

  • Generate initial integer solution, which is improved via MIP search
  • Can specify partial solution, to be completed by solver


  • Controlled via Start variable attribute (or load an .MST MIP start file)
  • Supports multiple start values via NumStart model attribute and StartNumber parameter

Variable Hints

  • Guide MIP search toward anticipated values, affects entire solution process
  • Can specify hints for subset of variables, to be used by solver


  • Controlled via VarHintVal variable attribute. Optionally, express confidence with VarHintPri
  • Supports only one hint per variable

Start Heuristics


No Relaxation Heuristic

  • Doesn’t solve any node LPs of the original model
  • Controlled with NoRelHeurTime / NoRelHeurWork parameters
  • Runs before optimization to find feasible solutions. When reaching time or work limit, switch to regular branch-and-cut
    • Initially designed for cases where solving the LP is hard
    • Can help find feasible solutions fast in other cases too

Other Start Heuristics (quite expensive, generally poor quality solutions)

  • ZeroObjNodes - Number of nodes for running the zero objective heuristic
  • MinRelNodes - Number of nodes for running the minimum relaxation heuristic
  • PumpPasses - Number of passes of the feasibility pump heuristic

Hidden Gems

  • Modeling
    • General Constraints
    • Multiple Objectives
    • Multiple Scenarios
    • Multiple Solutions
    • Python Matrix-Friendly API
    • Callbacks
  • Analysis and improvements
    • Infeasibility Analysis
    • Start Features
    • gurobi-logtools

Analyze Gurobi log files with gurobi-logtools

  • Open-source Python package developed by Gurobi Experts
  • Compare different parameter sets, multiple model instances…
  • Analyze runtimes, MIPGap, …

There is more to discover!



Gurobi Optimizer Features:

GitHub:

©️ Gurobi Optimization