Variable bin packing problem with Prolog (CLP)

Question

I am trying to find an algorithm for the NP-hard 2D Variable Size Bin Packing Problem (2DVSBPP) in (Swi-)Prolog using Constraint Logic Programming (CLP).

The problem could be explained like this: some ordered Products need to be packed as efficiently as possible into some Boxes (bins). The products have some given Width and Length (squares or rectangles, eg 2x3). There are four different size of boxes, each with a given cost for the shipper (eg $4 for the 5x5 box, $5 for 5x7 box). The goal is to minimize the total cost from the boxes.

I've been looking for an answer to this problem for a while now and have read numerous papers and similar examples in other languages. However, I can't find any working solution. I'm especially struggling with how to handle the unknown number of Boxes (bins).

---

To be able to find a solution to this problem I've tried to adapt a similar problem but really have no idea how to handle the variable amount of boxes. The following code can choose the cheapest possible box to fit all the products as long as there is only one box needed to fit them all. From the moment we need multiple boxes, the program just fails.

The boxes and products:

:- use_module(library(clpfd)).
:- use_module(library(clpr)).
:- expects_dialect(sicstus).


%% These are the possible productsizes that could need packing
% product (id, width, length)
product(1, 2, 2). 
product(2, 1, 2). 
product(2, 2, 1). % repeating product n2 because it can lay horizontal or vertical
product(3, 1, 3). 
product(3, 3, 1). % idem
product(4, 3, 3). % is square so does not need it
product(5, 2, 3). 
product(5, 3, 2). % iden
product(6, 4, 2). 
product(6, 2, 4). % idem

% because it can lay virtically or horizontally in a box
product_either_way(Number, Width, Length) :-
    product(Number, Width, Length).
product_either_way(Number, Width, Length) :-
    product(Number, Length, Width).


%% These are the so called bins from the 2DVSBPP problem
%% There are 4 sizes, but there is an unlimited supply
% box(Width, Length, Cost)
box(4,4,4).
box(4,6,6).
box(5,5,7).
box(9,9,9).

The constraints:

area_box_pos_combined(W_total*H_total,prod(N),X+Y,f(X,Width,Y,Height)) :-
    product_either_way(N, Width, Height), % Getting the width and height (length) of a product
    % Constraint: the product should 'fit' inside the choosen box
    % thus limiting its coordinates (XY)
    X #>= 1,
    X #=< W_total-Width+1,
    Y #>= 1,
    Y #=< H_total-Height+1.

positions_vars([],[]).
positions_vars([X+Y|XYs],[X,Y|Zs]) :-
    positions_vars(XYs,Zs).

area_boxes_positions_(ProductList,Ps,Zs) :-
    box(W, H, Cost), % finding a suitable box with a W & H
    %% minimize(Cost),
    maplist(area_box_pos_combined(W*H),ProductList,Ps,Cs), % Setting up constraints for each product
    disjoint2(Cs), % making sure they dont overlap with other product inside the box
    positions_vars(Ps,Zs).

A possible query which asks to pack 4 products (numbers 2, 1, 3 and 5)

area_boxes_positions_([prod(2),prod(1),prod(3),prod(5)],Positions,Zs),
labeling([ffc],Zs).

Gives the following as output, one possible way to pack the products:
Positions = [3+1, 1+1, 4+1, 1+3],
Zs = [3, 1, 1, 1, 4, 1, 1, 3] .

But how do I model multiple boxes, when we would have an order with more products that would not fit inside one box?

Any help or examples are really appreciated!

Answer 1

There are some redundancies in your partial solution, maybe caused by premature optimization.

First, since you have a product_either_way/3, you should not change your input specification, adding products with same id and dimensions swapped. After all, width and height are properties you cannot swap arbitrarly in the real world, and you already have produced a predicate that takes care of this, so I have started removing such duplicates.

Second, the purpose of disjoint/2 is to compute a placement of a set of rectangles, so your area_box_pos_combined/4 and positions_vars/2 are pretty much useless.

Here is how I would approach this problem. First, write a predicate that given a list of products and a box, puts as many as possible into it, and 'returns' those that didn't fit. For instance

fill_box([P|Ps],W,H,Placed,Rs) :-
    (   product(P,W_i,H_i)
    ;   product(P,H_i,W_i)
    ),
    W_p #= W - W_i,
    H_p #= H - H_i,
    X_i in 0..W_p,
    Y_i in 0..H_p,
    U=[p(X_i, W_i, Y_i, H_i)|Placed],
    disjoint2(U),
    fill_box(Ps,W,H,U,Rs).
fill_box(Rs,_,_,_,Rs).

It's somewhat buggy, because it will stop at the first product it cannot place, but there could be more placeable after this. But what's important, now we can start to test if it's working, given the interaction with key concepts of CLP(FD). disjoint/2 works on bounded variables, so the domain declaration of X_i and Y_i is needed.

?- fill_box([1,1],4,2,[],R).
R = [] .

?- fill_box([1,1],3,2,[],R).
R = [1] .

Now we can provide a driver, maybe as simple as

products_placed_cost([],0).
products_placed_cost(Ps,C) :-
    box(W,H,C0),
    fill_box(Ps,W,H,[],Rs),
    Ps\=Rs,
    products_placed_cost(Rs,C1),
    C #= C0+C1.

and then let Prolog generate as many solutions as it can, just order them by cost, by means of library(solution_sequences):

?- order_by([asc(C)],products_placed_cost([1,1],C)).
C = 4 ;
C = 4 ;
C = 4 ;
C = 4 ;
C = 6 ;
...

But we don't know which placements have been generated. We have to add arguments that carry back the information. Then

products_placed_cost([],[],0).
products_placed_cost(Ps,[box(W,H,C0,Q)|Qs],C) :-
    box(W,H,C0),
    fill_box(Ps,W,H,[],Rs,Q),
    Ps\=Rs,
    products_placed_cost(Rs,Qs,C1),
    C #= C0+C1.

fill_box([P|Ps],W,H,Placed,Rs,[P|Qs]) :-
    (   product(P,W_i,H_i)
    ;   product(P,H_i,W_i)
    ),
    W_p #= W - W_i,
    H_p #= H - H_i,
    X_i in 0..W_p,
    Y_i in 0..H_p,
    U=[p(X_i, W_i, Y_i, H_i)|Placed],
    disjoint2(U),
    fill_box(Ps,W,H,U,Rs,Qs).
fill_box(Rs,_,_,_,Rs,[]).

To be true, library(clpfd) is used just as commodity, but mixed with the searching capabilities of (pure) Prolog gives us a short and declarative solution.

See the specific documentation of library(clpBNR) for a better approach.

Answer 2

I'm especially struggling with how to handle the unknown number of Boxes (bins).

You can put an upper bound on the number of boxes: For N indivisible elements you will never need more than N boxes. Furthermore, we can define a special "unused" kind of box with 0 size but 0 cost. Then we can ask for a solution with an assignment of items to exactly N (or any other number of) boxes, some of which can remain unused.

Here is a description of a single box, relating its kind, size, and cost using disjunctive and conjunctive constraints:

kind_width_length_cost(Kind, Width, Length, Cost) :-
    % unused box
    (Kind #= 0 #/\ Width #= 0 #/\ Length #= 0 #/\ Cost #= 0) #\/
    % small box
    (Kind #= 1 #/\ Width #= 4 #/\ Length #= 4 #/\ Cost #= 4) #\/
    % medium box
    (Kind #= 2 #/\ Width #= 4 #/\ Length #= 6 #/\ Cost #= 6) #\/
    % large box
    (Kind #= 3 #/\ Width #= 5 #/\ Length #= 5 #/\ Cost #= 7) #\/
    % X-large box
    (Kind #= 4 #/\ Width #= 9 #/\ Length #= 9 #/\ Cost #= 9),
    % make sure all variables have finite domains, the above disjunction is
    % not enough for the system to infer this
    Kind in 0..4,
    Width in 0..9,
    Length in 0..9,
    Cost in 0..9.

A collection of N boxes can be represented as a term boxes(Numbers, Kinds, Widths, Lengths, Costs) where Numbers are [1, 2, ..., N] and the I-th element of each of the other lists is the length/width/cost of box number I:

n_boxes(N, boxes(Numbers, Kinds, Widths, Lengths, Costs)) :-
    numlist(1, N, Numbers),
    length(Kinds, N),
    maplist(kind_width_length_cost, Kinds, Widths, Lengths, Costs).

For example, three boxes are:

?- n_boxes(3, Boxes).
Boxes = boxes([1, 2, 3], [_G9202, _G9205, _G9208], [_G9211, _G9214, _G9217], [_G9220, _G9223, _G9226], [_G9229, _G9232, _G9235]),
_G9202 in 0..4,
_G9202#=4#<==>_G9257,
_G9202#=3#<==>_G9269,
_G9202#=2#<==>_G9281,
_G9202#=1#<==>_G9293,
_G9202#=0#<==>_G9305,
... a lot more constraints

Note that this uses a term containing lists rather than the more "usual" representation with a list containing terms box(Num, Width, Length, Cost). The reason for this is that we will want to index into these lists of FD variables using element/3. This predicate cannot be used to index into lists of other terms.

Turning to products, here is the FD version of your disjunctive product_either_way predicate:

product_either_way_fd(Number, Width, Length) :-
    product_width_length(Number, W, L),
    (Width #= W #/\ Length #= L) #\/ (Width #= L #/\ Length #= W),
    % make sure Width and Length have finite domains
    Width #>= min(W, L),
    Width #=< max(W, L),
    Length #>= min(W, L),
    Length #=< max(W, L).

The placement of an item is expressed with a term box_x_y_w_l containing the number of the box, the X and Y coordinates inside the box, and the item's width and length. The placement must be compatible with the dimensions of the chosen box:

product_placement(Widths, Lengths, Number, Placement) :-
    product_either_way_fd(Number, W, L),
    Placement = box_x_y_w_l(_Box, _X, _Y, W, L),
    placement(Widths, Lengths, Placement).

placement(Widths, Lengths, box_x_y_w_l(Box, X, Y, W, L)) :-
    X #>= 0,
    X + W #=< Width,
    Y #>= 0,
    Y + L #=< Length, 
    element(Box, Widths, Width),
    element(Box, Lengths, Length).

This is where we use the Widths and Lengths lists of FD variables. The number of the chosen box is itself an FD variable that we use as an index to look up the box's width and length using the element/3 constraint.

Now we must model non-overlapping placements. Two items placed in different boxes are automatically non-overlapping. For two items in the same box we must check their coordinates and sizes. This binary relation must be applied to all unordered pairs of items:

placement_disjoint(box_x_y_w_l(Box1, X1, Y1, W1, L1),
                   box_x_y_w_l(Box2, X2, Y2, W2, L2)) :-
    Box1 #\= Box2 #\/
    (Box1 #= Box2 #/\
     (X1 #>= X2 + W2 #\/ X1 + W1 #< X2) #/\
     (Y1 #>= Y2 + L2 #\/ Y1 + L1 #< Y2)).

alldisjoint([]).   
alldisjoint([Placement | Placements]) :-
    maplist(placement_disjoint(Placement), Placements),
    alldisjoint(Placements).

Now we are ready to put everything together. Given a list of products and a number N of boxes (some of which may be unused), the following predicate computes constraints on placements in boxes, the kinds of boxes used, their costs, and a total cost:

placements_(Products, N, Placements, BoxKinds, Costs, Cost) :-
    n_boxes(N, boxes(_BoxNumbers, BoxKinds, Widths, Lengths, Costs)),
    maplist(product_placement(Widths, Lengths), Products, Placements),
    alldisjoint(Placements),
    sum(Costs, #=, Cost).

This constructs a term representing N boxes, computes placement constraints for each product, ensures the placements are disjoint, and sets up the computation of the total cost. That is all!

I'm using the following products copied from the question. Note that I have removed duplicates with swapped widths/lengths since this swapping is done by product_either_way_fd when needed.

product_width_length(1, 2, 2).
product_width_length(2, 1, 2).
product_width_length(3, 1, 3).
product_width_length(4, 3, 3).
product_width_length(5, 2, 3).
product_width_length(6, 4, 2).

We're ready for testing. To reproduce your example of placing items 2, 1, 3, and 5 in a single box:

?- placements_([2, 1, 3, 5], 1, Placements, Kinds, Costs, Cost).
Placements = [box_x_y_w_l(1, _G17524, _G17525, _G17526, _G17527), box_x_y_w_l(1, _G17533, _G17534, 2, 2), box_x_y_w_l(1, _G17542, _G17543, _G17544, _G17545), box_x_y_w_l(1, _G17551, _G17552, _G17553, _G17554)],
Kinds = [_G17562],
Costs = [Cost],
_G17524 in 0..8,
_G17524+_G17526#=_G17599,
_G17524+_G17526#=_G17611,
_G17524+_G17526#=_G17623,
...

With labeling:

?- placements_([2, 1, 3, 5], 1, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), labeling([], Variables).
Placements = [box_x_y_w_l(1, 0, 0, 1, 2), box_x_y_w_l(1, 7, 7, 2, 2), box_x_y_w_l(1, 4, 6, 3, 1), box_x_y_w_l(1, 2, 3, 2, 3)],
Kinds = [4],
Costs = [9],
Cost = 9,
Variables = [0, 0, 1, 2, 7, 7, 4, 6, 3|...] .

(You might want to check this carefully for correctness!) Everything was placed in box number 1, which is of kind 4 (size 9x9) with cost 9.

Is there a way to fit these items in a cheaper box?

?- Cost #< 9, placements_([2, 1, 3, 5], 1, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), labeling([], Variables).
false.

Now, how about putting all the products in (up to) 6 boxes?

?- placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), labeling([], Variables).
Placements = [box_x_y_w_l(1, 0, 0, 2, 2), box_x_y_w_l(1, 3, 3, 1, 2), box_x_y_w_l(1, 5, 6, 1, 3), box_x_y_w_l(2, 0, 0, 3, 3), box_x_y_w_l(2, 4, 4, 2, 3), box_x_y_w_l(3, 0, 0, 2, 4)],
Kinds = [4, 4, 1, 0, 0, 0],
Costs = [9, 9, 4, 0, 0, 0],
Cost = 22,
Variables = [1, 0, 0, 1, 3, 3, 1, 2, 1|...] .

The first solution found uses three boxes and left the other three unused. Can we go cheaper?

?- Cost #< 22, placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), labeling([], Variables).
Cost = 21,
Placements = [box_x_y_w_l(1, 0, 0, 2, 2), box_x_y_w_l(1, 3, 3, 1, 2), box_x_y_w_l(1, 5, 6, 1, 3), box_x_y_w_l(2, 0, 0, 3, 3), box_x_y_w_l(3, 0, 0, 2, 3), box_x_y_w_l(4, 0, 0, 2, 4)],
Kinds = [4, 1, 1, 1, 0, 0],
Costs = [9, 4, 4, 4, 0, 0],
Variables = [1, 0, 0, 1, 3, 3, 1, 2, 1|...] .

Yes! This solution uses more boxes, but ones that are overall slightly cheaper. Can we do even better?

?- Cost #< 21, placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), labeling([], Variables).
% ... takes far too long

We need to be a bit more sophisticated. Playing around with the number of boxes it's clear that cheaper solutions with fewer boxes are available:

?- Cost #< 21, placements_([1, 2, 3, 4, 5, 6], 2, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), labeling([], Variables).
Cost = 18,
Placements = [box_x_y_w_l(1, 0, 0, 2, 2), box_x_y_w_l(1, 3, 3, 1, 2), box_x_y_w_l(1, 5, 6, 1, 3), box_x_y_w_l(2, 0, 6, 3, 3), box_x_y_w_l(2, 6, 4, 3, 2), box_x_y_w_l(2, 4, 0, 2, 4)],
Kinds = [4, 4],
Costs = [9, 9],
Variables = [1, 0, 0, 1, 3, 3, 1, 2, 1|...] .

Maybe directing the search to label box kinds first is useful, since the up strategy will essentially try to use as few boxes as possible:

?- Cost #< 21, placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), time(( labeling([], Kinds), labeling([ff], Variables) )).
% 35,031,786 inferences, 2.585 CPU in 2.585 seconds (100% CPU, 13550491 Lips)
Cost = 15,
Placements = [box_x_y_w_l(5, 2, 4, 2, 2), box_x_y_w_l(6, 8, 7, 1, 2), box_x_y_w_l(6, 5, 6, 3, 1), box_x_y_w_l(6, 2, 3, 3, 3), box_x_y_w_l(6, 0, 0, 2, 3), box_x_y_w_l(5, 0, 0, 2, 4)],
Kinds = [0, 0, 0, 0, 2, 4],
Costs = [0, 0, 0, 0, 6, 9],
Variables = [5, 2, 4, 6, 8, 7, 1, 2, 6|...] .

This really does need ff or ffc, the default leftmost strategy doesn't return results in a reasonable time frame.

Can we do even better?

?- Cost #< 15, placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), time(( labeling([], Kinds), labeling([ff], Variables) )).
% 946,355,675 inferences, 69.984 CPU in 69.981 seconds (100% CPU, 13522408 Lips)
false.

No! The solution with cost 15 is optimal (but not unique).

However, I find 70 seconds to be too slow for this very small problem size. Are there some some symmetries we can exploit? Consider:

?- Cost #= 15, placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), time(( labeling([], Kinds), labeling([ff], Variables) )).
% 8,651,030 inferences, 0.611 CPU in 0.611 seconds (100% CPU, 14163879 Lips)
Cost = 15,
Placements = [box_x_y_w_l(5, 2, 4, 2, 2), box_x_y_w_l(6, 8, 7, 1, 2), box_x_y_w_l(6, 5, 6, 3, 1), box_x_y_w_l(6, 2, 3, 3, 3), box_x_y_w_l(6, 0, 0, 2, 3), box_x_y_w_l(5, 0, 0, 2, 4)],
Kinds = [0, 0, 0, 0, 2, 4],
Costs = [0, 0, 0, 0, 6, 9],
Variables = [5, 2, 4, 6, 8, 7, 1, 2, 6|...] .

?- Kinds = [4, 2, 0, 0, 0, 0], Cost #= 15, placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), time(( labeling([], Kinds), labeling([ff], Variables) )).
% 11,182,689 inferences, 0.790 CPU in 0.790 seconds (100% CPU, 14153341 Lips)
Kinds = [4, 2, 0, 0, 0, 0],
Cost = 15,
Placements = [box_x_y_w_l(1, 7, 7, 2, 2), box_x_y_w_l(1, 6, 5, 1, 2), box_x_y_w_l(2, 3, 3, 1, 3), box_x_y_w_l(2, 0, 0, 3, 3), box_x_y_w_l(1, 4, 2, 2, 3), box_x_y_w_l(1, 0, 0, 4, 2)],
Costs = [9, 6, 0, 0, 0, 0],
Variables = [1, 7, 7, 1, 6, 5, 1, 2, 2|...] .

These aren't permutations of the same solution, but they are permutations of the same boxes and therefore have identical costs. We don't need to consider both of them! In addition to labeling Kinds a bit more intelligently than in the beginning, we can also require the Kinds list to be monotonically increasing. This excludes lots of redundant solutions and gives much faster termination, and even with better solutions first:

?- placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), chain(Kinds, #=<), time(( labeling([], Kinds), labeling([ff], Variables) )).
% 34,943,765 inferences, 2.865 CPU in 2.865 seconds (100% CPU, 12195550 Lips)
Placements = [box_x_y_w_l(5, 2, 4, 2, 2), box_x_y_w_l(6, 8, 7, 1, 2), box_x_y_w_l(6, 5, 6, 3, 1), box_x_y_w_l(6, 2, 3, 3, 3), box_x_y_w_l(6, 0, 0, 2, 3), box_x_y_w_l(5, 0, 0, 2, 4)],
Kinds = [0, 0, 0, 0, 2, 4],
Costs = [0, 0, 0, 0, 6, 9],
Cost = 15,
Variables = [5, 2, 4, 6, 8, 7, 1, 2, 6|...] .

?- Cost #< 15, placements_([1, 2, 3, 4, 5, 6], 6, Placements, Kinds, Costs, Cost), term_variables(Placements, Variables, [Cost | Costs]), chain(Kinds, #=<), time(( labeling([], Kinds), labeling([ff], Variables) )).
% 31,360,608 inferences, 2.309 CPU in 2.309 seconds (100% CPU, 13581762 Lips)
false.

More tweaks are possible and probably necessary for larger problem sizes. I found that adding bisect in the final labeling helps a bit. So does removing the logically redundant Box1 #= Box2 constraint in placement_disjoint/2. Finally, given the use of chain/2 to restrict Kinds, we can remove the preliminary labeling of Kinds entirely to get a nice speedup! I'm sure there's more, but for a prototype I think it's reasonable enough.

Thank you for this interesting problem!