Sunday, April 13, 2014

Playing with job shop problem ft10 (1)

As discussed in http://yetanothermathprogrammingconsultant.blogspot.com/2009/09/job-shop-scheduling.html a MIP solver like Gurobi can solve this difficult problem quite easily.

Here we show the performance with 1 and 4 threads:

Rplot01

The big advantage of using a MIP solver is that any time (after finding the first integer solution) we have both a best found (upper bound) and best possible or best bound (lower bound) available. This not only gives us a proven optimal solution at the end but a sense of the quality of the solution during the solution process. We can interrupt before the end using this information (e.g. by stopping on a small gap).

There are a few interesting observations:

  • The 4 thread integer solutions are not uniformly better than their 1 thread counterparts: quite often the red line (top part of the graph) is below the blue line
  • But in the end 4 threads wins big
  • Most progress is achieved in the beginning w.r.t. to finding good integer solutions. So when in a hurry we can stop early.
  • This also means that building a heuristic that finds a good solution quickly is of limited value: the MIP solver finds good solutions quickly on its own. I.e. such a heuristic helps where the MIP solver does not need much help.
  • The plot was produced with R and ggplot. I think it looks quite nice.
GAMS model

$ontext

   
Job Shop Scheduling

   
Manne formulation

   
Alan Manne, "On the job-shop scheduling problem", Operations Research
   
vol 8, no 2 March-April 1960.

   
Data set ft10 is from:
   
H. Fisher, G. L. Thompson, Probabilistic learning combinations of
   
local job-shop scheduling rules, in: J. F. Muth, G. L. Thompson (eds.),
   
Industrial Scheduling, Prentice Hall, Englewood Cliffs, N.J., 1963.

 
job1:   0 29 ; 1 78 ; 2  9 ; 3 36 ; 4 49 ; 5 11 ; 6 62 ; 7 56 ; 8 44 ; 9 21
 
job2:   0 43 ; 2 90 ; 4 75 ; 9 11 ; 3 69 ; 1 28 ; 6 46 ; 5 46 ; 7 72 ; 8 30
 
job3:   1 91 ; 0 85 ; 3 39 ; 2 74 ; 8 90 ; 5 10 ; 7 12 ; 6 89 ; 9 45 ; 4 33
 
job4:   1 81 ; 2 95 ; 0 71 ; 4 99 ; 6  9 ; 8 52 ; 7 85 ; 3 98 ; 9 22 ; 5 43
 
job5:   2 14 ; 0  6 ; 1 22 ; 5 61 ; 3 26 ; 4 69 ; 8 21 ; 7 49 ; 9 72 ; 6 53
 
job6:   2 84 ; 1  2 ; 5 52 ; 3 95 ; 8 48 ; 9 72 ; 0 47 ; 6 65 ; 4  6 ; 7 25
 
job7:   1 46 ; 0 37 ; 3 61 ; 2 13 ; 6 32 ; 5 21 ; 9 32 ; 8 89 ; 7 30 ; 4 55
 
job8:   2 31 ; 0 86 ; 1 46 ; 5 74 ; 4 32 ; 6 88 ; 8 19 ; 9 48 ; 7 36 ; 3 79
 
job9:   0 76 ; 1 69 ; 3 76 ; 5 51 ; 2 85 ; 9 11 ; 6 40 ; 7 89 ; 4 26 ; 8 74
 
job10:  1 85 ; 0 13 ; 2 61 ; 6  7 ; 8 64 ; 9 76 ; 5 47 ; 3 52 ; 4 90 ; 7 45


 
Description: jobs need to follow a sequence of operations. The numbers indicate
 
machine number and processing time.

 
This is a well-known benchmark problem and the optimal obj=930.

$offtext

set
  j
'tasks' / job1*job10/
  t
'stages'  /t1*t10/
;


table proctime(j,t) 'processing times for each stage'
       
t1    t2    t3    t4    t5    t6    t7    t8    t9   t10
job1    29    78     9    36    49    11    62    56    44   21
job2    43    90    75    11    69    28    46    46    72   30
job3    91    85    39    74    90    10    12    89    45   33
job4    81    95    71    99     9    52    85    98    22   43
job5    14     6    22    61    26    69    21    49    72   53
job6    84     2    52    95    48    72    47    65     6   25
job7    46    37    61    13    32    21    32    89    30   55
job8    31    86    46    74    32    88    19    48    36   79
job9    76    69    76    51    85    11    40    89    26   74
job10   85    13    61     7    64    76    47    52    90   45
;

table machine(j,t) 'machine numbers for each stage (zero based)'
       
t1  t2  t3  t4  t5  t6  t7  t8  t9  t10
job1    0    1   2   3   4   5   6   7   8   9
job2    0    2   4   9   3   1   6   5   7   8
job3    1    0   3   2   8   5   7   6   9   4
job4    1    2   0   4   6   8   7   3   9   5
job5    2    0   1   5   3   4   8   7   9   6
job6    2    1   5   3   8   9   0   6   4   7
job7    1    0   3   2   6   5   9   8   7   4
job8    2    0   1   5   4   6   8   9   7   3
job9    0    1   3   5   2   9   6   7   4   8
job10   1    0   2   6   8   9   5   3   4   7
;


alias (j,j2),(t,t2);

set ovl(j,t,j2,t2) 'this subset of (j,t),(j2,t2) needs to be checked against overlap'
;
ovl(j,t,j2,t2)$(
ord(j)<ord(j2) and machine(j,t)=machine(j2,t2)) = yes
;

variables

  x(j,t)       
'start time of subtask'
  y(j,t,j2,t2) 
'binary variable to implement non-overlap: (j,t) before or after (j2,t2)'
  z            
'objective variable'
;
binary variable y;
positive variable
x;


scalar TMAX 'max time horizon'
;
TMAX =
sum
((j,t), proctime(j,t));

equations

   NoOverlap1(j,t,j2,t2) 
'machine overlap'
   NoOverlap2(j,t,j2,t2) 
'machine overlap'
   Precedence(j,t)       
'orders need to follow a certain sequence of machines'
   zmax(j,t)             
'make span'
;


*
* each task within a job should start after previous one has finished
*
Precedence(j,t)$(
ord(t)<card(t)).. x(j,t+1) =G= x(j,t) + proctime(j,t);

*

* tasks on a machine cannot overlap:
* execute either before or after
*
NoOverlap1(ovl(j,t,j2,t2)).. x(j2,t2) =G= x(j,t) + proctime(j,t) - TMAX*y(j,t,j2,t2);
NoOverlap2(ovl(j,t,j2,t2)).. x(j,t) =G= x(j2,t2) + proctime(j2,t2) - TMAX*(1-y(j,t,j2,t2));


*
* minimize completion time of last task
*
zmax(j,t).. z =G= x(j,t) + proctime(j,t);

option optcr=0;
model jobshop /all/
;
solve
jobshop minimizing z using mip;

R Code to create the plot


library(ggplot2)   
setwd("C:/projects/tmp")
# read csv files
d1<-read.csv("1thread.csv")
d4<-read.csv("4thread.csv")
# stack them into a single data frame
d<-rbind(d1,d4)
# add a Threads column indicating the data set
d$Threads <- c(rep("1 thread",nrow(d1)),rep("4 threads",nrow(d4)))
# convert factor to double ("na" will be converted to NA)
d$bestFound <- as.double(as.character(d$bestFound))
# actual plotting
ggplot(data=d)+
geom_line(aes(group=Threads,x=seconds,y=bestBound,color=Threads))+
geom_line(aes(group=Threads,x=seconds,y=bestFound,color=Threads))+
ggtitle("Gurobi laptop performance on ft10.gms")+
xlab("Time (seconds)")+
ylab("Best Found (top), Best Bound (bottom)")




Created by Pretty R at inside-R.org