# How to write a mixed-integer linear programming formulation in Python using Gurobi?

I can remember that I spent some time in understanding how to formulate my first model. So I aimed at presenting a complete model here, wishing to save some time for students or researchers needing it.

The model is a flow shop scheduling problem, presented in Wilson (1989), as following:

$$$$\label{eq1} z = \min (s_{m,n} + \sum_{j=1}^{n}{p_{m,j} z_{j,n}})$$$$

subject to

$$$$\sum_{j=1}^{n}{z_{j,i}} = 1, \quad 1 \leqslant i \leqslant n$$$$

$$$$\sum_{i=1}^{n}{z_{j,i}} = 1, \quad 1 \leqslant j \leqslant n$$$$

$$$$s_{1,1}=0$$$$

$$$$s_{1,i} + \sum_{j=1}^{n}{p_{1,j} z_{j,i} = s_{1,i+1}}, \quad 1 \leqslant i \leqslant n-1$$$$

$$$$s_{r,1} + \sum_{j=1}^{n}{p_{r,j} z_{j,1} = s_{r+1,1}}, \quad 1 \leqslant r \leqslant m-1$$$$

$$$$s_{r,i} + \sum_{j=1}^{n}{p_{r,j} z_{j,i}\leqslant s_{r+1,i}}, \quad 1 \leqslant r \leqslant m-1, \quad 2 \leqslant i \leqslant n$$$$

$$$$s_{r,i} + \sum_{j=1}^{n}{p_{r,j} z_{j,i}\leqslant s_{r,i+1}}, \quad 2 \leqslant r \leqslant m, \quad 1 \leqslant i \leqslant n-1$$$$

$$$$z_{j,i} \in \{0,1\}, \quad 1 \leqslant j \leqslant n, \quad 1 \leqslant i \leqslant n$$$$

$$$$s_{r,i} \geqslant 0, \quad 1 \leqslant r \leqslant m, \quad 1 \leqslant i \leqslant n$$$$

Note that $$s_{r,i}$$ is the starting time of job in position $$i$$ on machine $$r$$, and $$z_{j,i}$$ is equal to 1 if job $$j$$ is assigned to position $$i$$. Also, $$p_{r,j}$$ is the processing time of job $$j$$ on machine $$r$$. I don't go to the details of the model as in not the purpose of this post.

So, the question is how to formulate this model in Python, using the Gurobi solver. i.e. using the module gurobipy.

Details of model can be found in: Wilson JM. Alternative formulations of a flow-shop scheduling problem. Journal of the Operational Research Society (1989) 40:395–399.

Here is the complete implementation for the above-mentioned model.

from gurobipy import *
import numpy as np
import csv

# Parameters needed are:
(1) the total number of jobs (n). Here I denote it by "NumofJobs"
(2) the total number of machines (m). Here I denote it by "NumofMachines"
(3) the processing times.  Here I use a numpy matrix: "Tasktime[r, j]" : p_{r,j}

# Building the model:
m = Model ("Wilson")

# Generating variables:
z = {}
for j in range(NumofJobs):
for i in range(NumofJobs):

s = {}
for r in range(NumofMachines):
for j in range(NumofJobs):

m.update()

# Generating constraints:
for i in range(NumofJobs):
m.addConstr(quicksum(z[j, i] for j in range(NumofJobs)) == 1)

for j in range(NumofJobs):
m.addConstr(quicksum(z[j, i] for i in range(NumofJobs)) == 1)

for i in range(NumofJobs - 1):
m.addConstr(s[0, i] + quicksum((Tasktime[0, j]) * z[j, i] for j in range(NumofJobs)) == s[0, i + 1])

for r in range(NumofMachines - 1):
m.addConstr(s[r, 0] + quicksum((Tasktime[r, j]) * z[j, 0] for j in range(NumofJobs)) == s[r + 1, 0])

for r in range(NumofMachines - 1):
for i in range(1, NumofJobs):
m.addConstr(s[r, i] + quicksum(Tasktime[r, j] * z[j, i] for j in range(NumofJobs)) <= s[r + 1, i])

for r in range(1, NumofMachines):
for i in range(NumofJobs - 1):
m.addConstr(s[r, i] + quicksum(Tasktime[r, j] * z[j, i] for j in range(NumofJobs)) <= s[r, i + 1])

# Setting the objective function:
m.setObjective(s[NumofMachines - 1, NumofJobs - 1] + quicksum(Tasktime[NumofMachines - 1, j] * z[j, NumofJobs - 1] for j in range(NumofJobs)), GRB.MINIMIZE)

m.optimize()

• My unsolicited advice: I would recommend changing the import statement from from gurobipy import * to import gurobipy as grb (this is better design imo and also helps when you have bigger models, multiple files, etc.). Also, variables can be added more efficiently without for loops: z=m.addVars(NumofJobs, NumofJobs, vtype=grb.GRB.BINARY). – David M. Jun 27 '19 at 13:15
• Also, lines 6-8 should be commented out, and I don’t think you use the csv module, so no need to import it! – David M. Jun 27 '19 at 13:21
• Instead of using the Gurobi solver, could the OP have also just used any other solver like CONOPT or CPLEX? Just curious from someone new to the field. – GrayLiterature Jun 27 '19 at 13:22
• Since you already wrote the solution, it'd be great to show how you retrieve the objective function and variables values too. Also, after Gurobi 6.5 or 7.0, you don't need to call model update() in most programs (including your small example here) – EhsanK Jun 27 '19 at 13:40
• @D.Gray I am mostly using Gurobi for my own works. But I am planning to do the same for other solvers and post them as well. BTW I think Gurobi is the easiest to implement, and a powerful solver. – Mostafa Jun 28 '19 at 11:04

Substantively similar to the OP's answer, but with some Python tweaks. Mostly, I just put things in their own functions, tweaked some Gurobi API calls to be cleaner and more efficient, and provided an example of how to check solve output after you solve.

Part I: Eliminated extra imports + added comment characters

import gurobipy as grb

# Parameters needed are:
# (1) the total number of jobs (n). Here I denote it by "NumofJobs"
# (2) the total number of machines (m). Here I denote it by "NumofMachines"
# (3) the processing times.  Here I use a numpy matrix: "Tasktime[r, j]" : p_{r,j}


Part II: Wrote a separate main function to separate the building from the solving.

def main(): # Always wrap things in functions!
# Separated the building of the model from the solving of the model.
NumofJobs = 10
NumofMachines = 14
Tasktime = ... # Whatever data you want

model.optimize()
# Check to see what happened!
if (model.status != grb.GRB.OPTIMAL):
print('Oh no, something went wrong!')


Part III: Build the actual model.

def buildModel(NumofJobs, NumofMachines, Tasktime)
# Building the model:
m = grb.Model("Wilson")

# Generating variables:

s = m.addVars(NumofMachines, NumofMachines) # Default is continuous

# Don't need to call m.update() these days.

# Generating constraints:
for i in range(NumofJobs):
m.addConstr(grb.quicksum(z[j, i] for j in range(NumofJobs)) == 1)

for j in range(NumofJobs):
m.addConstr(grb.quicksum(z[j, i] for i in range(NumofJobs)) == 1)

for i in range(NumofJobs - 1):
for j in range(NumofJobs)) == s[0, i + 1])

for r in range(NumofMachines - 1):
for j in range(NumofJobs)) == s[r + 1, 0])

for r in range(NumofMachines - 1):
for i in range(1, NumofJobs):
for j in range(NumofJobs)) <= s[r + 1, i])

for r in range(1, NumofMachines):
for i in range(NumofJobs - 1):
for j in range(NumofJobs)) <= s[r, i + 1])

# Setting the objective function:
m.setObjective(s[NumofMachines - 1, NumofJobs - 1] +
grb.quicksum(Tasktime[NumofMachines - 1, j] * z[j, NumofJobs - 1]
for j in range(NumofJobs)), sense=grb.GRB.MINIMIZE)

return m

if __name__=='__main__': # Standard way to call Python main function
main()


Of course there are lots more things that could be done here, but this is a template of how I write prototypes. As has been discussed, the OR community can always do better with our code!

To complete the answer, I will add some other things that may be useful.

1. To change the parameters of the solver gurobi, e.g. setting the time-limit to 600 seconds: m.setParam('TimeLimit', 600)

2. To retrieve the objective function of the problem: Objective_Of_The_Problem = m.objVal

3. To retrieve the status of the problem: Status_Of_The_Problem = m.status

4. To retrieve the run-time of the solver: Time_To_Solve = m.Runtime

5. To retrieve the value of a variable $$z[j, i]$$: Wilson_Variable = z[j, i].x

Note that we defined m = Model ("Wilson") previously.