一、三维装箱简介
1 前言
三维装箱问题（Three Dimensional Container Loading Poblem）指在满足容积限制、外形几何限制和稳定性限制等条件的情况下，把一定数量体积较小的物品放入具有较大容量的一个或多个箱子，达到所用箱子数量最少、空间利用率最高、稳定性最好、装载价值最高、容重比最高等目的组合优化的问题。三维装箱问题广泛存在于运输和包装等行业的物资装载过程中，因此，如何针对实际应用需求设计具有较高操作性和装载效率的三维装箱求解算法，对于提高社会经济效益具有积极作用。
自20世纪60年代开始，三维装箱问题便因切割库存问题而提出，引起了学者关注，并被证明是具有较高复杂度的NP难问题。经过多年发展，三维装箱问题的求解不仅在枚举法、分支定界法等确定性算法上取得了进步。随着元启发式算法、遗传算法、粒子群算法等现代启发式算法的崛起，三维装箱问题的研究更加深入，逐步出现了带有多种现实约束和求解目标的复杂算法。其中，近年来具有代表性的相关研究如下所述。
Lodi等[9]受到一维装箱问题求解算法的启发，针对三维装箱问题提出了基于“先高后面积”（Height First-Area Second）的分层策略的结构式禁忌搜索算法，实验证明算法在分层堆放的情况下能得到较好的结果。Crainic等考虑了装载盒子时箱子剩余空间的描述和备选定位点的重要作用，为寻找每次装载时的合适空间位置和最优定位点，提出了极点（Extreme Points,EPS）的概念和计算方法，并设计了基于极点的三维装箱启发式算法。Karabulut[12]和Kang[13]等先后以规则长方体箱装载最大量规则长方体盒为研究目标，在设计深度、底部和左部方向优先装载（Deepest Bottom Left with Fill）物资装载策略的基础上，提出了基于遗传算法的三维装箱策略求解算法，2个研究的区别在于后者在前者的基础上，对装载剩余空间进行了更为细致的描述。Ramos等重点考虑了运输行业货物装载时静止时的纵向稳定性和行进时的横向稳定性，并提出了具有较高操作性的稳定性评价指标和基于顺序的货物装载启发式算法。
尽管三维装箱问题研究的深度和广度在不断延伸，相关研究的模型建立、算法设计和实际应用也取得了众多成果，但是一些现实约束仍难以满足，解决方法也存在较大的局限性。其中，现有求解三维装箱问题较为常用的遗传算法，在交叉算子设计、变异算子设计和进化性能等方面表现不足。
1 问题描述
1.1 基本问题描述
除球体、非规则立方体等形状盒子的装箱特例外，经典的三维装箱问题（见图1）可描述为：用一定数量的长宽高分别为Li,Wi和Hi的箱子（Container），按照一定顺序和装载规则根据不同转向逐个装入一定数量的长、宽、高分别为lj,wj和hj的盒子（Box）中，实现箱子容量的利用率最大、盒子堆放的稳定性最好、装载价值最高或容重比最恰当等其一或多种目标，并满足外形约束、容积约束、稳定性约束、盒子属性相关性约束、载重约束等实际限制。
图1 三维装箱问题示意
根据文献对三维装箱问题的分类法，该研究以SBSBPP类（Single Bin-Size Bin Packing Problem，有限数量单类型箱体三维装载问题）为研究对象，研究如何在单个规则长方体箱子中装载尽量多的盒子（这些盒子具有相似度较小且分散较大的特点），以最大限度地提高箱子的空间利用率。此研究旨在解决每个盒子应以何种顺序、何种方向放置在具体空间位置的问题，即提供具体的盒子装载方案。此过程需满足如下现实约束。
1）外形约束。进行装载时，所有盒子的几何顶点不能位于箱子之外，装进箱子的盒子在空间上不能发生重叠。
2）容积约束。装载盒子的总体积不能超过箱子的容积。
3）稳定性约束。装载盒子的堆放需满足重心要求，不能出现悬空或重心偏移的情况。
基本假设：盒子质量均匀分布，其质心位于其几何中心；不考虑盒子总质量超载的问题；盒子装载时应与箱子的边长方向平行（即不斜放）。
1.2 箱体最大剩余空间描述
最大剩余空间（Allowable Packing Area，简写为APA）指进行盒子装载前的最大长方体形闲置空间，每个最大剩余空间可用最靠近原点的长方体顶点三维坐标和最远离原点的长方体顶点三维坐标合并表示，全部的最大剩余空间即构成最大剩余空间集（简写为APAs）。盒子的装载是在箱子的最大剩余空间中按照一定规则进行的，最大剩余空间的刻画关系到盒子的放置和空间容纳能力，对于盒子的装载至关重要。该研究采用文献中的I-DBLF算法，完成最大剩余空间的描述。每装入一个箱子，其最大剩余空间将按照盒子的边界和x,y,z轴进行划分和更新。如图2a所示，空置箱子的最大剩余空间只有1个，表示为(0,0,0),(Li,Wi,Hi)；在临原点位置装入一个长、宽、高分别为lj,wj和hj的盒子后，最大剩余空间更新为3个（见图2b—d），最大剩余空间集为{((Li-lj,0,0),(Li,Wi,Hi)),((0,0,Hi-hj),(Li,Wi,Hi)),((0,Wi-wj,0),(Li,Wi,Hi))}。最大剩余空间的计算过程较为复杂，因篇幅所限，此处不做赘述，具体可见文献。
图2 剩余空间示意
二、部分源代码
close all;
clear;
clc;
%% set data
Dbox = [10, 10, 10];
Dobj = [1, 1.5, 1;
1, 1.5, 1];
PobjI = [7, 5, 5, 0, 0, 0;
3, 5, 5, 0, 0, 0];
PobjD = [3, 3, -3, 0, 0, 0;
8, 8, 13, 0, 0, 0];
objN = size(Dobj, 1);
Rcost = 10;
%% shape & bounding box & start point
[x, f] = genFunction(Rcost, PobjD);
Ar = [1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1];
syms vx vy vz;
Vcoor = [vx; vy; vz];
for i=1:objN
br = [Dobj(i,1); Dobj(i,1); Dobj(i,2); Dobj(i,2); Dobj(i,3); Dobj(i,3)];
Gc{i} = Ar * Vcoor - br;
end
G = zeros(0);
for i=1:objN
G = [G; genGi(i, Dbox, Dobj, x)];
end
%%
tic;
[X0,~] = zoutendijk1(f, x, G, Gc, Dbox, Dobj, PobjI, PobjD);
toc
figure;
plotResult(X0, Dbox, Dobj, PobjI, PobjD);
function [ f ] = genCollision( Gc, X, x0, Dobj, o1, o2 )
syms vx vy vz;
x = [vx vy vz];
Gc1 = Gc{o1};
Gc2 = Gc{o2};
a1 = Dobj(o1,1);
b1 = Dobj(o1,2);
c1 = Dobj(o1,3);
a2 = Dobj(o2,1);
b2 = Dobj(o2,2);
c2 = Dobj(o2,3);
R1V = [a1 b1 c1; -a1 -b1 c1; -a1 b1 c1; a1 -b1 c1;
a1 b1 -c1; -a1 -b1 -c1; -a1 b1 -c1; a1 -b1 -c1];
R2V = [a2 b2 c2; -a2 -b2 c2; -a2 b2 c2; a2 -b2 c2;
a2 b2 -c2; -a2 -b2 -c2; -a2 b2 -c2; a2 -b2 -c2];
% R1V(1,:) = [-a1, -b1];
% R1V(2,:) = [a1, b1];
% R1V(3,:) = [-a1, b1];
% R1V(4,:) = [a1, -b1];
%
% R2V(1,:) = [-a2, -b2];
% R2V(2,:) = [a2, b2];
% R2V(3,:) = [-a2, b2];
% R2V(4,:) = [a2, -b2];
newR1V = R1V;
newR2V = R2V;
for i=1:size(R1V,1)
newR1V(i,:) = rt(R1V(i,:), transpose(x0(6*(o1-1)+1:6*(o1-1)+3)), transpose(x0(6*(o2-1)+1:6*(o2-1)+3)), x0(6*(o1-1)+4:6*(o1-1)+6), x0(6*(o2-1)+4:6*(o2-1)+6));
newR2V(i,:) = rt(R2V(i,:), transpose(x0(6*(o2-1)+1:6*(o2-1)+3)), transpose(x0(6*(o1-1)+1:6*(o1-1)+3)), x0(6*(o2-1)+4:6*(o2-1)+6), x0(6*(o1-1)+4:6*(o1-1)+6));
end
temp = zeros(size(newR2V,1),3);
minEdged = zeros(0,3);
for i=1:size(Gc1,1)
for j=1:size(newR2V,1)
temp(j,:) = [subs(Gc1(i), x, newR2V(j,:)), i, j];
end
[~, vec] = sort(temp);
temp1 = temp(vec(:,1),:);
for k=1:size(newR2V,1)
if abs(temp1(k,1)-temp1(k+1,1)) < 0.1
continue
else
break
end
end
minEdged = [minEdged; temp1(1:k,:)];
end
[~, vec] = sort(minEdged, 'descend');
temp1 = minEdged(vec(:,1),:);
for k=1:size(temp1,1)
if abs(temp1(k,1)-temp1(k+1,1)) < 0.1
continue
else
break
end
end
maxEdged1 = temp1(1:k,:);
minEdged = zeros(0,3); % value % edge num % v num
for i=1:size(Gc1,1)
for j=1:size(newR2V,1)
temp(j,:) = [subs(Gc2(i), x, newR1V(j,:)), i, j];
end
[~, vec] = sort(temp);
temp1 = temp(vec(:,1),:);
for k=1:size(temp1,1)
if abs(temp1(k,1)-temp1(k+1,1)) < 0.1
continue
else
break
end
end
minEdged = [minEdged; temp1(1:k,:)];
end
[~, vec] = sort(minEdged, 'descend');
temp1 = minEdged(vec(:,1),:);
for k=1:size(temp1,1)
if abs(temp1(k,1)-temp1(k+1,1)) < 0.1
continue
else
break
end
end
maxEdged2 = temp1(1:k,:);
f = zeros(0,0);
if abs(maxEdged1(1,1)-maxEdged2(1,1)) < 0.1
for i=1:size(maxEdged1,1)
f = [f; -subs(Gc1(maxEdged1(i,2)), x, rt(R2V(maxEdged1(i,3),:), X(6*(o2-1)+1:6*(o2-1)+3), X(6*(o1-1)+1:6*(o1-1)+3), X(6*(o2-1)+4:6*(o2-1)+6), X(6*(o1-1)+4:6*(o1-1)+6)))];
end
for i=1:size(maxEdged2,1)
f = [f; -subs(Gc2(maxEdged2(i,2)), x, rt(R1V(maxEdged2(i,3),:), X(6*(o1-1)+1:6*(o1-1)+3), X(6*(o2-1)+1:6*(o2-1)+3), X(6*(o1-1)+4:6*(o1-1)+6), X(6*(o2-1)+4:6*(o2-1)+6)))];
end
elseif (maxEdged1(1,1) < maxEdged2(1,1))
for i=1:size(maxEdged2,1)
f = [f; -subs(Gc2(maxEdged2(i,2)), x, rt(R1V(maxEdged2(i,3),:), X(6*(o1-1)+1:6*(o1-1)+3), X(6*(o2-1)+1:6*(o2-1)+3), X(6*(o1-1)+4:6*(o1-1)+6), X(6*(o2-1)+4:6*(o2-1)+6)))];
end
else
for i=1:size(maxEdged1,1)
f = [f; -subs(Gc1(maxEdged1(i,2)), x, rt(R2V(maxEdged1(i,3),:), X(6*(o2-1)+1:6*(o2-1)+3), X(6*(o1-1)+1:6*(o1-1)+3), X(6*(o2-1)+4:6*(o2-1)+6), X(6*(o1-1)+4:6*(o1-1)+6)))];
end
end
return
end
三、运行结果
四、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
[3]周品.MATLAB 神经网络设计与应用[M].清华大学出版社，2013.
[4]陈明.MATLAB神经网络原理与实例精解[M].清华大学出版社，2013.
[5]方清城.MATLAB R2016a神经网络设计与应用28个案例分析[M].清华大学出版社，2018.
