$title Box Packing Problem (boxpacking,SEQ=434)

$onText
This model tries to find the optimal configuration of boxes in a container
given the 3-dimensional measures of both objects while maintaining the
geometrical restrictions. The model uses different methods to determine the order
in which the boxes should be packed. GAMS Connect is used to read the data from the br1.csv file.
The data is from Bischoff/Ratcliff paper on container packing
and is available at the ORLIB instance collection, http://people.brunel.ac.uk/~mastjjb/jeb/orlib/thpackinfo.html

This formulation is described in detail in:
Ocloo, Valentina & Fügenschuh, Armin & Pamen, Olivier. (2020).
A New Mathematical Model for a 3D Container Packing Problem. 10.26127/btuopen-5088. 

Keywords: mixed integer linear programming, container packing problem,
          3D-bin packing problem, box packing problem
$offText

$eolCom !!

$if not set DEBUG        $set DEBUG         0        !! DEBUG=1 activates additional output
$if not set METHOD       $set METHOD        batch    !! [standard, greedy, batch]
$if not set FIXLOCDIM    $set FIXLOCDIM     0        !! FIXLOCDIM=1 fixes location of placed boxes in the container
$if not set MAXCOPIES    $set MAXCOPIES     50       !! maximum number of copies per box type
$if not set BATCH_SIZE   $set BATCH_SIZE    5        !! number of boxes per batch
$if not set TIMELIMIT    $set TIMELIMIT     30       !! total time limit

* define the 3D measure of the container
$if not set X_CONT $set X_CONT 587
$if not set Y_CONT $set Y_CONT 233
$if not set Z_CONT $set Z_CONT 220

* raw input data
Set box_id      'input data box ids'
    box_id_hdr  'input data header'   / x, rot_x, y, rot_y, z, rot_z, nb, weight/
    mult_box    'input data box copy id index for multiple boxes of the same type' / cp1*cp%MAXCOPIES% /;

Parameter
    box_data(box_id<, box_id_hdr)
    dim_raw(box_id,box_id_hdr);

* Reading the data
$onEmbeddedCode Connect:
- CSVReader:
    file: br1.csv
    name: box_data
    header: True
    indexColumns: 1
    valueColumns: "2:lastCol"
- GAMSWriter:
    symbols: all
$offEmbeddedCode

* reduce the number of boxes,
* the full scale problem becomes computationally challenging
box_data(box_id, 'nb') = ceil(box_data(box_id, 'nb')/2.5);

display$%DEBUG% box_data

$eval MAXB card(box_id)*%MAXCOPIES%

* raw data consistency checks
abort$[smax(box_id, box_data(box_id, 'nb')) > %MAXCOPIES%] 'Input data contains more copies of a box than supported. Increase MAXCOPIES';

* internal data structures
Set bb                        super set of boxes       / b1*b%MAXB%, vbox 'virtual box' /
    b_exist(bb)               set of existing boxes
    b(bb)                     dynamic subset of boxes
    bmap(box_id,mult_box,bb)  mapping of input data box ids with copy id to internal box ids
    bmap_red(bb,box_id)       mapping of input data box ids to internal box ids    
    d(box_id_hdr)             dimensions                            / x, y, z /
    r(box_id_hdr)             possibility of rotating a box         / rot_x, rot_y, rot_z /    
    o                         six possible orientations of a box    / o1*o6 /                         
    rmap(o,d,d)               rotations that define mapping of input dimensions to orientation dimensions
                            / o1.(x.x,y.y,z.z) 'input                                   x->x y->y z->z'                                              
                              o2.(x.x,y.z,z.y) 'rotate input 90 degrees along x axis    x->x y->z z->y (not possible if rot_x = 0)'               
                              o3.(x.z,y.y,z.x) 'rotate input 90 degrees along y axis    x->z y->y z->x (not possible if rot_y = 0)'               
                              o4.(x.y,y.x,z.z) 'rotate input 90 degrees along z axis    x->y y->x z->z (not possible if rot_z = 0)'               
                              o5.(x.z,y.x,z.y) 'rotate input 90 degrees along x+z axis  x->z y->x z->y (not possible if rot_x = 0 or rot_z = 0)'  
                              o6.(x.y,y.z,z.x) 'rotate input 90 degrees along x+y axis  x->y y->z z->x (not possible if rot_x = 0 or rot_y = 0)' /
    ro(r,o)                   rotation restriction  / rot_x.(o2,o5,o6), rot_y.(o3,o6), rot_z.(o4,o5) /
    bo(bb,o)                  available orientations   / #bb.#o /;

Alias (bb,bb1,bb2), (b,b1,b2), (d,d1);

Parameter
    dim_cont(d)   'container dimesnions'                    / x %X_CONT%, y %Y_CONT%, z %Z_CONT% /
    dim_o(bb,o,d) 'dimensions of box bb for orientation o'
    rot(bb, r)    'allowed rotations'
    Box_vol(bb)   'volume of each box';

* map raw box ids to internal box ids
b('b1') = yes;
loop((box_id,mult_box)$[ord(mult_box) <= box_data(box_id,'nb')],
  bmap(box_id,mult_box,b) = yes;
  b(bb) = b(bb-1);
);
option clear=b;
option bmap_red<bmap, b_exist<bmap;

* compute box dimensions for all rotations
dim_o(bb,o,d) = sum((bmap_red(bb,box_id),rmap(o,d1,d)), box_data(box_id,d1));

* compute available rotations 
rot(bb, r)    = sum(bmap_red(bb,box_id), box_data(box_id, r));

* compute volume for each box, scaled to m3
Box_vol(bb)   = prod(d, dim_o(bb,'o1', d)/100);

* disable the impossible orientations of a box given its rotation restriction
loop(r, bo(bb,o)$(ro(r,o) and not rot(bb,r)) = no);

Binary Variable
    OMEGA(bb)       'box in container (1) or not (0)'
    ALPHA(bb,o)     'box in container with orientation o (1) or not (0)'
    RELPOS(bb,bb,d) 'relative position: 1 if first box location + dimension <= second box location';

NonNegative Variable
    LOC(bb,d)    '(x,y,z) location of bottom left back corner of bb in container'
    DIM(bb,d)    '(x,y,z) dimension of box bb in container considering its orientatin';

Variable
    VOL          'Objective value - total volume of all boxes in container';

Equation
    eq_def_DIM(bb,d)              'define dimension of box in container considering orientation'
    eq_couple_ALPHA_OMEGA(bb)     'select orientation only for boxes in container' 
    eq_inside_container(bb,d)     'respect container dimensions'
    eq_deactivate_RELPOS(bb,bb,d) 'RELPOS should only take non-zero value if first box in container'
    eq_no_overlap(bb,bb, d)       'boxes in container must not overlap' 
    eq_def_RELPOS(bb,bb)          'define RELPOS, i.e. ist first box left/behind/below second box (1), otherwise (0)'
    eq_def_VOL                    'Objective function - maximize the volume utilization of container';

* (23/24/25) [The numbers correspond to the equation number in the paper]
eq_def_DIM(b,d)$[not sameas(b,'vbox')]..  DIM(b,d) =e= sum(o, dim_o(b,o,d)*ALPHA(b,o)); 

* (26)
eq_couple_ALPHA_OMEGA(b)$[not sameas(b,'vbox')]..  sum(bo(b,o), ALPHA(b,o)) =e= OMEGA(b);

* (27 a/b/c)
eq_inside_container(b,d)..  LOC(b,d) + DIM(b,d) =l= dim_cont(d) * OMEGA(b);

* (28 a/b/c)
eq_deactivate_RELPOS(b1(bb1),b2(bb2),d)$[ord(bb1)<ord(bb2)]..  RELPOS(b1,b2,d) + RELPOS(b2,b1,d) =l= OMEGA(b1);

* (29)
eq_def_RELPOS(b1(bb1),b2(bb2))$[ord(bb1)<ord(bb2)]..  sum(d, RELPOS(b1,b2,d) + RELPOS(b2,b1,d)) =g= OMEGA(b1) + OMEGA(b2) - 1;

* (30 a/b/c)
eq_no_overlap(b1,b2,d)$[not sameas(b1,b2)]..   LOC(b1,d) + DIM(b1,d) =l= LOC(b2,d) + dim_cont(d)*(1-RELPOS(b1,b2,d));

* (31)
eq_def_VOL.. VOL =e= sum(b, Box_vol(b)*OMEGA(b));

Model cpp / all /;

cpp.reslim = %TIMELIMIT%;

$ifThenI %system.mip%==cplex
* create Cplex option file
    $$onEcho > cplex.opt
    solvefinal 0
    mipstart   1
    $$offEcho
* use option file
    cpp.optfile = 1;    
$elseIfI %system.mip%==gurobi
* create gurobi option file
    $$onEcho > gurobi.opt
    solvefixed      0
    mipstart        1
    norelheurtime   %TIMELIMIT%
    $$offEcho
* use option file
    cpp.optfile = 1;    
$endIf

Scalar num_of_bat;
num_of_bat = ceil(card(b_exist)/%BATCH_SIZE%);

option limrow=0, limcol=0, solprint=off, solvelink=5;
$ifE %DEBUG%==1 option limrow=1e6, limcol=1e6, solprint=on;

* solve for all existing boxes at once
$ifThenI.method %METHOD%==standard
    b(b_exist) = yes;
    solve cpp max VOL use MIP;

* iteratively solve for each box
$elseIfI.method %METHOD%==greedy
loop(bb$b_exist(bb),
    b(bb) = yes;
    solve cpp max VOL use MIP;
    OMEGA.fx(bb)   = round(OMEGA.l(bb));
    $$ifE %FIXLOCDIM%=1 ALPHA.fx(bb,o) = round(ALPHA.l(bb,o)); LOC.fx(bb,d) = LOC.l(bb,d); DIM.fx(bb,d) = DIM.l(bb,d);    
  );

* put n boxes in a batch and solve model
$elseIfI.method %METHOD%==batch
Scalar i;
for(i = 1 to num_of_bat,
    b(bb)$(ord(bb) <= %BATCH_SIZE%*i) = yes;
    solve cpp max VOL use MIP;
    OMEGA.fx(b)   = round(OMEGA.l(b));
    $$ifE %FIXLOCDIM%=1 ALPHA.fx(b,o) = round(ALPHA.l(b,o)); LOC.fx(b,d) = LOC.l(b,d); DIM.fx(b,d) = DIM.l(b,d);  
);
$endIf.method

Parameter rep(*) report;
rep('total number of boxes')                       = card(b_exist);
rep('number of boxes in container')                = sum(b_exist, OMEGA.l(b_exist));
rep('volume of container [m3]')                    = prod(d, dim_cont(d)/100);
rep('volume smallest box [m3]')                    = smin(b_exist, Box_vol(b_exist));
rep('volume biggest box [m3]')                     = smax(b_exist, Box_vol(b_exist));
rep('total volume of all boxes [m3]')              = sum(b_exist, Box_vol(b_exist));
rep('total volume of boxes in container [m3]')     = sum(b_exist, Box_vol(b_exist)*OMEGA.L(b_exist));
rep('ratio of boxes in container')                 = (rep('number of boxes in container') / rep('total number of boxes'));
rep('ratio of total volume of boxes in container') = (rep('total volume of boxes in container [m3]') / rep('total volume of all boxes [m3]'));
rep('ratio of container volume utilized')          = (rep('total volume of boxes in container [m3]') / rep('volume of container [m3]'));
option rep:2:0:1; display rep;