14
$\begingroup$

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.

$\endgroup$
12
$\begingroup$

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()
$\endgroup$
8
  • 6
    $\begingroup$ 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). $\endgroup$
    – David M.
    Jun 27 '19 at 13:15
  • 1
    $\begingroup$ Also, lines 6-8 should be commented out, and I don’t think you use the csv module, so no need to import it! $\endgroup$
    – David M.
    Jun 27 '19 at 13:21
  • 2
    $\begingroup$ 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. $\endgroup$ Jun 27 '19 at 13:22
  • 2
    $\begingroup$ 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) $\endgroup$
    – EhsanK
    Jun 27 '19 at 13:40
  • 2
    $\begingroup$ @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. $\endgroup$
    – Mostafa
    Jun 28 '19 at 11:04
6
$\begingroup$

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.

$\endgroup$
5
$\begingroup$

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!

$\endgroup$
3
  • $\begingroup$ What kind of data is used for Tasktime = ... # Whatever data you want? Is it a 10x14 matrix? Can I get an example? $\endgroup$ Mar 15 at 10:35
  • $\begingroup$ @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 $\endgroup$
    – TheSimpliFire
    Mar 15 at 21:54
  • $\begingroup$ @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. $\endgroup$
    – Mostafa
    Mar 16 at 5:20
1
$\begingroup$

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

New contributor
YcK is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Mostafa
    2 days ago
  • $\begingroup$ 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. $\endgroup$
    – YcK
    2 days ago
  • $\begingroup$ Yep, I will add your idea to the implementation, Yck. Thanks for the suggestion. $\endgroup$
    – Mostafa
    2 days ago

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.