Solving MINLPs with gurobipy

Robert Luce

Principal Developer

Simon Bowly

Software Developer

2025-01-15

Agenda

A little background on MINLP
Outline nonlinear modelling mechanics in Gurobi 12
Cover three brief examples of NL modelling

Mixed Integer Nonlinear Programming

Why?

Application	Phenomenon	Nonlinearity
Finance	Risk	quadratic (convex)
Truss topology	Physical Forces	quadratic (convex)
Pooling (petrochemical, mining, agriculture)	Mixing Products	quadratic nonconvex
electricity distribution (ACOPF)	Alternating Current	\(\sin\) and \(\cos\) (can be made quadratic)
machine learning	—	logistic function, \(\tanh\)
chemical engineering	Chemical Reactions	anything
mechanical enginering design	Equations of Motion	ratios, angles, derivatives

With Gurobi 12 you can …

Leverage MIP technology to solve your NLPs to global optimality
Extend your discrete decision making problems with nonlinearities

… all without leaving your Pythonic armchair!

Mathematical Form

In general, we aim to solve the following problem form:

\[ \begin{aligned} & \min f(x)\\ & \text{s.t:}\\ & g_i(x) \le 0, i = 1,\ldots, m \\ & x_j \in \mathbb Z, j \in J\\ & l \le x \le u \end{aligned} \]

\(f(x)\), \(g_i(x)\) given in closed form (not via callbacks)
Mix of integer, binary, and continuous variables
Bounded search space

New Modelling Objects in `gurobipy` 12.0

New objects added to capture nonlinear expressions
You can define nonlinear constraints using Python operators

NLExpr

High level modelling class, analogous to LinExpr, QuadExpr
Constructed by arithmetic operations with other objects, or through nonlinear functions

import gurobipy as gp

env = gp.Env()
model = gp.Model(env=env)

x = model.addVar(name="x")
y = model.addVar(name="y")
model.update()

Linear expressions

x + 3.0 * y

<gurobi.LinExpr: x + 3.0 y>

Quadratic expressions

x * y + 2.0 * y

<gurobi.QuadExpr: 2.0 y + [ x * y ]>

Nonlinear expressions

x / y + 2*y

<gurobi.NLExpr: x / y + 2.0 * y>

Expression Types

gurobipy automatically promotes expressions to NLExpr whenever it is necessary to store them.

import gurobipy as gp
from gurobipy import GRB

env = gp.Env()
model = gp.Model(env=env)
x = model.addVar(name="x")
y = model.addVar(name="y")
z = model.addVar(name="z")
model.update()

2 * x                    # Linear expressions stored as LinExpr

<gurobi.LinExpr: 2.0 x>

2 * x * y                # Quadratic expressions stored as QuadExpr

<gurobi.QuadExpr: 0.0 + [ 2.0 x * y ]>

2 * x * y * z            # Higher degree polynomials promoted to NLExpr

<gurobi.NLExpr: 2.0 * x * y * z>

2 * x * y / (x + y)      # Division by non-constant: NLExpr

<gurobi.NLExpr: (2.0 * x * y) / (x + y)>

Supported Operations

Any Python operator:

add: +
subtract: -
multiply: *
divide: /
power: **

On any combination of types:

Constant
Var
LinExpr
QuadExpr
NLExpr

((2 ** x) + y / z) / (x * y - 2.0)

<gurobi.NLExpr: (2.0 ** x + y / z) / (-2.0 + x * y)>

Supported Functions

Nonlinear functions live in the nlfunc module:

sin
cos
tan
exp
sqrt

log
log2
log10
logistic
square

The result of an nlfunc is always a nonlinear expression.

from gurobipy import nlfunc

nlfunc.sin(x)

<gurobi.NLExpr: sin(x)>

nlfunc.sin(x + y)

<gurobi.NLExpr: sin(x + y)>

nlfunc.log(nlfunc.sqrt(x)) + y / z

<gurobi.NLExpr: log(sqrt(x)) + y / z>

Matrix-friendly objects

MNLExpr is to NLExpr as MLinExpr is to LinExpr

X = model.addMVar((2, 3), name="X")
Y = model.addMVar((3,), name="X")
model.update()

Follows numpy broadcasting rules:

X / Y

<MNLExpr (2, 3)>
array([[ X[0,0] / X[0],  X[0,1] / X[1],  X[0,2] / X[2]],
       [ X[1,0] / X[0],  X[1,1] / X[1],  X[1,2] / X[2]]])

Nonlinear functions are applied element-wise:

nlfunc.exp(X)

<MNLExpr (2, 3)>
array([[ exp(X[0,0]),  exp(X[0,1]),  exp(X[0,2])],
       [ exp(X[1,0]),  exp(X[1,1]),  exp(X[1,2])]])

Nonlinear constraints are general constraints

In Gurobi, nonlinear constraints are implemented as general constraints. They always take the form

\[ y = f(x) \]

with \(y\) being the resultant variable, and \(x\) being the argument variables.

To add a nonlinear inequality constraint \(l \le f(x) \le u\) do: \[ y = f(x), \; l \le y \le u \]

To minimize a nonlinear function do: \[ \min y \quad \text{s.t.} \; y = f(x), \; -\infty \le y \le \infty \]

Example: MINLPLib

Nonlinear modelling syntax
Points of interest in the solver log

HS62

\[ \begin{aligned} & \min && -32.174 (\\ & && 255 \log\left(\frac{0.03 + x_2 + x_3 + x_4}{0.03 + 0.09 x_2 + x_3 + x_4}\right) +\\ & && 280 \log\left(\frac{0.03 + x_3 + x_4}{0.03 + 0.07 x_3 + x_4}\right) +\\ & && 290 \log\left(\frac{0.03 + x_4}{0.03 + 0.13 x_4}\right)\\ & && )\\ & \text{s.t.} && 20 (x_2 + x_3 + x_4 - 1)^2 = 0\\ & && x_2, x_3, x_4 \ge 0\\ \end{aligned} \]

# Add variables to a new model
model = gp.Model(env=env)
x2 = model.addVar(name="x2", lb=0.0)
x3 = model.addVar(name="x3", lb=0.0)
x4 = model.addVar(name="x4", lb=0.0)

Objective Function

Introduce an unbounded variable to model the objective using a constraint:

objvar = model.addVar(lb=-float('inf'), name="objvar")
nle = -32.174 * (
    255 * nlfunc.log((0.03 + x2 + x3 + x4) / (0.03 + 0.09 * x2 + x3 + x4))
    + 280 * nlfunc.log((0.03 + x3 + x4) / (0.03 + 0.07 * x3 + x4))
    + 290 * nlfunc.log((0.03 + x4) / (0.03 + 0.13 * x4))
)
model.addGenConstrNL(objvar, nle)
model.setObjective(objvar, sense=GRB.MINIMIZE)

Constraints

Introduce a fixed variable to model \(f(x) = c\)

\[ \begin{aligned} & \text{s.t.} && 20 (x_2 + x_3 + x_4 - 1)^2 = 0\\ \end{aligned} \]

resvar = model.addVar(lb=0, ub=0, name="resvar")
nle = 20 * nlfunc.square((-1) + x2 + x3 + x4)
model.addGenConstrNL(resvar, nle)

Note

This could also be written directly as a quadratic constraint (recommended).

model.addConstr(((-1) + x2 + x3 + x4) ** 2 == 0)

Gurobi will handle either form via presolve reductions.

LP File Format

Use LP file format to inspect models when necessary
Internally, nonlinear expressions are expression trees
Readable form expressions are shown as a comment

model.write("HS62.lp")

\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
  0 x2 + 0 x3 + 0 x4 + objvar + 0 resvar
Subject To
Bounds
 objvar free
 resvar = 0
General Constraints
\ objvar = -32.174 * (((255 * log((0.03 + x2 + x3 + x4) / (0.03 + (0.09 * 
\ x2) + x3 + x4))) + (280 * log((0.03 + x3 + x4) / (0.03 + (0.07 * x3) + 
\ x4)))) + (290 * log((0.03 + x4) / (0.03 + (0.13 * x4)))))
 GC0: objvar = NL : ( MULTIPLY , -1 , -1 ) ( CONSTANT , -32.174 , 0 ) 
  ( PLUS , -1 , 0 ) ( PLUS , -1 , 2 ) ( MULTIPLY , -1 , 3 ) 
  ( CONSTANT , 255 , 4 ) ( LOG , -1 , 4 ) ( DIVIDE , -1 , 6 ) 
  ( PLUS , -1 , 7 ) ( CONSTANT , 0.03 , 8 ) ( VARIABLE , x2 , 8 ) 
  ( VARIABLE , x3 , 8 ) ( VARIABLE , x4 , 8 ) ( PLUS , -1 , 7 ) 
  ( CONSTANT , 0.03 , 13 ) ( MULTIPLY , -1 , 13 ) ( CONSTANT , 0.09 , 15 ) 
  ( VARIABLE , x2 , 15 ) ( VARIABLE , x3 , 13 ) ( VARIABLE , x4 , 13 ) 
  ( MULTIPLY , -1 , 3 ) ( CONSTANT , 280 , 20 ) ( LOG , -1 , 20 ) 
  ( DIVIDE , -1 , 22 ) ( PLUS , -1 , 23 ) ( CONSTANT , 0.03 , 24 ) 
  ( VARIABLE , x3 , 24 ) ( VARIABLE , x4 , 24 ) ( PLUS , -1 , 23 ) 
  ( CONSTANT , 0.03 , 28 ) ( MULTIPLY , -1 , 28 ) ( CONSTANT , 0.07 , 30 ) 
  ( VARIABLE , x3 , 30 ) ( VARIABLE , x4 , 28 ) ( MULTIPLY , -1 , 2 ) 
  ( CONSTANT , 290 , 34 ) ( LOG , -1 , 34 ) ( DIVIDE , -1 , 36 ) 
  ( PLUS , -1 , 37 ) ( CONSTANT , 0.03 , 38 ) ( VARIABLE , x4 , 38 ) 
  ( PLUS , -1 , 37 ) ( CONSTANT , 0.03 , 41 ) ( MULTIPLY , -1 , 41 ) 
  ( CONSTANT , 0.13 , 43 ) ( VARIABLE , x4 , 43 )
\ resvar = 20 * sqr(-1 + x2 + x3 + x4)
 GC1: resvar = NL : ( MULTIPLY , -1 , -1 ) ( CONSTANT , 20 , 0 ) 
  ( SQUARE , -1 , 0 ) ( PLUS , -1 , 2 ) ( CONSTANT , -1 , 3 ) 
  ( VARIABLE , x2 , 3 ) ( VARIABLE , x3 , 3 ) ( VARIABLE , x4 , 3 )
End

Solver Log

model.optimize()

Optimize a model with 0 rows, 5 columns and 0 nonzeros
Model fingerprint: 0x4a8801ad
Model has 2 general nonlinear constraints (7 nonlinear terms)
Variable types: 5 continuous, 0 integer (0 binary)

...

Added 15 variables to disaggregate expressions.
Presolve time: 0.00s
Presolved: 83 rows, 27 columns, 165 nonzeros
Presolved model has 12 bilinear constraint(s)
Presolved model has 6 nonlinear constraint(s)

You can model NL expressions naturally in gurobipy. Presolve will transform them to suit Gurobi.

...

Warning: Model contains variables with very large bounds participating
         in nonlinear terms.
         Presolve was not able to compute smaller bounds for these variables.
         Consider bounding these variables or reformulating the model.

Important

Mind the warnings! Bounds are always advisable in MINLP.

Branch and Bound Log

...

Solving non-convex MINLP

Variable types: 27 continuous, 0 integer (0 binary)

Root relaxation: objective -7.354055e+04, 25 iterations, 0.00 seconds (0.00 work units)

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

     0     0 -73540.546    0   14          - -73540.546      -     -    0s
     0     0 -68838.434    0   11          - -68838.434      -     -    0s
     0     0 -56834.224    0   15          - -56834.224      -     -    0s
H    0     0                    -26272.51449 -56834.224   116%     -    0s
     0     2 -56834.224    0   15 -26272.514 -56834.224   116%     -    0s

Explored 597 nodes (1851 simplex iterations) in 0.04 seconds (0.02 work units)

...

Optimal solution found (tolerance 1.00e-04)
Best objective -2.627251448732e+04, best bound -2.627276807630e+04, gap 0.0010%

Example: Optimal Design

Application from mechanical engineering
Using a global solver gives us some extra tricks

Compound Gear Train Design

Design a compound gear train with a gear ratio of \(6.931\)
Each gear must have between 12 and 60 teeth
Driveshafts and gear body must not interfere

Equations of Motion

Relate velocities at gear pitch diameter:

\[ \begin{aligned} \omega_d T_d & = - \omega_a T_a & \\ \omega_b T_b & = - \omega_f T_f & \\ \omega_a & = \omega_b \\ \end{aligned} \]

Find the compound gear ratio:

\[ \begin{aligned} \frac{\omega_d}{\omega_f} & = \frac{\omega_a T_a}{-T_d} \frac{-T_f}{\omega_b T_b} \\ & = \frac{T_f T_a}{T_b T_d} = 6.931 \\ \end{aligned} \]

This is a nonlinear relationship in our integer design variables

Physical Design Restrictions

Model each gear by its number of teeth

Integer variables, bounded in \([12, 60]\)

Add gear diameter restrictions to avoid collisions

\[ T_a \le T_b + T_f \]

`gurobipy` model

gears = ["d", "a", "b", "f"]
min_teeth = 12
max_teeth = 60
target_ratio = 6.931

model = gp.Model(env=env)

T = model.addVars(
    gears,
    lb=min_teeth, ub=max_teeth,
    vtype=GRB.INTEGER,
    name="T",
)

gear_ratio = (T["f"] * T["a"]) / (T["b"] * T["d"])  # NLExpr

z = model.addVar()
model.addGenConstrNL(z, (gear_ratio - target_ratio) ** 2)
model.setObjective(z, sense=GRB.MINIMIZE)

model.addConstr(T["a"] <= T["b"] + T["f"])

Global Optimization: Find Multiple Designs

model.Params.PoolSolutions = 20  # Find 20 solutions
model.Params.PoolSearchMode = 2  # Find the 20 best solutions

model.optimize()

Added 4 variables to disaggregate expressions.
Presolve time: 0.00s
Presolved: 16 rows, 10 columns, 44 nonzeros
Presolved model has 4 bilinear constraint(s)

...
Explored 37728 nodes (86875 simplex iterations) in 0.20 seconds (0.05 work units)
Thread count was 11 (of 11 available processors)

Solution count 20: 6.23269e-09 6.23269e-09 6.23269e-09 ... 2.29353e-07
No other solutions better than 2.29353e-07

The objective here is squared error
We have 20 solutions with the correct gear ratio to 3 decimal places

Solutions

Example: AC Optimal Power Flow

Complex model to solve (covered in previous webinars)
Brief example of NL modelling using matrix variables

ACOPF: Power Flow Equations

Power system modelled as a connected network with:

Generators (supply nodes)
Loads (demand nodes)
Intermediate buses
Power flows following laws of physics

\[ P_i = v_i \sum_{j=1}^N v_j \left( G_{ij} \cos(\theta_i - \theta_j) + B_{ij} \sin(\theta_i - \theta_j) \right) \]

\[ Q_i = v_i \sum_{j=1}^N v_j \left( G_{ij} \sin(\theta_i - \theta_j) + B_{ij} \cos(\theta_i - \theta_j) \right) \]

Matrix Variables

import gurobipy as gp
from gurobipy import GRB
from gurobipy import nlfunc

env = gp.Env()
model = gp.Model("OptimalPowerFlow", env=env)

P = model.addMVar(N, name="P", lb=P_lb, ub=P_ub)  # Real power for buses
Q = model.addMVar(N, name="Q", lb=Q_lb, ub=Q_ub)  # Reactive power for buses
v = model.addMVar(N, name="v", lb=v_lb, ub=v_ub)  # Voltage magnitude at buses

# Voltage angle at buses
theta = model.addMVar(N, name="theta", lb=theta_lb, ub=theta_ub).reshape(N, 1)

theta - theta.T

<MLinExpr (4, 4)>
array([[ theta[0] + -1.0 theta[0],  theta[0] + -1.0 theta[1],
         theta[0] + -1.0 theta[2],  theta[0] + -1.0 theta[3]],
       [ theta[1] + -1.0 theta[0],  theta[1] + -1.0 theta[1],
         theta[1] + -1.0 theta[2],  theta[1] + -1.0 theta[3]],
       [ theta[2] + -1.0 theta[0],  theta[2] + -1.0 theta[1],
         theta[2] + -1.0 theta[2],  theta[2] + -1.0 theta[3]],
       [ theta[3] + -1.0 theta[0],  theta[3] + -1.0 theta[1],
         theta[3] + -1.0 theta[2],  theta[3] + -1.0 theta[3]]])

Power Balance Equations

\[ G_{ij} \cos(\theta_i - \theta_j) \]

G * nlfunc.cos(theta - theta.T)

<MNLExpr (4, 4)>
array([[ 1.7647 * cos(theta[0] + -1.0 * theta[0]),
         -0.5882 * cos(theta[0] + -1.0 * theta[1]),
         0.0 * cos(theta[0] + -1.0 * theta[2]),
         -1.1765 * cos(theta[0] + -1.0 * theta[3])],
       [ -0.5882 * cos(theta[1] + -1.0 * theta[0]),
         1.5611 * cos(theta[1] + -1.0 * theta[1]),
         -0.3846 * cos(theta[1] + -1.0 * theta[2]),
         -0.5882 * cos(theta[1] + -1.0 * theta[3])],
       [ 0.0 * cos(theta[2] + -1.0 * theta[0]),
         -0.3846 * cos(theta[2] + -1.0 * theta[1]),
         1.5611 * cos(theta[2] + -1.0 * theta[2]),
         -1.1765 * cos(theta[2] + -1.0 * theta[3])],
       [ -1.1765 * cos(theta[3] + -1.0 * theta[0]),
         -0.5882 * cos(theta[3] + -1.0 * theta[1]),
         -1.1765 * cos(theta[3] + -1.0 * theta[2]),
         2.9412 * cos(theta[3] + -1.0 * theta[3])]])

Power Balance Equations

\[ \sum_{j=1}^N v_j G_{ij} \cos(\theta_i - \theta_j) \]

nle = (G * nlfunc.cos(theta - theta.T)) @ v
print(nle[0].item())

1.7647 * cos(theta[0] + -1.0 * theta[0]) * v[0] + -0.5882 * cos(theta[0] + -1.0 * theta[1]) * v[1] + 0.0 * cos(theta[0] + -1.0 * theta[2]) * v[2] + -1.1765 * cos(theta[0] + -1.0 * theta[3]) * v[3]

\[ P_i = v_i \sum_{j=1}^N v_j \left( G_{ij} \cos(\theta_i - \theta_j) + B_{ij} \sin(\theta_i - \theta_j) \right) \]

constr_P = model.addGenConstrNL(
    P,
    v * ((G * nlfunc.cos(theta - theta.T) + B * nlfunc.sin(theta - theta.T)) @ v),
    name="constr_P",
)

Single call adds power balance constraints for the entire network!

Final Words

Nonlinear modelling is dead easy in gurobipy
Global NLP solver integrates with solution pools
Vectorization / matrix notation for larger models

Add mixed integer nonlinear global optimization to your toolbox!

Solving MINLPs with gurobipy

Agenda

Mixed Integer Nonlinear Programming

Why?

With Gurobi 12 you can …

Mathematical Form

New Modelling Objects in gurobipy 12.0

NLExpr

Expression Types

Supported Operations

Supported Functions

Matrix-friendly objects

Nonlinear constraints are general constraints

Example: MINLPLib

HS62

Objective Function

Constraints

LP File Format

Solver Log

Branch and Bound Log

Example: Optimal Design

Compound Gear Train Design

Equations of Motion

Physical Design Restrictions

gurobipy model

Global Optimization: Find Multiple Designs

Solutions

Example: AC Optimal Power Flow

ACOPF: Power Flow Equations

Matrix Variables

Power Balance Equations

Power Balance Equations

Final Words

Questions

New Modelling Objects in `gurobipy` 12.0

`gurobipy` model