由于要数学建模中需要解决一个三维装箱的问题，我经过搜寻选定了张德富教授等在计算机学报上发表的《求解三维装箱问题的混合模拟退火算法》这篇论文作为解决问题的理论基础。
该论文的摘要如下：
提出了一个高效求解三维装箱问题(Three Dimensional Container Loading Problem 3D-CLP)的混合模拟 退火 算法 .三维装箱问题要求装载给定箱子集合的一个子 集到容器中 , 使得被装载 的箱子总体 积最大 .文中介绍 的 混合 模拟退火算法基于三个重要算法 :(1)复合块生成算法 , 与传统算法 不同的是文中提出的复合块 不只包含单 一 种类的箱子 , 而是可以在一定的限制条件下包含 任意种类的箱子 .(2)基础启发式算法 , 该算法 基于块装 载 , 可以 按 照指定装载序列生成放置方案.(3)模拟退火算法,以复合块生成和基础启发式算法为基础, 将装载序列作为可行 放置 方案的编码 , 在编码空间中采用模拟退火算法进行搜 索以寻找问题的近似最 优解 .文 中采用 1 50 0 个弱异构 和 强异构的装箱问题数据对算法进行测试 .实验结 果表明 , 混合模拟退火算法的填充率超过了目前已知的 优秀算法 .
DOI 号 :10.3724/ SP .J.1016 .2009.02147
本文算法实现其计算结果
该实现的核心函数MATLAB代码
function \[solution, newPlan, rate\] = container\_packeting(bcon, plan, blockTable, blockNeed)
%UNTITLED2 根据大小容器的规格以及小容器的数量信息计算出打包的方案
% bigContainer 大容器的规格向量 \[L W H\]
% smallContainers 小容器的规格矩阵\[L W H;\]
% smallContainersNumber 小容器的数量向量 \[n\]
solution = {};
newPlan = \[\];
searchDeep = 40960;
clf;
draw\_container(\[0,0,0\], bcon);
% 剩余空间栈
lspace\_s = get\_stack();
% 空隙栈
iterval\_s = get\_stack();
itervalSpace\_s = get\_stack();
% 初始化剩余空间为大容器空间
lspace\_s = stack\_push(lspace\_s, get\_space(\[0,0,0\], bcon));
psCount = 0;
while lspace\_s.count > 0
psCount = psCount + 1;
% 从剩余空间栈中取得可用的剩余空间
lspace = stack\_top(lspace\_s);
lspace\_s = stack\_pop(lspace\_s);
lspace\_size = lspace.size;
lspace\_base = lspace.base\_point;
if\_input = 0;
\[blockList, blockListNeed\] = get\_block\_list(blockTable, blockNeed, lspace\_size(1), lspace\_size(2), lspace\_size(3));
theBlockNeed = NaN;
find\_idx = 1;
if size(blockList, 1) == 0
find\_avaliable = 0;
else
find\_avaliable = 1;
if plan(psCount) <= size(blockListNeed, 1)
theBlockNeed = blockListNeed(plan(psCount),:);
find\_idx = plan(psCount);
else
theBlockNeed = blockListNeed(size(blockListNeed, 1),:);
find\_idx = size(blockListNeed, 1);
end
newPlan = \[newPlan, find\_idx\];
end
%% locate container
if find\_avaliable == 1
newLocatedBlock = {lspace\_base, blockList(find\_idx, :), theBlockNeed{1}, theBlockNeed{2}};
solution(size(solution, 1) + 1, :) = newLocatedBlock;
spaceInfo = newLocatedBlock{4};
spaceLastInfo = spaceInfo(size(spaceInfo,1), :);
locateInfo = newLocatedBlock{2};
hold on;
draw\_container(lspace\_base, locateInfo(1:3));
itervalSpace\_s = stack\_push(itervalSpace\_s, (1-spaceLastInfo(1)) \* prod(locateInfo(1:3)));
if\_input = 1;
% 剩余空间计算
% H方向剩余空间
location = lspace\_base;
location(2) = location(2) + locateInfo(2);
iterval\_space = get\_space(location, \[locateInfo(1) , lspace\_size(2) - locateInfo(2), locateInfo(3)\]);
lspace\_s = stack\_push(lspace\_s, iterval\_space);
% L方向剩余空间
location = lspace\_base;
location(1) = location(1) + locateInfo(1);
left\_space = get\_space(location, \[lspace\_size(1) - locateInfo(1), lspace\_size(2), lspace\_size(3)\]);
lspace\_s = stack\_push(lspace\_s, left\_space);
% W方向剩余空间
location = lspace\_base;
location(3) = location(3) + locateInfo(3);
iterval\_space = get\_space(location, \[locateInfo(1) , lspace\_size(2), lspace\_size(3) - locateInfo(3)\]);
lspace\_s = stack\_push(lspace\_s, iterval\_space);
end
if if\_input == 1
else
% 压入空隙栈
iterval\_s = stack\_push(iterval\_s, lspace);
% lspace\_s = stack\_push(lspace\_s, lspace);
end
end
rate = space\_used\_rate(bcon, solution);
end
function rate = space\_used\_rate(container, solution)
spaceHasTotal = 0.0;
for i = 1:1:size(solution, 1)
mergeMessage = solution{i, 4};
for k = 1:1:size(mergeMessage,1)
spaceHasTotal = spaceHasTotal + prod(mergeMessage(k, 3:5));
end
end
rate = spaceHasTotal / prod(container);
end
箱体产生函数
function \[blockTable, boxNeedTable\] = simple\_block\_generate\_indepent(container, box, num)
simBlockCount = 0;
simBlockTable = zeros(256, 8);
for type = 1:1:size(box, 1)
for nx = 1:1:num(type)
for ny = 1:1:num(type)/nx
for nz = 1:1:num(type)/ny/nx
if box(type,3) \* nx < container(3) && box(type, 1) \* ny < container(1) && box(type, 2) \* nz < container(2)
newSimBlock = \[box(type, :), nx, ny, nz, type, prod(box(type, :)) \* nx \* ny \* nz\];
simBlockTable(simBlockCount + 1, :) = newSimBlock;
simBlockCount = simBlockCount + 1;
end
end
end
end
end
simBlockTable = simBlockTable(1:simBlockCount, :);
simBlockTable = sortrows(simBlockTable, 8);
simBlockTable = simBlockTable(:, 1:7);
% 处理成启发式算法可以处理的格式
simBlockTable(:, 13) = ones(size(simBlockTable, 1), 1);
simBlockTable(:, 11) = simBlockTable(:,4) .\* simBlockTable(:,3);
simBlockTable(:, 12) = simBlockTable(:,5) .\* simBlockTable(:,1);
% ly = ay
simBlockTable(:, 8) = simBlockTable(:, 12);
% lx = ax
simBlockTable(:, 10) = simBlockTable(:, 11);
% lz = lz0 \* nz
simBlockTable(:, 9) = simBlockTable(:, 2) .\* simBlockTable(:, 6);
boxNeedTable = cell(size(simBlockTable, 1), 2);
for i = 1:1:size(simBlockTable, 1)
boxNeedTable{i, 1} = zeros(1, size(box ,1));
boxNeedTable{i, 1}(simBlockTable(i, 7)) = simBlockTable(i, 4) \* simBlockTable(i, 5) \* simBlockTable(i, 6);
boxNeedTable{i, 2} = \[simBlockTable(i, 7), 0, simBlockTable(i, 8:10)\];
end
blockTable = simBlockTable(:, 8:13);
blockTable(:,7) = blockTable(:, 1) .\* blockTable(:, 2) .\* blockTable(:, 3);
\[blockTable, idx\] = sortrows(blockTable, \[ 7 \]);
boxNeedTable = boxNeedTable(idx, :);
end
箱体组合产生函数
function \[blockTable, boxNeedTable\] = complex\_block\_genrate(container, box, num, level)
simBlockTable = simple\_block\_generate(container, box, num);
simBlockTable(:, 13) = ones(size(simBlockTable, 1), 1);
simBlockTable(:, 11) = simBlockTable(:,4) .\* simBlockTable(:,3);
simBlockTable(:, 12) = simBlockTable(:,5) .\* simBlockTable(:,1);
spaceUsedRateMin = 0.90;
% ly = ay
simBlockTable(:, 8) = simBlockTable(:, 12);
% lx = ax
simBlockTable(:, 10) = simBlockTable(:, 11);
% lz = lz0 \* nz
simBlockTable(:, 9) = simBlockTable(:, 2) .\* simBlockTable(:, 6);
boxNeedTable = cell(size(simBlockTable, 1), 2);
for i = 1:1:size(simBlockTable, 1)
boxNeedTable{i, 1} = zeros(1, size(box ,1));
boxNeedTable{i, 1}(simBlockTable(i, 7)) = simBlockTable(i, 4) \* simBlockTable(i, 5) \* simBlockTable(i, 6);
boxNeedTable{i, 2} = \[simBlockTable(i, 7), 0, simBlockTable(i, 8:10)\];
end
blockTable = simBlockTable(:, 8:13);
for level = 1:1:level
newBlockTable = zeros(128,6);
newBlockCount = 1;
boxNeedCount = size(boxNeedTable, 1);
for a = 1:1:size(blockTable,1)
for b = a:1:size(blockTable, 1)
if b == a
continue;
end
if blockTable(a, 6) == blockTable(b, 6)
% x 方向合并
if blockTable(a, 1) == blockTable(a, 5) && blockTable(a, 3) == blockTable(a, 4) && blockTable(a, 2) == blockTable(b, 2)
tempBlockMerged = \[max(blockTable(a,1), blockTable(b,1)), max(blockTable(a,2), blockTable(b, 2)), blockTable(a,3) + blockTable(b,3), blockTable(a,4) + blockTable(b,4), min(blockTable(a,5), blockTable(b,5)), max(blockTable(a, 6), blockTable(b, 6)) + 1\];
rate = space\_used\_rate(tempBlockMerged(1:3), blockTable(a,1:3), blockTable(b,1:3));
if rate > spaceUsedRateMin
newBlockTable(newBlockCount, :) = tempBlockMerged;
boxNeedCount = boxNeedCount + 1;
newBlockCount = newBlockCount + 1;
boxNeedTable{boxNeedCount, 1} = boxNeedTable{a,1} + boxNeedTable{b,1};
boxNeedTable{boxNeedCount, 2} = \[boxNeedTable{a, 2}; rate, 1, blockTable(b, 1:3)\];
end
end
% y 方向合并
if blockTable(a, 5) == blockTable(a, 1) && blockTable(b, 5) == blockTable(b, 1) && blockTable(a, 2) == blockTable(b, 2)
tempBlockMerged = \[blockTable(a,1) + blockTable(b,1), max(blockTable(a, 2), blockTable(b, 2)), max(blockTable(a, 3), blockTable(b, 3)), min(blockTable(a, 4), blockTable(b, 4)), blockTable(a, 5) + blockTable(b, 5), max(blockTable(a, 6), blockTable(b, 6)) + 1\];
rate = space\_used\_rate(tempBlockMerged(1:3), blockTable(a,1:3), blockTable(b,1:3));
if rate > spaceUsedRateMin
newBlockTable(newBlockCount, :) = tempBlockMerged;
newBlockCount = newBlockCount + 1;
boxNeedCount = boxNeedCount + 1;
boxNeedTable{boxNeedCount, 1} = boxNeedTable{a,1} + boxNeedTable{b,1};
boxNeedTable{boxNeedCount, 2} = \[boxNeedTable{a, 2}; rate, 2, blockTable(b, 1:3)\];
end
end
% z 方向合并
if blockTable(a,4) >= blockTable(b,3) && blockTable(a, 5) >= blockTable(b, 1)
tempBlockMerged = \[max(blockTable(a,1), blockTable(b, 1)), blockTable(a, 2) + blockTable(b, 2), max(blockTable(a, 3), blockTable(b, 3)), blockTable(b, 4), blockTable(b, 5), max(blockTable(a, 6), blockTable(b, 6)) + 1\];
rate = space\_used\_rate(tempBlockMerged(1:3), blockTable(a,1:3), blockTable(b,1:3));
if rate > spaceUsedRateMin
newBlockTable(newBlockCount, :) = tempBlockMerged;
newBlockCount = newBlockCount + 1;
boxNeedCount = boxNeedCount + 1;
boxNeedTable{boxNeedCount, 1} = boxNeedTable{a,1} + boxNeedTable{b,1};
boxNeedTable{boxNeedCount, 2} = \[boxNeedTable{a, 2}; rate, 3, blockTable(b, 1:3)\];
end
end
end
end
end
blockTable = cat(1, blockTable, newBlockTable(1:newBlockCount-1, :));
% 消除等价复杂块
blockTableTemp = blockTable(:, 1:5);
\[~, ia\] = unique(blockTableTemp, 'stable', 'rows');
blockTable = blockTable(ia, :);
boxNeedTable = boxNeedTable(ia, :);
end
blockTable(:,7) = blockTable(:, 1) .\* blockTable(:, 2) .\* blockTable(:, 3);
\[blockTable, idx\] = sortrows(blockTable, \[ 7 \]);
boxNeedTable = boxNeedTable(idx, :);
end
function \[rate\] = space\_used\_rate(container, box1, box2)
rate = (prod(box1) + prod(box2)) /prod(container);
end
Main函数
ts = 1;
tf = 0.005;
dt = 0.98;
length = 50;
maxSeq = 64;
maxChoice = 1024;
click = 1;
tempRate = \[\];
tempRateX = \[\];
numList = \[64,64,64,64,64,64,64,128,128,128\];
container = StandardISOContainer;
% \[blockTable, blockNeed\] = complex\_block\_genrate(StandardISOContainer, Packets, numList, 1);
ps = ones(1, maxSeq);
ps(1) = 1;
\[solution, plan, rate\] = container\_packeting(container, ps, blockTable, blockNeed);
best = plan;
best\_rate = rate;
best\_solution = solution;
t = ts;
while t > tf
for i = 1:1:length
k = randi(size(ps));
nps = ps;
nps(k) = randi(maxChoice);
\[solution, nplan, nrate\] = container\_packeting(container, nps, blockTable, blockNeed);
% disp(\['TempPlan: ', mat2str(nplan)\]);
% fprintf("TempRate: %f\\n",nrate);
click = click + 1;
if nrate > rate
ps = nps;
plan = nplan;
rate = nrate;
tempRateX = \[tempRateX, click\];
tempRate = \[tempRate, nrate\];
elseif rand(0, 1) > exp((nrate - rate) \* 10 / t)
ps = nps;
plan = nplan;
rate = nrate;
end
if nrate > best\_rate
best = nplan;
best\_rate = nrate;
best\_solution = solution;
disp(\['BestPlan: ', mat2str(best)\]);
fprintf("BestRate: %f\\n",best\_rate);
end
end
t = (1 - t \* dt) \* t;
end
disp(\['BestPlan: ', mat2str(best)\]);
fprintf("BestRate: %f\\n",best\_rate);
figure(1)
draw\_solution(container, best\_solution);
figure(2)
plot(tempRateX, tempRate);
save Solution best\_solution best\_rate tempRateX tempRate