gurobipy: Exercices

1. The gurobipy card game

2. Creating an easy-to-use project

The goal for this exercise is to create a trivial Python script that I, as your examiner, will be able to run automatically.

The result of this exercise is the URL to a (public!) GitHub repo that should have the following files in the root of the repository:

A text file .python-version (in UTF-8 encoding) exists that contains the Python version that you used to create this project. If you use pyenv and ran pyenv local …, this file was created for you. Otherwise, create it manually. For example,

python -c "
import sys;
print(f'{sys.version_info.major}.'
      f'{sys.version_info.minor}')
" > .python-version

A text file requirements.txt (in UTF-8 encoding) exists with the Python libraries that are required for the script below to run correctly. See e.g. pip freeze documentation.
A file easy.py exists in this project, with this content:

import gurobipy as gp

parameters = {
    "OutputFlag": 0
}

with gp.Env(params=parameters) as env, gp.Model(env=env) as m:
    m.optimize()
print(gp.GRB.VERSION_MAJOR)

Send to xnodet@uco.fr the https URL to your repo with your name, and I will run the following script to check that your submission is valid.

#!/bin/bash -e

url=$1
name=$2
branch=$3

git clone -b $branch $url repo-$name
cd repo-$name
python -m venv venv --prompt "repo-$name"
source venv/bin/activate
pip install --requirement requirements.txt
echo 12 > correct-output.txt
python easy.py > output.txt
diff correct-output.txt output.txt
echo SUCCESS

3. Knapsack

Given a set of $n$ items, each with a value $r_i$ and weight $c_i$, find the optimal allocation of items to a knapsack with capacity $C$ to maximize the total value collected.

\[ x_i = \begin{cases} 1 & \text{if item } i \text{ is selected} \\ 0 & \text{otherwise} \end{cases} \] \[ \text{maximize} \quad \sum_{i=1}^n r_i x_i \] \[ \text{subject to} \quad \sum_{i=1}^n c_i x_i \leq C \] \[ x_i \in \{0, 1\}, \quad i = 1, \dots, n \]

import numpy as np
import gurobipy as gp
from gurobipy import GRB

def generate_knapsack(num_items):
    # Fix seed value
    rng = np.random.default_rng(seed=0)
    # Item values, weights
    values = rng.uniform(low=1, high=25, size=num_items)
    weights = rng.uniform(low=5, high=100, size=num_items)
    # Knapsack capacity
    capacity = 0.7 * weights.sum()

    return values, weights, capacity


def solve_knapsack_model(values, weights, capacity):
    num_items = len(values)
    # Turn values and weights numpy arrays to dict
    # ...

    with gp.Env() as env:
        with gp.Model(name="knapsack", env=env) as model:
            # Define decision variables using the Model.addVars() method
            # ...

            # Define objective function using the Model.setObjective() method
            # Build the LinExpr using the tupledict.prod() method
            # ...

            # Define capacity constraint using the Model.addConstr() method
            # ...

            model.optimize()

data = generate_knapsack(10000)
solve_knapsack_model(*data)

4. Portfolio Optimization

In a portfolio optimization problem, there are $n$ assets. Each asset $i$ is associated with an expected return $\mu_i$ and each pair of assets $(i,j)$ has a covariance (risk) $\sigma_{ij}$. The goal is to find the optimal fraction of the portfolio invested in each asset to minimize the risk of investment such that 1) the total expected return of the investment exceeds the minimum target return $\mu_0$ and 2) the portfolio invests in at most $k \le n$ assets.

$x_i$: Relative investment in asset
$y_i$: Binary variable controlling whether asset $i$ is traded

The mathematical formulation of the model is therefore:

\[\begin{align} minimize & \sum_{i, j}x_i \sigma_{ij}x_j \\ s.t. & \sum_{i, j}\mu_i x_i \geq \mu_0 \\ & \sum_{i, j}x_i = 1 \\ & \sum_{i, j}y_i \leq k \\ & x_i \leq y_i, \forall i \\ \end{align}\]

import json
import pandas as pd
import numpy as np
import gurobipy as gp
from gurobipy import GRB

with open("data/portfolio-example.json", "r") as f:
    data = json.load(f)

n = data["num_assets"]
sigma = np.array(data["covariance"])
mu = np.array(data["expected_return"])
mu_0 = data["target_return"]
k = data["portfolio_max_size"]


with gp.Model("portfolio") as model:
    # Name the modeling objects to retrieve them
    # ...

    model.optimize()

    # Write the solution into a DataFrame
    portfolio = [var.X for var in model.getVars() if "x" in var.VarName]
    risk = model.ObjVal
    expected_return = model.getRow(model.getConstrByName("return")).getValue()
    df = pd.DataFrame(
        data=portfolio + [risk, expected_return],
        index=[f"asset_{i}" for i in range(n)] + ["risk", "return"],
        columns=["Portfolio"],
    )
    print(df)

5. Multi-Period Lot-Sizing with Setups and Minimum Batch

In this production-planning problem, a single product faces known demands over $T$ periods. Each period has a linear production cost and a fixed setup cost if production occurs. When you produce in period $t$, you must make at least $Q_{\min}$ units and at most $Q_{\max}$ units. Inventory can be carried between periods with a holding cost. The goal is to meet demand on time at minimum total cost.

$x_t$: Production quantity in period $t$
$y_t$: Binary variable indicating whether production occurs in period $t$
$I_t$: End-of-period inventory after meeting demand $d_t$

Parameters:

$d_t$: Demand in period $t$
$c_t$: Variable production cost (€/unit)
$f_t$: Fixed setup cost (€/setup)
$h_t$: Holding cost (€/unit carried to $t+1$)
$Q_{\min}$: Minimum batch size if producing
$Q_{\max}$: Maximum production capacity
$I_0$: Initial inventory

The mathematical formulation of the model is therefore:

\[\begin{align} \min \quad & \sum_{t=0}^{T-1} \left(c_t x_t + f_t y_t + h_t I_t\right) \\[6pt] \text{s.t.} \quad & I_{t-1} + x_t - d_t = I_t && \forall t = 0,\dots,T-1 \quad (\text{with } I_{-1} = I_0) \\[6pt] & x_t \le Q_{\max} y_t && \forall t = 0,\dots,T-1 \\[4pt] & x_t \ge Q_{\min} y_t && \forall t = 0,\dots,T-1 \\[6pt] & x_t \ge 0,\; I_t \ge 0 && \forall t = 0,\dots,T-1 \\[4pt] & y_t \in \{0,1\} && \forall t = 0,\dots,T-1 \end{align}\]

import json
import gurobipy as gp
from gurobipy import GRB
from pathlib import Path

# ----- Load data from JSON -----
with open("data/lot_sizing_data.json", "r") as f:
    data = json.load(f)

name = data["name"]
H    = int(data["H"])
d    = [float(val) for val in data["demand"]]
c    = [float(val) for val in data["var_cost"]]
f    = [float(val) for val in data["setup_cost"]]
h    = [float(val) for val in data["hold_cost"]]
Qmin = float(data["Qmin"])
Qmax = float(data["Qmax"])
I0   = float(data["I0"])

# Basic validation
assert len(d) == H and len(c) == H and len(f) == H and len(h) == H
assert 0 <= Qmin <= Qmax

# ----- Build model -----
with gp.Env() as env, gp.Model(name, env=env) as model:

    # ...

    # Optimize
    model.optimize()

    if model.SolCount:
        assert model.ObjVal == 1198.5
        print(f"Total cost = {model.ObjVal:.2f}")
        for t in range(H):
            print(f"t={t:2d}: y={int(y[t].X)} x={x[t].X:.1f} I={I[t].X:.1f}")

6. Bucket Design with Fixed Surface

We want to design a bucket in the shape of a truncated cone (frustum) and maximize its volume. The bucket is open at the top, with:

bottom radius $r \ge 0$,
top radius $R \ge 0$,
height $h \ge 0$.

The bucket has a bottom disk but no lid. Its volume is:

\[V = \frac{\pi h}{3}\,(R^2 + Rr + r^2).\]

The amount of material available is limited: the sum of the bottom area and the lateral surface area must be exactly 1 square meter.

Bottom area:

\[A_{\text{bot}} = \pi r^2.\]

Lateral area:

\[A_{\text{lat}} = \pi (R + r)\,\sqrt{(R-r)^2 + h^2}.\]

Material constraint:

\[A_{\text{bot}} + A_{\text{lat}} = 1.\]

7. Custom Termination Criteria

Use a callback function to terminate the optimization of the model mkp.mps.bz2 (in the data.zip file alongside the list of exercises) if the MIPGap does not decrease more than epsilon (1e-4) for a given time period (15 seconds) after finding the first incumbent solution.

Hint: Check the possible what values for the callback code GRB.Callback.MIP.

from functools import partial

import gurobipy as gp
from gurobipy import GRB


class CallbackData:
    def __init__(self):
        self.last_gap_change_time = -GRB.INFINITY
        self.last_gap = GRB.INFINITY


def callback(model, where, *, cbdata):
    if where != GRB.Callback.MIP:
        return
    if model.cbGet(GRB.Callback.MIP_SOLCNT) == 0:
        return

    # Use model.terminate() to end the search when required...
    # ...


with gp.read("data/mkp.mps.bz2") as model:
    # Global variables used in the callback function
    max_time_between_gap_updates = 15
    epsilon_to_compare_gap = 1e-4

    # Initialize data passed to the callback function
    callback_data = CallbackData()
    callback_func = partial(callback, cbdata=callback_data)

    model.optimize(callback_func)

8. Unit Commitment Problem

Unit Commitment is a classic problem in Power System that has been extensively studied. Lately, the increased interest in renewable generation (solar and wind) and the development of large storage systems have contributed to its use by system operators on grids of different sizes.

Sets

$T$ Time intervals (in hours)
$I$ Generator units

Parameters

$\bar{\bar{p_i}}$ Maximum Power for unit $i$ (MW)
$\bar{p_i}$ Minimum Power for unit $i$ (MW)
$\alpha_i$ Fixed Cost for unit $i$ ($/h)
$\beta_i$ Linear Cost for unit $i$ ($/MWh)
$\gamma_i$ Quadratic Cost for unit $i$ ($/MWh/MWh)
$\delta_i$ Start Up Cost for unit $i$ ($)
$\zeta_i$ Shutdown Cost for unit $i$ ($)
$L_t$ Load Forecast at time interval $t$ (MW)
$S_t$ Solar Energy Forecast at time interval $t$ (MW)

Variables

$p_{i,t} \in \mathbb{R^+}$ Output power for unit $i$ at time interval $t$ (MW )
$u_{i,t} \in \{0,1\}$ Commitment Status for thermal unit $i$ at time interval $t$ (1: Online, 0: Otherwise))
$v_{i,t} \in \{0,1\}$ Startup Status for thermal unit $i$ at time interval $t$ (1: Start Up, 0: Otherwise)
$w_{i,t} \in \{0,1\}$ Shutdown Status for thermal unit $i$ at time interval $t$ (1: Shutdown, 0: Otherwise)

Objective

The objective function of the unit commitment problem is to minimize the total costs of operating thermal generators over a time horizon:

\[\begin{align} minimize \sum_{i \in I}\sum_{t \in T}(\gamma_i p^2_{i,t} + \beta_i p_{i,t} + \alpha_i u_{i,t} + \delta_i v_{i,t} + \zeta_i w_{i,t}) \end{align}\]

Power balance

All generating units have to fulfill the load forecast not covered by the solar power for all time intervals:

\[\begin{align} \sum_{i \in I}p_{i,t} + S_t = L_t, \forall t \in T \end{align}\]

Thermal Units Operating Limits

Thermal units are subject to physical operating constraints when the unit is online:

\[\begin{align} u_{i,t} \bar{p_i} \leq p_{i,t} \leq u_{i,t} \bar{\bar{p_i}}, \forall i \in I, \forall t \in T \end{align}\]

Thermal Units Logical Constraints

Thermal units on/off statuses must relate logically to start up and shutdown statuses:

\[\begin{align} u_{i,t} - u_{i,t-1} = v_{i,t} - w_{i,t}, \forall i \in I, \forall t \in T \end{align}\]

A thermal unit cannot start up and shutdown on the same time interval:

\[\begin{align} v_{i,t} + w_{i,t} \leq 1, \forall i \in I, \forall t \in T \end{align}\]

import gurobipy as gp
from gurobipy import GRB


# 24 Hour Load Forecast (MW)
load_forecast = [
     4,  4,  4,  4,  4,  4,   6,   6,
    12, 12, 12, 12, 12,  4,   4,   4,
     4, 16, 16, 16, 16,  6.5, 6.5, 6.5,
]

# solar energy forecast (MW)
solar_forecast = [
    0,   0,   0,   0,   0,   0,   0.5, 1.0,
    1.5, 2.0, 2.5, 3.5, 3.5, 2.5, 2.0, 1.5,
    1.0, 0.5, 0,   0,   0,   0,   0,   0,
]

# global number of time intervals
nTimeIntervals = len(load_forecast)

# thermal units
thermal_units = ["gen1", "gen2", "gen3"]

# thermal units' costs  (a + b*p + c*p^2), (startup and shutdown costs)
thermal_units_cost, a, b, c, sup_cost, sdn_cost = gp.multidict(
    {
        "gen1": [5.0, 0.5, 1.0, 2, 1],
        "gen2": [5.0, 0.5, 0.5, 2, 1],
        "gen3": [5.0, 3.0, 2.0, 2, 1],
    }
)

# thernal units operating limits
thermal_units_limits, pmin, pmax = gp.multidict(
    {"gen1": [1.5, 5.0], "gen2": [2.5, 10.0], "gen3": [1.0, 3.0]}
)

# thermal units dynamic data (initial commitment status)
thermal_units_dyn_data, init_status = gp.multidict(
    {"gen1": [0], "gen2": [0], "gen3": [0]}
)


def show_results():
    obj_val_s = model.ObjVal
    print(f" OverAll Cost = {round(obj_val_s, 2)}	")
    print("\n")
    print("%5s" % "time", end=" ")
    for t in range(nTimeIntervals):
        print("%4s" % t, end=" ")
    print("\n")

    for g in thermal_units:
        print("%5s" % g, end=" ")
        for t in range(nTimeIntervals):
            print("%4.1f" % thermal_units_out_power[g, t].X, end=" ")
        print("\n")

    print("%5s" % "Solar", end=" ")
    for t in range(nTimeIntervals):
        print("%4.1f" % solar_forecast[t], end=" ")
    print("\n")

    print("%5s" % "Load", end=" ")
    for t in range(nTimeIntervals):
        print("%4.1f" % load_forecast[t], end=" ")
    print("\n")


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

    # add variables for thermal units (power and statuses for commitment, startup and shutdown)
    thermal_units_out_power = model.addVars(
        # ...
        name="thermal_units_out_power"
    )
    thermal_units_startup_status = model.addVars(
        # ...
        name="thermal_unit_startup_status",
    )
    thermal_units_shutdown_status = model.addVars(
        # ...
        name="thermal_unit_shutdown_status",
    )
    thermal_units_comm_status = model.addVars(
        # ...
        name="thermal_unit_comm_status"
    )

    # define objective function as an empty quadratic construct and add terms
    obj_fun_expr = gp.QuadExpr(0)
    for t in range(nTimeIntervals):
        for g in thermal_units:
                    # ...
    model.setObjective(obj_fun_expr)

    # Power balance equations
    for t in range(nTimeIntervals):
        model.addConstr(
            # ...
            name="power_balance_" + str(t),
        )

    # Thermal units logical constraints
    for t in range(nTimeIntervals):
        for g in thermal_units:
            if t == 0:
                model.addConstr(
                    # ...
                    name="logical1_" + g + "_" + str(t),
                )
            else:
                model.addConstr(
                    # ...
                    name="logical1_" + g + "_" + str(t),
                )

            model.addConstr(
                # ...
                name="logical2_" + g + "_" + str(t),
            )

    # Thermal units physical constraints, using indicator constraints
    for t in range(nTimeIntervals):
        for g in thermal_units:
            model.addGenConstrIndicator(
                # ...
            )
            model.addGenConstrIndicator(
                # ...

            )
            model.addGenConstrIndicator(
                # ...
            )

    model.optimize()
    show_results()

9. Unit Commitment Problem (matrix API)

Change the code from the previous exercise to use the matrix API:

all variables should be created with model.addMVar
with one exception, all loops to post constraints are eliminated. The exception is that posting the indicator constraints for the physical constraints still requires a code that looks like the following:

for g in range(len(thermal_units)):
    model.addGenConstrIndicator(...)
    model.addGenConstrIndicator(...)
    model.addGenConstrIndicator(...)

The matrix API works with integer indices only, so the data needs some massaging…

10. Robotic Arm

A planar robot arm with two revolute joints has link lengths $L_1$ and $L_2$. The base is fixed at the origin. Joint angles are $\theta_1$ (shoulder) and $\theta_2$ (elbow, relative to link 1). The end-effector position $(x, y)$ is:

\[\begin{aligned} x &= L_1 \cos\theta_1 + L_2 \cos(\theta_1+\theta_2),\\ y &= L_1 \sin\theta_1 + L_2 \sin(\theta_1+\theta_2). \end{aligned}\]

We want the end-effector to reach a target point $(x^*, y^*)$ while respecting joint limits and avoiding a circular obstacle (e.g., a post) of radius $r$ centered at $(x_o, y_o)$.

We approximate collision avoidance by keeping the midpoint of the first link outside the obstacle.

import math
import gurobipy as gp
from gurobipy import GRB, nlfunc as nl


# Parameters
L1, L2 = 1.0, 0.8             # Lengths of the links
x_star, y_star = 1.20, 0.60   # Point to reach
xo, yo, r = 0.50, 0.00, 0.20  # Disk to avoid

# Joint limits
theta1_min, theta1_max = -math.pi, math.pi
theta2_min, theta2_max = -0.75*math.pi, 0.75*math.pi


# Build model
m = gp.Model("robot_arm_nlfunc")

# Decision variables (angles)
theta1 = m.addVar(lb=theta1_min, ub=theta1_max, name="theta1")
theta2 = m.addVar(lb=theta2_min, ub=theta2_max, name="theta2")

...

m.optimize()

sol = None
if m.Status == GRB.OPTIMAL:
    sol = {
        "theta1": theta1.X, "theta2": theta2.X,
        "x": x.X, "y": y.X, "xm": xm.X, "ym": ym.X,
        "obj": m.ObjVal,
    }
    print("Optimal objective:", m.ObjVal)
    print(sol)
else:
    print("Optimization status:", m.Status)


import matplotlib.pyplot as plt

def draw_arm(ax, L1, L2, th1, th2, xo, yo, r, x_star, y_star, title):
    x1 = L1*math.cos(th1); y1 = L1*math.sin(th1)
    x2 = x1 + L2*math.cos(th1 + th2); y2 = y1 + L2*math.sin(th1 + th2)

    ax.plot([0, x1], [0, y1], linewidth=3)
    ax.plot([x1, x2], [y1, y2], linewidth=3)
    ax.scatter([0, x1, x2], [0, y1, y2], s=40)

    # obstacle
    t = [i*2*math.pi/300 for i in range(301)]
    cx = [xo + r*math.cos(tt) for tt in t]
    cy = [yo + r*math.sin(tt) for tt in t]
    ax.plot(cx, cy, linewidth=2)

    # target
    ax.scatter([x_star], [y_star], marker='x', s=80)

    ax.set_aspect('equal', adjustable='box')
    ax.set_xlim(-0.5, L1+L2+0.2)
    ax.set_ylim(-0.5, L1+L2+0.2)
    ax.grid(True, linestyle=':')
    ax.set_title(title)

if sol is not None:
    fig, ax = plt.subplots(figsize=(6,6))
    draw_arm(ax, L1, L2, sol['theta1'], sol['theta2'], xo, yo, r, x_star, y_star,
             title=f"Robot Arm (nlfunc)\nobj={sol['obj']:.4g}")
    # Save as PNG instead of showing
    plt.savefig("images/robot-arm.png", dpi=100, bbox_inches="tight")
    plt.close(fig)
else:
    print("No solution available to plot yet.")