%% A MIP-based slicing heuristic for three-dimensional bin packing
% ECE 602 Final project
% 
% Yinghui Zhao 20742812
%
% 26/04/2019
%
% Professor: Oleg Michailovich

%% Introduction and problem formulation
% The mixed-case palletization problem describes the situation when a set of 
% boxes of rectangular shape with different sizes and stability have to be packed 
% into pallets of given dimensions for ease of transportation, storage, and 
% display at stores. The problem, in essence, a three-dimensional bin 
% packing problem(3DBPP) with side constraints. In the paper, the focus is 
% on the modelling and solution of the open dimension bin packing problem, 
% which refers to the particular case where only one bin is used, and the 
% objective is to minimize the final height of the packed boxes.
%
% A mix integer programming formulation, which maximizes the density of the
% bottom layers and the compactness of the pallet to ensure the stability of
% the top layers, is proposed. A relative-position concept, where the 
% relative positions of boxes are the main decision variables, is used. 
% 
%% Proposed solution
% A relative-position formulation [R P] is firstly proposed. The formulation is:
%
% <html>
% <ul>
% <li>
% <img
% src="http://latex.codecogs.com/gif.latex?c_{i}^{s}:%20\quad%20\text{the%
% 20front,%20lower,%20left%20coordinate%20of%20box%20i%20along%20direction%
% 20s%20=%20x,%20y,%20z}" height="12">
% </li>
% <li>
% <img
% src="http://latex.codecogs.com/gif.latex?$$P_{ij}^{s}%20=%20\begin{cases}%201,%20\quad%20\text{if%20box%20i%20precedes%20box%20j%20along%20direction%20s.}\\%200,%20\quad%20\text{otherwise.}%20\end{cases}$$"
% height="36">
% </li>
% </ul>
% </html>
%
% $$[R\quad P]: min \quad h $$
%
% $$s.t. \sum_{s=x,y,z}{(P_{ij}^{s} + P_{ji}^{s}) \geq 1} \qquad \forall j > i \qquad (1)$$
%
% $$P_{ij}^{s} + P_{ji}^{s} \leq 1 \qquad \forall j > i,s \qquad (2)$$
%
% $$c_{i}^{s} + l_{i}^{s} \leq c_{j}^{s} + L^{s}(1 - p_{ij}^{s}) \qquad i \neq j,s \qquad (3)$$
%
% $$0 \leq c_{i}^{s} \leq L^{s} - l_{i}^{s} \qquad \forall i,s \qquad (4)$$
%
% $$c_{i}^{z} + l_{i}^{z} \leq h \qquad \forall i \qquad (5)$$
%
% $$p_{ji}^{s} \in \{0, 1\} \qquad \forall i \neq j,s \qquad (6)$$
%
% $$h \geq 0 \qquad (7)$$
%
% However, the authors of the paper claim that some boxes are either suspended
% in the air or lack adequate support when testing the [R P]. To tackle this 
% issue, the bin is cut into horizontal slices, and bottom slices are enforced
% to have a higher fill/density than top slices. The object is modified to 
% penalize unfilled space with more weight given to bottom slices and to push
% the boxes towards the vertical axis. New variables are introduced. The new 
% relative-position formulation with slicing [R P S] is: 
%
%
% <html>
% <ul>
% <li>
% <img
% src="http://latex.codecogs.com/gif.latex?B_{i}^{k}%20=%20\begin{cases}%20&%201,%20\quad%20\text{if%20part%20of%20box%20i%20is%20contained%20in%20slice%20k}\\%20&%200,%20\quad%20\text{otherwise.}%20\end{cases}" height="36">
% </li>
% <li>
% <img
% src="http://latex.codecogs.com/gif.latex?S_{k}:%20\quad%20\text{the%20total%20area%20occupied%20by%20boxes%20in%20slice%20k}"height="12">
% </li>
% </ul>
% </html>
%
% $$[R\quad P\quad S]: min \quad h-\sum_{k=1}^{K}{(\frac{K-k}{L^{x}L^{y}})S_{k} + \sum_{i
% \in N}{(c_{i}^{x} + c_{i}{y})} }$$
%
% $$s.t. (1), (2), (3), (4), (5), (6), (7) \quad and$$
%
% $$-c_{i}^{z} + (s_{k} - l_{i}^{z})B_{i}^{k} \leq 0 \qquad \forall i,k \qquad (8)$$
%
% $$c_{i}^{z} + L^{z}B_{i}^{k} \leq s_{k} - \epsilon + L^{z} \qquad \forall i,k \qquad (9)$$
%
% $$\sum_{i \in N}{B_{i}^{k}(l_{i}^{x}l_{i}^{y}) = S_{k}} \qquad \forall k \qquad (10)$$
%
% $$-S_{k} + S_{k + 1} \leq 0 \qquad \forall k = 1......K - 1 \qquad (11)$$
%
% $$\sum_{k = 1}^{K}{B_{i}^{k} = l_{i}^{z}} \qquad \forall i \qquad (12)$$
% 
% $$B_{i}^{k} \in \{0, 1\} \qquad \forall i,k \qquad (13)$$
% 
%% Data sources
% An illustration of the results of model [R P] and [ R P S] on an example 
% with 30 boxes and a pallet of size $6\times 6\times15$ is conducted in the paper. The 
% box vary in size with length, width and height taking value in {1, 2, 3, 4}
% , {1, 2, 3}, and {1, 2, 3}, respectively. A computational time limit of 200s
% is used. the _rand_ function in MATLAB is used to randomly generate the length, 
% width,  and height for the boxes. 

%% Solution using Cplex Solver
%
% In the paper, the numerical test is implemented in MATLAB but using Cplex 
% 12.6 as the solver. In this project, the modelling part is done in Oplide,
% the Cplex IDE using Optimization Programming Language (OPL). The OPL code 
% can be found in the appendix at the end of this document. MATLAB is used to
% read the model developed in Oplide, to call the Cplex solver, and to visualize
% the final packing result. The MATLAB code for the visualization has also been 
% attached in the appendix.
%
% Both the [R P] and [R P S] formulation are implemented. The Matlab codes,
% results and the corresponding visualizations are shown below.
%%
% 
% * packing by the relative-position formulation [R P]
cplex=Cplex();
cplex.readModel('cplexbpp.lp');
cplex.Param.timelimit.Cur=200;
cplex.solve();
N=30;
csolution=vec2mat(cplex.Solution.x(2702:2791),3);
% csolution=vec2mat(cplex.Solution.x(2714:2803),3);
l =[1,	2	3;2 1	4;3	4	1;4	3	2;3	4	1;4	3	2;
    2	3	4;3	2	1;4	3	1;1	2	3;2	1	4;3	4	1;
    4	3	2;3	4	1;4	3	2;2	3	4;3	2	1;4	3	1;
    1	1	1;2	2	2;3	3	3;1	2	2;2	1	2;1	1	3;
    1	3	1;1	1	1;2	2	2;3	3	3;1	2	2;2	1	2];
figure(1)
for i=1:N
    clr=[rand rand rand];
    plotcube(l(i,:),csolution(i,:),1,clr);
    hold on
end
disp(figure(1));
%%
%
% * packing by the relative-position formulation [R P S]
%%
% 
% 
cplex2=Cplex();
cplex2.readModel('binpacking2.lp');
cplex2.Param.timelimit.Cur=200;
cplex2.solve();
N=30;
% csolution=vec2mat(cplex.Solution.x(2702:2791),3);
csolution2=vec2mat(cplex2.Solution.x(2714:2803),3);
l =[1,	2	3;2 1	4;3	4	1;4	3	2;3	4	1;4	3	2;
    2	3	4;3	2	1;4	3	1;1	2	3;2	1	4;3	4	1;
    4	3	2;3	4	1;4	3	2;2	3	4;3	2	1;4	3	1;
    1	1	1;2	2	2;3	3	3;1	2	2;2	1	2;1	1	3;
    1	3	1;1	1	1;2	2	2;3	3	3;1	2	2;2	1	2];

%%
figure(2)
for i=1:N
    clr=[rand rand rand];
    plotcube(l(i,:),csolution2(i,:),1,clr);
    hold on
end
disp(figure(2));
%% Solution using MATLAB solver

% By curiosity, an equivalent MATLAB code using MATLAB optimization solver is
% developed. The code is shown below. However, due to the limitation of the 
% MATLAB solver,  the MIP can not converge when trying to pack 30 boxes using 
% the [R P] formulation. Only the problem with a smaller number of boxes, 
% for example, ten boxes,  can be solved using MATLAB solver, and the processing 
% time is much slower compared to the Cplex solver.


n=10; % the number of the boxes
lx=ceil(4*rand(1,n));
ly=ceil(3*rand(1,n));
lz=ceil(3*rand(1,n));
l=[lx;ly;lz]; % the dimension of the boxes,lx=length,ly=width,lz=hight
LX=6;
LY=6;
LZ=15;

bpp=optimproblem;

h=optimvar('h','type','continuous','lowerbound',0);
cname=["cx","cy","cz"];
c=optimvar('c',cname,n,'type','continuous','lowerbound',0);
pname=["px","py","pz"];
p=optimvar("p",pname,n,n,'type','integer','lowerbound',0,'upperbound',1);

bpp.Objective=h;

bpp.Constraints.maxhight=c('cz',:)+lz(1,:)<=h;
bpp.Constraints.palletboundary=c<=[LX;LY;LZ]-l;
spatialrelation=optimexpr(1,n/2*(n-1));
orderindirection=optimexpr(3,n/2*(n-1));
nonoverlapping=optimexpr(3,n^2-n);

m=1;
for i=1:n-1
    for j=i+1:n
        spatialrelation(m)=p('px',i,j)+p('px',j,i)+p('py',i,j)+p('py',j,i)+p('pz',i,j)+p('pz',j,i);
        orderindirection(1,m)=p('px',i,j)+p('px',j,i);
        orderindirection(2,m)=p('py',i,j)+p('py',j,i);
        orderindirection(3,m)=p('pz',i,j)+p('pz',j,i);
        m=m+1;
    end 
end 
m=1;
for i=1:n
    for j=1:n
        if j~=i
          nonoverlapping(1,m)=c('cx',i)+lx(1,i)-c('cx',j)-LX*(1-p('px',i,j));
          nonoverlapping(2,m)=c('cy',i)+ly(1,i)-c('cy',j)-LY*(1-p('py',i,j));
          nonoverlapping(3,m)=c('cz',i)+lz(1,i)-c('cz',j)-LZ*(1-p('pz',i,j));
          m=m+1;
        end
    end 
end 


bpp.Constraints.spatialrelation=spatialrelation>=1;
bpp.Constraints.orderindirection=orderindirection<=1;
bpp.Constraints.nonoverlapping=nonoverlapping<=0;

writeproblem(bpp,'bpp10.lp');
options = cplexoptimset('cplex');
options.Display = 'on';
options.mip.tolerances.mipgap =  0;
options.parameter2058 = 1;
options.mip.pool.relgap= 0;
options.parameter2104 = 2;
opt.exportmodel = 'myModel.lp';
[sol,fval] = solve(bpp);
figure(3)
for i=1:n
    clr=[rand rand rand];
    plotcube(transpose(l(:,i)),transpose(sol.c(:,i)),1,clr);
    hold on
end
disp(figure(3));
%% Analysis and conclusion
% The implementations of the MIP formulation, both the [R P] and [R P S], 
% appear to have successfully packed the boxes into the bin without overlapping. 
% 
% Consistent with the description in the article, the [R P] formulation converge 
% to an optimum solution. However, there exists some boxes that are suspended 
% in the air.
%
% The packing using the [R P S] formulation has the same minimized height and
% is relatively more compact in bottom slices. However, the result is not as 
% ideal as the result of the original paper. There still exists unfilled space
% in the bottom slices, and boxes without support. The reason is that the [R P S] 
% formulation aims to enforce the fill rate of the bottom slices. There is no
% guarantee that the unfilled space will be completely eliminated and that 
% each box will get full support. Another reason is that I am implementing 
% the MIP on my Macbook air. The result may not be as ideal as the result 
% obtained from the computer with an Intel Core i7 processor under the same 
% time limit. 
% 
% From this project, we can also intuitively feel the difference between 
% MATLAB and CPLEX solver in computing speed for MIP by comparing the previous
% two sections.
%% References
%
% Elhedhli, S., Gzara, F., & Yan, Y. F. (2017). A MIP-based slicing heuristic 
% for three-dimensional bin packing. Optimization Letters, 11(8), 1547-1563.
%% Appendix
%
% * Visualization function
%
%%
% 
% 
% function plotcube(varargin)
% PLOTCUBE - Display a 3D-cube in the current axes
%
%    PLOTCUBE(EDGES,ORIGIN,ALPHA,COLOR) displays a 3D-cube in the current axes
%    with the following properties:
%   * EDGES : 3-elements vector that defines the length of cube edges
%    * ORIGIN: 3-elements vector that defines the start point of the cube
%    * ALPHA : scalar that defines the transparency of the cube faces (from 0
%              to 1)
%    * COLOR : 3-elements vector that defines the faces color of the cube
% 
%  Example:
%    >> plotcube([5 5 5],[ 2  2  2],.8,[1 0 0]);
%    >> plotcube([5 5 5],[10 10 10],.8,[0 1 0]);
%    >> plotcube([5 5 5],[20 20 20],.8,[0 0 1]);
%  Default input arguments
% inArgs = { ...
%   [10 56 100] , ... % Default edge sizes (x,y and z)
%   [10 10  10] , ... % Default coordinates of the origin point of the cube
%   .7          , ... % Default alpha value for the cube's faces
%   [1 0 0]       ... % Default Color for the cube
%   };
% % Replace default input arguments by input values
% inArgs(1:nargin) = varargin;
% % Create all variables
% [edges,origin,alpha,clr] = deal(inArgs{:});
% XYZ = { ...
%   [0 0 0 0]  [0 0 1 1]  [0 1 1 0] ; ...
%   [1 1 1 1]  [0 0 1 1]  [0 1 1 0] ; ...
%   [0 1 1 0]  [0 0 0 0]  [0 0 1 1] ; ...
%   [0 1 1 0]  [1 1 1 1]  [0 0 1 1] ; ...
%   [0 1 1 0]  [0 0 1 1]  [0 0 0 0] ; ...
%   [0 1 1 0]  [0 0 1 1]  [1 1 1 1]   ...
%   };
% XYZ = mat2cell(cellfun( @(x,y,z) x*y+z, XYZ, repmat(mat2cell(edges,1,[1 1 1]),6,1), repmat(mat2cell(origin,1,[1 1 1]),6,1),'UniformOutput',false),6,[1 1 1]);
% cellfun(@patch,XYZ{1},XYZ{2},XYZ{3},...
%   repmat({clr},6,1),...
%   repmat({'FaceAlpha'},6,1),...
%   repmat({alpha},6,1)...
%   );
% view(3);
% end
%%
%
% * Modeling using Cplex Optimization Programming Language (OPL) 
%
%  *********************************************
%
%  * OPL 12.9.0.0 Model
%
%  * Author: zhaoyinghui
%
%  * Creation Date: Apr 24, 2019 at 6:49:15 PM
%
%  *********************************************/
% //Parameters 
%  int n=...;
%  range Box=1..n;
%  range dimension=1..3;
%  int K=...;
%  range slice=1..K;
%  int s[slice]=...;
%  float e=1e-2;
%  
%  int l[Box][dimension]=...;
%  int L[dimension]=...;
%  
%  // variables
%  dvar boolean p[dimension][Box][Box];
%  dvar float h;
%  dvar boolean b[slice][Box];
%  dvar int S[slice];
%  dvar float c[Box][dimension];
%  
%  execute
%  {
%     cplex.tilim = 200;
% 
%  }
%  // objectives
%  minimize h-(sum(k in 1..K )S[k]*(K-k)/L[1]/L[2])+(sum(i in Box)(c[i][1]+c[i][2]));
%  
% 
%  //constraints
%  subject to {
%  	forall(j in Box , i in Box) 
%  	 if (j>i)
%  	 constraint_1:
%  	 sum (k in dimension) p[k][j][i] >= 1;
%  	
%  	forall(k in dimension, j in Box, i in Box)
%  	  if (j>i)
%  	  constraint_2:
%  	  p[k][j][i]+p[k][i][j]<=1;
%  	
%  	forall(k in dimension, j in Box, i in Box)
%  	  if (j != i)
%  	  constraint_3:
%  	  c[i][k]+l[i][k]<=c[j][k]+L[k]*(1-p[k][i][j]);
%  	  
%  	forall(i in Box, k in dimension)
%  	   constraint_4:
%  	   0<=c[i][k]<=L[k]-l[i][k];
%  	   
%  	forall(i in Box)
%  	  constraint_5:
%  	   c[i][3]+l[i][3]<=h;
%  	  
%  	  constraint_7:  
%  		h >=0;
%  		
%  	forall(i in Box, k in slice)
%  	  constraint_8:
%  	  -c[i][3]+(s[k]-l[i][3])*b[k][i]<=0;
%  	  
%    forall (i in Box, k in slice)
%      constraint_9:
%      c[i][3]+L[3]*b[k][i]<=s[k]-e+L[3];
%      
%    forall (k in slice)
%      constraints_10:
%      sum(i in Box) b[k][i]*(l[i][1]*l[i][2])== S[k];
%      
%    forall (k in 1..K-1)
%      constraint_11:
%      -S[k]+S[k+1]<=0;
%      
%    forall(i in Box)
%      constraint_12:
%      sum(k in 1..K) b[k][i]==l[i][3];
%   } 	
%  	
%  
%  // postprocessing
% 
%   execute{
%   if (cplex.getCplexStatus()==1){
% 	var f=new IloOplOutputFile("solution.csv");
% 	f.writeln(h);
% 	for (i in Box)
% 	{f.writeln(c[i]);}  	
%     f.close();}
%   else {writeln("error, solution not found");}
%   }