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

Question

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:

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

subject to

\begin{equation} \sum_{j=1}^{n}{z_{j,i}} = 1, \quad 1 \leqslant i \leqslant n \end{equation}

\begin{equation} \sum_{i=1}^{n}{z_{j,i}} = 1, \quad 1 \leqslant j \leqslant n \end{equation}

\begin{equation} s_{1,1}=0 \end{equation}

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

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

\begin{equation} 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 \end{equation}

\begin{equation} 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 \end{equation}

\begin{equation} z_{j,i} \in \{0,1\}, \quad 1 \leqslant j \leqslant n, \quad 1 \leqslant i \leqslant n \end{equation}

\begin{equation} s_{r,i} \geqslant 0, \quad 1 \leqslant r \leqslant m, \quad 1 \leqslant i \leqslant n \end{equation}

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.

Answer 1

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:

NumofMachines = Tasktime.shape[0]
NumofJobs = Tasktime.shape[1]

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


# Generating variables:
z = {}
for j in range(NumofJobs):
    for i in range(NumofJobs):
        z[j, i] = m.addVar(vtype=GRB.BINARY)

s = {}
for r in range(NumofMachines):
    for j in range(NumofJobs):
        s[r, j] = m.addVar(vtype=GRB.CONTINUOUS)

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)

m.addConstr(s[0, 0] == 0)

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()

Answer 2

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 = buildModel(NumofJobs, NumofMachines, Tasktime)

    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:
    z = m.addVars(NumofJobs, NumofJobs, vtype=grb.GRB.BINARY)

    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)

    m.addConstr(s[0, 0] == 0)

    for i in range(NumofJobs - 1):
        m.addConstr(s[0, i] + grb.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] + grb.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] + grb.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] + grb.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] + 
        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!

Answer 3

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.

Answer 4

@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.