# How to model 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

# 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}

# Based on the matrix Tasktime, we can set:

# 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). Jun 27, 2019 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! Jun 27, 2019 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. Jun 27, 2019 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, 2019 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. Jun 28, 2019 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!

• What kind of data is used for Tasktime = ... # Whatever data you want? Is it a 10x14 matrix? Can I get an example? Mar 15, 2021 at 10:35
• @JonasMård Apologies for mishandling your case earlier today. I have converted your post into a comment, which I believe is for David (hopefully you get an answer). As the comment threshold is 50 reputation, you won't be able to post further comments, but you are welcome to ask a question or browse the list of unanswered questions to gain the rep needed Mar 15, 2021 at 21:54
• @JonasMård Since the model is defined with $Tasktime[r, j] : p_{r,j}$, where $r$ and $j$ are the indices for machines and jobs, the size of the matrix should be 14×10, if NumofJobs = 10 and NumofMachines = 14. Mar 16, 2021 at 5:20

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.

@Mostafa, j is the number of jobs while r is the number of machines, this means $$NumofJobs = Tasktime.shape[1]$$ and $$NumofMachines = Tasktime.shape[0]$$, and also, s must be $$s = m.addVars(NumofMachines,NumofJobs)$$, otherwise the optimality cannot be achieved. Do this modification and run the example by using the data in Table 1 found in your reference paper(Wright, 1989), you will easily see the optimal cost of 83.

• Nope, $r$ and $j$ are indices for machines and jobs, the total number of them are $m$ and $n$ in the model, and 'NumofMachines', 'NumofJobs' in the implementation. Sep 21, 2021 at 21:23
• Yes, I missed it, anyway, the code is correct in which shape[1] represents j whereas shape[0] refers to i in i x j flow data, which automates the process if you use i as the machine number and j as the job number.
– YcK
Sep 22, 2021 at 8:00
• Yep, I will add your idea to the implementation, Yck. Thanks for the suggestion. Sep 22, 2021 at 8:06