Sunday, August 19, 2012

Problems reproducing results

In quite a few cases I have been unable to reproduce published results. Here is an example. (Update: actually the presented results are correct, just a typo in the problem description: see the comments). In the paper http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.138.5011 an example is given of a job shop problem. The data is as follows:

image

I.e. there are 10 jobs each having 3 job steps. The job steps have to be executed in a certain given order (it might have been easier just to reorder the job step numbers and say all job steps have to be executed in the given order). The authors mention the following:

image

However when I try to reproduce this with a MIP model I can achieve a makespan of 7. Below is the GAMS model:

$ontext

  
We try to reproduce the results from

  
KHALED MESGHOUNI, SLIM HAMMADI, PIERRE BORNE
  
EVOLUTIONARY ALGORITHMS FOR JOB-SHOP SCHEDULING
  
Int. J. Appl. Math. Comput. Sci., 2004, Vol. 14, No. 1, 91–103

  
erwin@amsterdamoptimization.com


$offtext


sets
   j 
'jobs' /j1*j10/
   m 
'machines' /m1*m10/
   op
'operations' /op1*op3/
   js
'job steps' /js1*js30/
;


table ordering(j,op)
   
Op1 Op2 Op3

J1   1   3   2
J2   2   1   3
J3   1   2   3
J4   1   2   3
J5   2   3   1
J6   1   2   3
J7   3   2   1
J8   3   1   2
J9   1   2   3
J10  3   2   1
;


table proctime(j,op,m) 'processing times'

       
M1    M2    M3    M4    M5    M6    M7    M8    M9   M10
J1.Op1   1     4     6     9     3     5     2     8     9     5
J1.Op2   3     2     5     1     5     6     9     5    10     3
J1.Op3   4     1     1     3     4     8    10     4    11     4
J2.Op1   4     8     7     1     9     6     1    10     7     1
J2.Op2   2    10     4     5     9     8     4    15     8     4
J2.Op3   6    11     2     7     5     3     5    14     9     2
J3.Op1   8     5     8     9     4     3     5     3     8     1
J3.Op2   9     3     6     1     2     6     4     1     7     2
J3.Op3   7     1     8     5     4     9     1     2     3     4
J4.Op1   5    10     6     4     9     5     1     7     1     6
J4.Op2   4     2     3     8     7     4     6     9     8     4
J4.Op3   7     3    12     1     6     5     8     3     5     2
J5.Op1   6     1     4     1    10     4     3    11    13     9
J5.Op2   7    10     4     5     6     3     5    15     2     6
J5.Op3   5     6     3     9     8     2     8     6     1     7
J6.Op1   8     9    10     8     4     2     7     8     3    10
J6.Op2   7     3    12     5     4     3     6     9     2    15
J6.Op3   4     7     3     6     3     4     1     5     1    11
J7.Op1   5     4     2     1     2     1     8    14     5     7
J7.Op2   3     8     1     2     3     6    11     2    13     3
J7.Op3   1     7     8     3     4     9     4    13    10     7
J8.Op1   8     3    10     7     5    13     4     6     8     4
J8.Op2   6     2    13     5     4     3     5     7     9     5
J8.Op3   5     7    11     3     2     9     8     5    12     8
J9.Op1   3     9     1     3     8     1     6     7     5     4
J9.Op2   4     6     2     5     7     3     1     9     6     7
J9.Op3   8     5     4     8     6     1     2     3    10    12
J10.Op1  9     2     4     2     3     5     2     4    10    23
J10.Op2  3     1     8     1     9     4     1     4    17    15
J10.Op3  4     3     1     6     7     1     2     6    20     6
;

*
* set up a mapping between the job steps (js) and the job/operations (j,o)
*
sets
   jsmap(js,j,op)
'job step mapping between js and (j,o)'
   jscur(js)
'temp' /js1/
;

loop((j,op),
  jsmap(jscur,j,op)=
yes
;
  jscur(js)=jscur(js-1);
);

display
jsmap;

*

* convert data to jobsteps
*
parameter proctime2(js,m) 'processing time';
proctime2(js,m) =
sum
(jsmap(js,j,op),proctime(j,op,m));
display
proctime2;

alias
(op,op2),(js,js2),(j,j2);

set prec1(j,op,op2) 'ordering as precedence relations'
;
scalar
s;
for(s=1 to card
(op)-1,
   prec1(j,op,op2)$(ordering(j,op)=s
and ordering(j,op2)=s+1)=yes
;
);

display
prec1;

set prec2(js,js2) 'ordering in terms of job steps'
;
prec2(js,js2) =
sum((j,op,op2)$(jsmap(js,j,op) and
jsmap(js2,j,op2)
                 
and
prec1(j,op,op2)),1);
display
prec2;



variables

  start(js)          
'start time of task'
  finish(js)         
'completion time'
  assign(js,m)       
'assign job step to machine'
  overlap(js,js2,m)  
'binary variable to implement non-overlap'
  z                  
'objective variable'
;
binary variable overlap,assign;
positive variable
start,finish;


scalar TMAX 'max time horizon'
;
TMAX =
sum(js,smin
(m,proctime2(js,m)));

equations

   onemach(js)
   calcend
   NoOverlap1(js,js2,m)  
'machine occupation'
   NoOverlap2(js,js2,m)  
'machine occupation'
   Precedence(js,js2) 
'orders need to follow a certain sequence of steps'
   zmax(js)          
'make span'
;


OneMach(js)..
sum(m, assign(js,m)) =e= 1;

calcEnd(js).. finish(js) =e= start(js) +
sum
(m, assign(js,m)*proctime2(js,m));

Precedence(prec2(js,js2)).. start(js2) =g= finish(js);


set ut(js,js2) 'upper triangular part'
;
ut(js,js2)$(
ord(js)<ord(js2))=yes
;

NoOverlap1(ut(js,js2),m)..
   start(js) =g= finish(js2) - TMAX*overlap(js,js2,m) - TMAX*(1-assign(js,m)) -TMAX*(1-assign(js2,m));
NoOverlap2(ut(js,js2),m)..
   start(js2) =g= finish(js)  - TMAX*(1-overlap(js,js2,m))- TMAX*(1-assign(js,m)) -TMAX*(1-assign(js2,m));

zmax(js).. z =G= finish(js);


option
optcr=0;
model jobshop /all/
;
solve
jobshop minimizing z using mip;
display
start.l,finish.l,assign.l;

parameters

   start2(j,op)
   finish2(j,op)
   assign2(j,op,m)
;
start2(j,op) =
sum(jsmap(js,j,op),start.l(js));
finish2(j,op) =
sum
(jsmap(js,j,op),finish.l(js));
assign2(j,op,m) =
sum
(jsmap(js,j,op),assign.l(js,m));
display
start2,finish2,assign2;

parameter
results(j,op,*);
results(j,op,
'start'
) = start2(j,op);
results(j,op,
'proc') = sum
(m,assign2(j,op,m)*proctime(j,op,m));
results(j,op,
'finish'
) = finish2(j,op);
results(j,op,
'machine') = sum(m,assign2(j,op,m)*ord
(m));
option
results:0;
display
results;


The model has a few interesting wrinkles:

  • I remap the job steps from tuples (job1,op1),(job1,op2),.. to a one dimensional set js1,js2,… This causes some headaches when converting the input data but makes the model equations simpler.
  • I map the solution back to the tuples (job1,op1),(job1,op2),.. after the solve.
  • A problem with these models is always to find a good big-M. For real problems you may find one by preceding the MIP model with a simple heuristic that constructs an initial solution.
  • I try to make sure we don’t compare the job steps twice for overlap if they are scheduled on the same machine. Hence the set UT.

Here are my results:

----    157 PARAMETER results 

              start        proc      finish     machine

j1 .op1           1           1           2           1
j1 .op2           6           1           7           4
j1 .op3           2           1           3           2
j2 .op1           4           1           5           4
j2 .op2           2           2           4           1
j2 .op3           5           2           7           3
j3 .op1                       1           1          10
j3 .op2           3           1           4           4
j3 .op3           6           1           7           7
j4 .op1                       1           1           7
j4 .op2           1           3           4           3
j4 .op3           4           2           6          10
j5 .op1           2           1           3           4
j5 .op2           4           2           6           9
j5 .op3                       1           1           9
j6 .op1                       2           2           6
j6 .op2           2           2           4           9
j6 .op3           4           1           5           7
j7 .op1           5           1           6           4
j7 .op2           1           2           3           8
j7 .op3                       1           1           1
j8 .op1           4           3           7           2
j8 .op2                       2           2           2
j8 .op3           2           2           4           5
j9 .op1           4           1           5           3
j9 .op2           5           1           6           7
j9 .op3           6           1           7           6
j10.op1           4           3           7           5
j10.op2           1           1           2           4
j10.op3                       1           1           3

 

image

Of course it is possible there is a transcription error in the data table, or may be I am misreading the model, or some other trivial problem. However the results look reasonable when looking at the GANTT charts, i.e. there is no overlap, and the selected machines for processing the job steps are efficient in doing so (we did not select slow machines).

This is actually a fairly small problem and a good MIP code has no problem with it. For instance Gurobi says:

Explored 0 nodes (287 simplex iterations) in 0.47 seconds

I would expect that genetic algorithms become more interesting for larger instances where a MIP solver may not get anywhere.