Caleb Ju

This post will go over how to setup a simple LP with Gurobi using a matrix-friendly approach, i.e., utilize numpy. This post is a ~5min read. For the full code, see this pastebin.

I am running Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64).

We will solve the linear program (LP), $$\begin{aligned} \min~& c^Tx\\ \text{s.t.}~& Ax \leq b \\ & x \in \mathbb{R}^n. \end{aligned}$$ Our focus is a simple economic dispatch problem where $n = 3$, and we have the demand and limit constraints, $$\begin{aligned} x_1 + x_2 &\geq d_1 \\ x_1 + x_3 &\geq d_2 \\ x_2 + x_3 &\geq d_3 \\ x_0 &\geq \mathrm{lb}_0 \\ x_i &\leq \mathrm{ub}_i, \ \forall i \in \{1,\ldots,n\}. \end{aligned}$$ Here, variables $x_i$ will be the thermal generator for region $i$. We assume region(s) can supply energy to the various demands, denoted by $d_i$. The lower bound is included for the sake of exploring Gurobi.

If the terminology above is alien, feel free to (safely) ignore the particular problem and focus on the general LP instead.

We begin by setting up and solving the model. First, note the use of addMVars() (where the M stands for matrix) over addVar() or addVars(). I like the matrix version since it returns a numpy array, which is useful for performing matrix and vector operations. Second, note the generous naming of variables and constraints. I believe this is (optional) good practice for when you want to extract or modify specific parts of the MILP. Third, the constants are arbitrary and are chosen to ensure the MILP is feasible.

import numpy as np
import gurobipy as gp

# define Gurobi model
m = gp.Model("Unit Commitment")

n = 3
ubs = [100, 90, 80]
demands = [185, 150, 150]
lb = 90

# define variables
x = m.addMVar(n, ub=ubs, name="x")

# add constraints and names
m.addConstr(sum(x[0], x[1]) >= demands[0], name="demand_0")
m.addConstr(sum(x[0], x[2]) >= demands[1], name="demand_1")
m.addConstr(sum(x[1], x[2]) >= demands[2], name="demand_2")
m.addConstr(x[0] >= lb, name="lb")

# define objective function
c = np.array([10, 8, 9])
m.setObjective(c@x, gp.GRB.MINIMIZE)

# solve
m.optimize()

To see how the LP performed, we want to first check if the LP was feasible or not (more statuses listed in here). From there, one can extract the primal and dual (assuming an LP and not mixed-integer LP) solutions and the objective function, as shown below.

# check status
if m.status == gp.GRB.INFEASIBLE:
    print("infeasible")
elif m.status == gp.GRB.OPTIMAL:
    print("feasible")

    x_sol = x.X
    y_sol = np.array([
        m.getAttr("Pi", [m.getConstrByName("demand_{}".format(i))])[0]
        for i in range(n)
    ])
    val = m.objVal

Now, let us say we want to two modifications: a new generator and relax the lower bound constraint.

This new generator is the variable $x_4$. Let's assume the fourth generator has a cost of $8 per unit of power and is connected to all regions. We will limit it to 20 units of power. We can either write a new LP from scatch or extend the previous one. In the latter case, we will extend the cost vector with an element of 8 and append a new column and row to the constraint matrix, as well as extend the RHS vector $b$ with an 8.

# gather existing constraints 
constrs_list = []
for cs in m.getConstrs():
    sense = cs.sense
    name = cs.ConstrName
    if "demand_" in name:
        constrs_list.append(cs)
print(constrs_list)
col = gp.Column(np.ones(len(constrs_list)), constrs_list)

# add new variable
x3 = m.addVar(obj=8, name='x[3]', column=col)
m.update()
x = m.getVars()

# add new upper bound constraints
m.addConstr(x[3] <= 20)

# modify RHS
m.setAttr("RHS", m.getConstrByName("lb"), 50)

Two things to note. First, we include update() to allow Gurobi to process the changes (update is useful for examining the model, but need not be called for optimizing). We do this so we can retrieve our updated list of objects. Note that when we used getVars(), we return a list of Gurobi variables rather than MVar. Therefore, we cannot run x.X like we did in Line 7 in the second code snippet to retrieve the array. Rather, we must go through each individual element, i.e., x[i].X.

When we solve a Gurobi model, there may be a superfluous amount of output for someone trying to just perform optimization. One may also seek performance improvements by using multiple threads. Here, we state four useful options. More information and other settings can be found here.

Finally, we show how to retrieve the LP itself. This may be useful when one wants to embed an LP into another LP or form the LP dual.

m.setParam('OutputFlag', 0)
m.setParam(gp.GRB.Param.Threads, 4)
m.setParam(gp.GRB.Param.TimeLimit, 60*60) # 1h time limit
m.setParam(gp.GRB.Param.Crossover, 0)
m.update()

# retrieve (c,A,b)
c = m.getAttr('Obj', m.getVars())
A = m.getA() # returns sparse array
b = m.getAttr('RHS', m.getConstrs())

A matrix-friendly approach to Gurobi