As I need to solve a three-dimensional packing problem in mathematical modeling, I have selected the paper “Hybrid Simulated Annealing Algorithm for Solving Three-Dimensional Packing Problems” published by Prof. Zhang Defu et al. in the Journal of Computing as the theoretical basis of the problem after searching.
The abstract of the paper is as follows:
A hybrid simulated annealing algorithm for efficiently solving the Three Dimensional Container Loading Problem 3D-CLP is proposed . The Three Dimensional Container Loading Problem (3D-CLP) requires loading a subset of a given set of boxes into a container such that the total volume of the loaded boxes is maximized. The hybrid simulated annealing algorithm presented in this paper is based on three important algorithms: (1) the composite block generation algorithm, which differs from the traditional algorithms in that the composite block proposed in this paper does not contain only a single type of crates, but can contain any type of crates under certain constraints. (2) Basic heuristic algorithm, which is based on block loading and can generate a placement scheme according to a specified loading sequence. (3) Simulated annealing algorithm, based on composite block generation and basic heuristic algorithm, encodes the loading sequence as a feasible placement scheme, and searches the encoding space using simulated annealing algorithm to find the approximate optimal solution of the problem. The algorithm is tested using 1,500 weakly and strongly heterogeneous data for the packing problem. Experimental results show that the fill rate of the hybrid simulated annealing algorithm exceeds that of the best known algorithms.
DOI number :10.3724/ SP .J.1016 .2009.02147
The algorithm in this paper realizes its computational results
MATLAB code for the core function of this realization
function [solution, newPlan, rate] = container_packeting(bcon, plan, blockTable, blockNeed)
%UNTITLED2 Calculates the packing plan based on the specifications of the big and small containers and the information about the number of small containers
% bigContainer bigContainer's specification vector [L W H]
% smallContainers Specification matrix of small containers [L W H;].
% smallContainersNumber Vector of the number of small containers [n].
solution = {};
newPlan = [] ;
searchDeep = 40960;
clf.
draw_container([0,0,0], bcon);
% Remaining space stack
lspace_s = get_stack();
% void stack
iterval_s = get_stack(); % void stack
itervalSpace_s = get_stack();
% Initialize the remaining space as large container space
lspace_s = stack_push(lspace_s, get_space([0,0,0], bcon)); % Initialize the remaining space as large container space.
psCount = 0;
while lspace_s.count > 0
psCount = psCount + 1;
% Get the available remaining space from the remaining space stack
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)
if size(blockList, 1) == 0
find_avaliable = 0; else
if size(blockList, 1) == 0 find_avaliable = 0; else
find_avaliable = 1; if plan(psCount) <= size(blockList, 1)
if plan(psCount) <= size(blockListNeed, 1)
theBlockNeed = blockListNeed(plan(psCount),:);
find_idx = plan(psCount);
else
theBlockNeed = blockListNeed(size(blockListNeed, 1),:); find_idx = size(block(psCount),:); else
find_idx = size(blockListNeed, 1); end
end
newPlan = [newPlan, find_idx]; end
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;
% Remaining space calculation
% H-direction space left
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-direction remaining space
location = lspace_base;
location(1) = location(1) + locateInfo(1); % L-direction residual space
left_space = get_space(location, [lspace_size(1) - locateInfo(1), lspace_size(2), lspace_size(3)]);
lspace_s = stack_push(lspace_s, left_space);
% remaining space in W direction
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
end if if_input == 1
% Push into the gap stack
iterval_s = stack_push(iterval_s, lspace);
% lspace_s = stack_push(lspace_s, lspace);
end
end
rate = space_used_rate(bcon, solution); end
end
function rate = space_used_rate(container, solution)
spaceHasTotal = 0.0; for i = 1:1:size(solution, 1); end function
for i = 1:1:size(solution, 1)
mergeMessage = solution{i, 4}; for k = 1:1:size(solution, 1)
for k = 1:1:size(mergeMessage,1)
spaceHasTotal = spaceHasTotal + prod(mergeMessage(k, 3:5));
end
end
rate = spaceHasTotal / prod(container);
end
Box Generation Functions
function [blockTable, boxNeedTable] = simple_block_generate_indepent(container, box, num)
simple_block_generate_indepent(container, box, num); simBlockCount = 0;
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);
% Process into a format that the heuristic algorithm can handle
simBlockTable(:, 13) = ones(size(simBlockTable, 1), 1); % Process into a format that can be handled by the heuristic algorithm.
simBlockTable(:, 11) = simBlockTable(:,4) .* . simBlockTable(:,3); .
simBlockTable(:, 12) = simBlockTable(:,5) .* . simBlockTable(:,1);
% ly = ay
simBlockTable(:, 8) = simBlockTable(:, 12); % ly = ay
% lx = ax
simBlockTable(:, 10) = simBlockTable(:, 11); % lz = lz0 * nz0; % ly = ay
% 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} = zeros(1, size(box ,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(:, 7) = blockTable(:, 1) .* blockTable(:, 2) .* blockTable(:, 3);
[blockTable, idx] = sortrows(blockTable, [ 7 ]);
boxNeedTable = boxNeedTable(idx, :);
end
Box combination generator function
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); simBlockTable(:,2) = simBlockTable(:,3) .
spaceUsedRateMin = 0.90;
% ly = ay
simBlockTable(:, 8) = simBlockTable(:, 12); % ly = ax
% lx = ax
simBlockTable(:, 10) = simBlockTable(:, 11); % lz = lz0 * lz0 * lz0 * lz0
% 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} = zeros(1, size(box ,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
for level = 1:1:level
newBlockTable = zeros(128,6);
newBlockCount = 1; boxNeedCount = size(box)
newBlockCount = 1; boxNeedCount = size(boxNeedTable, 1); for a = 1:1:size
for a = 1:1:size(blockTable,1)
for b = a:1:size(blockTable, 1)
if b == a
continue; end
if b == a continue; end
if blockTable(a, 6) == blockTable(b, 6)
% x Direction merge
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(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 > spaceUsedRate(tempBlockMerged(1:3))
if rate > spaceUsedRateMin
newBlockTable(newBlockCount, :)) = tempBlockMerged; boxNeedCount = boxNeedCount, :)); if rate > spaceUsedRateMin
boxNeedCount = boxNeedCount + 1; newBlockCount = newBlockCount; newBlockCount = newBlockCount(b:3); if rate > spaceUsedRateMin
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-direction merge
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 > spaceUsedRate
if rate > spaceUsedRateMin
newBlockTable(newBlockCount, :) = tempBlockMerged;
newBlockCount = newBlockCount + 1; boxNeedCount = boxNeedCount; if rate > spaceUsedRateMin
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-direction merge
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, 4)), blockTable(b, 4)) blockTable(a, 2) + blockTable(b, 2), max(blockTable(a, 3), blockTable(b, 3)), blockTable(b, 3)), 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 > spaceUsedRate
if rate > spaceUsedRateMin
newBlockTable(newBlockCount, :) = tempBlockMerged;
newBlockCount = newBlockCount + 1; boxNeedCount = boxNeedCount; if rate > spaceUsedRateMin
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 end
blockTable = cat(1, blockTable, newBlockTable(1:newBlockCount-1, :));
% Eliminate equivalent complex blocks
blockTableTemp = blockTable(:, 1:5); % Eliminate equivalent complex blocks.
[~, 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 function [rate] = space_used_rate(container, box2, box1)
end
Main function
ts = 1; tf = 0.005; tf = 0.005; tf = 0.005
tf = 0.005; dt = 0.98; dt = 0.98
dt = 0.98; length = 50; dt = 0.98
length = 50; maxSeq = 64; maxSeq = 0; maxSeq = 0
maxSeq = 64; maxChoice = 1024; maxChoice = 1024
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; t = ps; 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 = ps; nps = solution
nps = ps; nps(k) = randi(maxChoice); nps(k)
nps(k) = randi(maxChoice); [solution, nplan, nrate] = container_packeting(container_packeting, nrate)
[solution, nplan, nrate] = container_packeting(container, nps, blockTable, blockNeed);
% disp(['TempPlan: ', mat2str(nplan)]);
% fprintf("TempRate: %f\n",nrate);
click = click + 1; if nrate > rate
if nrate > rate
ps = nps; plan = nplan; if nrate > rate
plan = nplan; rate = nrate; if nrate > rate
rate = nrate.
tempRateX = [tempRateX, click]; if nrate > rate ps = nps; plan = nplan; rate = nrate
tempRate = [tempRate, nrate]; elseif rand(0, 1 > exp(nrate))
elseif rand(0, 1) > exp((nrate - rate) * 10 / t)
ps = nps; plan = nplan; elseif
plan = nplan; rate = nrate; elseif
rate = nrate; end
ps = nps; plan = nplan; rate = nrate; end
if nrate > best_rate
best = nplan; best_rate = nrate; end if nrate > best_rate
best_rate = nrate.
best_solution = solution.
disp(['BestPlan: ', mat2str(best)]);
fprintf("BestRate: %f\n",best_rate);
fprintf("BestRate: \n",best_rate); end
end
t = (1 - t * dt) * t; 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); figure(2)
save Solution best_solution best_rate tempRateX tempRate
