gurobipy: Exercices

The gurobipy card game

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 characteristics:

A file .python-version 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,

3.12

A file requirements.txt 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

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)

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)

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)

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} + \sigma_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()

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…