由于要数学建模中需要解决一个三维装箱的问题,我经过搜寻选定了张德富教授等在计算机学报上发表的《求解三维装箱问题的混合模拟退火算法》这篇论文作为解决问题的理论基础。

该论文的摘要如下:

提出了一个高效求解三维装箱问题(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

 

Scroll to Top