This post will go over how to setup a simple LP with Gurobi using
a matrixfriendly 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.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26  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 mixedinteger LP) solutions
and the objective function, as shown below.
1
2
3
4
5
6
7
8
9
10
11
12  # 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.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20  # 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.
1
2
3
4
5
6
7
8
9
10  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())
