The Partridge Packing Problem is a packing puzzle that was originally proposed by Robert T. Wainwright at G4G2 (the Second Gathering for Gardner conference) in 1996. In this post we will model and solve the Partridge Packing Problem using MiniZinc. The inspiration was Matt Parker’s fun video on the problem.

Packing problems are a classic use-case for combinatorial solvers. In fact, the original paper that introduced the idea of global constraints for constraint programming, “Introducing global constraints in CHIP” by Beldiceanu and Contejean 1994 included the so-called diffn constraint for packing problems. The constraint ensures that a set of (n-dimensional) boxes are not overlapping.¹

This post assumes some familiarity with MiniZinc. For some background on MiniZinc, see the previous posts in the collection. The puzzle will be explained fully, and no specific knowledge of packing problems is assumed.

The Problem#

The Partridge Packing Problem is a packing problem for squares in a larger square. For size $n$, the goal is to pack:

• 1 square of size $1 \times 1$
• 2 squares of size $2 \times 2$
• …
• $n$ squares of size $n \times n$

into a square of size $\frac{n(n+1)}{2} \times \frac{n(n+1)}{2}$.² The name comes from the song “The Twelve Days of Christmas,” where the first gift is a partridge in a pear tree, then two turtle doves, and so on going up to twelve drummers drumming.

The sum of the areas of all the smaller squares equals the area of the larger square

$\sum_{i=1}^{n} i \cdot i^2 = \frac{n(n+1)}{2} \cdot \frac{n(n+1)}{2}$

But just because the area matches does not mean that it is possible. It is known that sizes 2 to 7 have no solution, and sizes from 8 to 33 have at least one solution. The problem becomes increasingly difficult as $n$ grows larger, as the number of parts grows quadratically.

Let’s look at the first interesting size with a solution, size 8. Here are all the parts to pack.³

These parts can be packed in a square of size $36\times 36$, where $36$ comes from $\frac{8 \times 9}{2} = 36$, and here is one such solution.

This visualization shows how all the squares pack together perfectly to fill the 36×36 grid.

As mentioned, for sizes below 8 the problem is infeasible (except 1, which is the trivial case). Consider size 2, which includes 1 part of size $1\times 1$ and 2 parts of size $2\times 2$ that should be packed in a $3\times 3$ square. As can be seen below, while the sum of the areas of the parts equals the area to pack in, there is no way to put the two larger squares on the area without them overlapping.

The Base MiniZinc Model#

Following previous parts in this collection, we will split up the model in parts. In this section the first basic model will be presented, including the data, the viewpoint, the basic constraints, and the search and output.

In the next section, improvements to the model will be discussed. Several of the improvements were suggested by Mats Carlsson, and made the model a lot better and faster.

Representing the Problem#

The problem is parameterized by a single value $n$, which determines both the number of different square sizes and the size of the target square.

int: n;
set of int: N = 1..n;

% Triangular number of n is both the total number of parts and
% the board size length
int: triangular_n = (n * (n+1)) div 2;

enum Parts = P(1..triangular_n);
set of int: Pos = 1..triangular_n;

array[Parts] of N: sizes = array1d(Parts, reverse([
    size
  | size in N, copy in 1..size
]));

constraint assert(sum(s in sizes) (s * s) == triangular_n * triangular_n,
                  "The squares fill the board completely");

The computed value triangular_n is the triangular number of the size parameter n. This is both the total number of parts to pack as well as the side length of the board where the parts are to be placed. The enum Parts is used to separate the set of parts from the Pos positions to place them at.⁴

The sizes are generated in increasing order but are reversed, resulting in the larger boxes being first in the array. This is useful since many solvers will use the input-order as a tie-breaker for heuristics, promoting packing hard-to-pack boxes (i.e., the larger ones) first.

Similar to the LinkedIn Queens post, we can use instance files to set the parameter n. However, running the model from the MiniZinc IDE the user is prompted for all unknown values, and for a single integer this is very easy to supply.

The Viewpoint#

There are many ways that one can model a packing problem. The most common way for box packings is to set one corner as the reference point, and to use the position of that reference point as the position for the box. The most natural expression for this is to use two arrays representing the x and y coordinates of the bottom-left corner of each square.

% Main variables for placement of parts, the x and y coordinate
array[Parts] of var Pos: x;
array[Parts] of var Pos: y;

MiniZinc has a feature where records can be used to structure data, and using that, we could declare the variables like this instead.

% Main variables for placement of squares, the x and y coordinate
array[Parts] of record(var Pos: x, var Pos: y): box;

However, there are several places in the model where a constraint is formulated over the x variables only, and then over the y variables. Therefore, it is easier to use two arrays instead of a single one.⁵

Basic Packing Constraints#

The base variables allow placement of the reference point anywhere inside the packing area. However, the allowed positions need to be adjusted based on the size of a part. This is done by adjusting the upper bounds of the x and y value based on the size, ensuring that the point is also in the Pos set.

constraint :: "Parts fit in x direction"
forall(p in Parts) (
    x[p] + sizes[p] - 1 in Pos
);

constraint :: "Parts fit in y direction"
forall(p in Parts) (
    y[p] + sizes[p] - 1 in Pos
);

In the above (and the rest of the constraint here), constraints are named using the :: string annotation. These names, such as "Parts fit in x direction", are translated into the FlatZinc format and are useful for debugging and for tools such as findMUS.

The main constraint for a packing problem is that no parts should overlap. The classic way to ensure this is to use the no-overlap constraint, which for historic reasons is named the diffn constraint in MiniZinc.

constraint :: "No-overlap packing constraint"
diffn(x, y, sizes, sizes);

The arguments to diffn are the x and y positions of the rectangles, and their extent in the x and y direction (that is, the width and the height). Since the parts are squares, their extents are the same in both directions.

Search and output#

This is a satisfaction problem and we will leave the search strategy to the solver.

solve satisfy;

There are two output blocks for this model. The first block will print an ASCII-art representation of the packing to the standard output.

/**
* Get the unique singleton value in the supplied set, assert if it is not a singleton.
*/
function $$T: singleton_value(set of $$T: values) =
    assert(card(values) == 1, "Values must have exactly one element, was \(values)",
        min(values)
    );

/**
* Return a character representation of the value v.
*
* Support values v in 1..36.
*/
function string: to_char(int: v) =
    if v in 0..9 then
        "\(v)"
    else
        ["a", "b", "c", "d", "e", "f", "g", "h",
        "i", "j", "k", "l", "m", "n", "o", "p",
        "q", "r", "s", "t", "u", "v", "w",
        "x", "y", "z"][v-9]
    endif;

% Base command-line output mapping the placed parts to their sizes.
%
output [
    let {
        any: fx = fix(x),
        any: fy = fix(y),
        any: board = array2d(Pos, Pos, [
            let {
                Parts: part_id = singleton_value({p | p in Parts where
                    tx in fx[p]..(fx[p] + sizes[p]-1) /\
                    ty in fy[p]..(fy[p] + sizes[p]-1)
                })
            } in
            to_char(sizes[part_id])
          | tx in Pos, ty in Pos
        ])
    } in
    concat(tx in Pos) (
        concat(board[tx, ..]) ++ "\n"
    )
];

While long, this code is reasonably straightforward. First, there are two helper functions: singleton_value, which transforms a set that is known to be just one element to the element, and to_char, which transforms a size to a character that represents it in base 36 (0-9 and a-z).

Next, a matrix is constructed where for each position, the part that is covering that position is found, and the size of that part is used to get the character. Finally, this matrix is concatenated into a set of strings.

The second output-block uses a feature of the MiniZinc IDE where custom visualizations can be used. These work by starting a webserver serving a webpage that receives the solutions as they are produced. For this problem, the existing vis_geost_2d visualization is used.

output vis_geost_2d(
    % Internal x and y offset of each part, 0 since each part is its own shape
  [p:0 | p in Parts], [p:0 | p in Parts],
  % Size of each part in x and y direction
  sizes, sizes,
  % Map each shape to the corresponding single part
  [p:{p} | p in Parts],
  % Reference points for each shape
  x, y,
  % The kind of each part
  array1d(Parts, Parts)
);

The vis_geost_2d family of visualizations can show packing problems with shapes made out of rectangles using internal offsets to a common shape reference point, matching the input for the geost constraint. As each part is just a square, each kind of shape will be a single part, and the internal offsets are just 0. Note that the construction [p:0 | p in Parts] will create an array with Parts as the index set, skipping the p: part would create an array with 1..card(Parts) as the index set. An alternative way to write this is to coerce the base array to the right index set: array1d(Parts, [0 | p in Parts]).

Performance of base model#

In all the tests here, we will use OR-Tools CP-SAT 9.14 bundled with MiniZinc IDE 2.9.4 on a MacBook Pro M1 Max with 64 GiB of memory. The configuration is set to use 10 threads (same as the number of cores in the CPU), and use free search.

As mentioned, sizes 2 to 7 are unsatisfiable, so the smallest interesting problem with a solution is size 8. However, this base model is not efficient at all. Finding a solution took about 3 and a half hours in one run, which makes it not very practical.

777777744448888888888888888333666666
777777744448888888888888888333666666
777777744448888888888888888333666666
777777744448888888888888888333666666
777777744448888888888888888333666666
777777744448888888888888888333666666
777777744448888888888888888227777777
444433344448888888888888888227777777
444433322777777777777776666667777777
444433322777777777777776666667777777
444455555777777777777776666667777777
444455555777777777777776666667777777
444455555777777777777776666667777777
444455555777777777777776666667777777
444455555777777777777776666667777777
777777788888888888888886666667777777
777777788888888888888886666667777777
777777788888888888888886666667777777
777777788888888888888886666667777777
777777788888888888888886666667777777
777777788888888888888885555588888888
777777788888888888888885555588888888
666666188888888888888885555588888888
666666555555555577777775555588888888
666666555555555577777775555588888888
666666555555555577777775555588888888
666666555555555577777775555588888888
666666555555555577777775555588888888
888888888888888877777775555588888888
888888888888888877777775555588888888
888888888888888866666666666688888888
888888888888888866666666666688888888
888888888888888866666666666688888888
888888888888888866666666666688888888
888888888888888866666666666688888888
888888888888888866666666666688888888
----------
==========
%%%mzn-stat: nSolutions=1
%%%mzn-stat-end
%%%mzn-stat: boolVariables=1023
%%%mzn-stat: failures=88389736
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: propagations=3870695549
%%%mzn-stat: solveTime=12697.9
%%%mzn-stat-end
Finished in 3h 31m 38s.

While the ASCII art is nice, the visualization is much easier to understand. Below you can see first the visualization from MiniZinc, and then the visualization for this post where squares of equal size get the same color and all squares are marked with their size.

Improved model#

The above model is the base, with just the constraints that are needed for a correct solution. In this part, we will add additional constraints that improve the model significantly. These constraints are of two types, implied constraints and symmetry breaking constraints. An implied constraint is a constraint that strengthens the model by adding additional constraints that are true in every solution. The goal is to add additional propagation that makes more deductions. A symmetry breaking constraint is used to reduce the number of solutions, by limiting the symmetries of solutions.

Symmetries often arise from modeling decisions, but sometimes also from the problem itself. For example, in the classic 8-queens problem there is a symmetry from the problem definition: the chessboard for a single solution can be rotated and mirrored diagonally to create 8 different solutions. If the model were to name the queens, then that would introduce a symmetry for which queen is placed where. This symmetry would occur because of modeling decisions, not from the problem itself where queens are indistinguishable.⁶

We will use a feature of MiniZinc to mark constraints with their type by enclosing the constraint in calls to implied_constraint and symmetry_breaking_constraint. While not useful for many solvers, some (such as Constraint-Based Local Search solvers) can use this information to decide what constraints to soften and what constraints to use for moves.

For each improvement, we will test it to see the effects. Note that the configuration that is used, OR-Tools CP-SAT with 10 threads, is not a deterministic system. One single run might not be indicative for all runs, but in most cases it will be a good indication.

Cumulative profile#

A classic implied constraint for packing problems is to add a cumulative profile constraint for the x and y direction. Cumulative is a classic scheduling constraint, and is typically used for tasks that use some set of resources while they are active. Below is an example of 8 tasks that are scheduled, with a capacity limit of 8 and varying amounts of usage at different points.

Note that the tasks do not have a fixed y-position; they only have a start, an end, and a resource usage (height). This means that tasks like the green task 4 and the purple task 6 are not shown as a rectangle but staggered based on the amount of other tasks. For the packing case, looking along one dimension, the orthogonal dimension can be seen as a resource, and the squares as tasks to be scheduled. This is a classic implied constraint that can strengthen the propagation, and OR-Tools CP-SAT even has several parameters that can be set to include cumulative-style reasoning. Here, the cumulative constraint is instead added as a MiniZinc constraint so that it can be used with all different solvers.

constraint :: "Cumulative profile of parts along the x axis." implied_constraint(
cumulative(x, sizes, sizes, card(Pos))
);
constraint :: "Cumulative profile of parts along the y axis." implied_constraint(
cumulative(y, sizes, sizes, card(Pos))
);

Running this, however, does not give better results at all. The simple model took three and a half hours, but this model takes more than an hour more!

%%%mzn-stat: boolVariables=2184
%%%mzn-stat: failures=99470613
5 collapsed lines
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: propagations=7359764734
%%%mzn-stat: solveTime=17031.9
%%%mzn-stat-end
Finished in 4h 43m 52s.

Unfortunately, this type of behavior is not uncommon when a learning system with automatic heuristics and randomization is combined with changes to a model. This shows the importance of benchmarking and testing all changes to see how the model behaves. Even well-known improvements might make it worse.

Exact fill#

The cumulative constraint above adds to the reasoning, but it is also a lot weaker than it could have been. The Partridge Packing Problem is a tight packing, where the board is fully covered. The cumulative constraint “just” says that too much area can’t be used. Consider instead a constraint that, for each row and column, checks which parts overlap it and requires that the sum of the sizes of overlapping parts equals the board size exactly.

% The sizes of the parts that overlap rc in the xy direction
% must equal the number of positions exactly.
predicate exact_fill(array[Parts] of var Pos: xy, Pos: rc) =
    let {
        % on_rc[p] is true iff the part overlaps the row/column rc
        array[Parts] of var bool: on_rc = [
            rc-sizes[p] < xy[p] /\ xy[p] <= rc
          | p in Parts
        ]
    } in
    sum(p in Parts) (
        sizes[p] * on_rc[p]
    ) = card(Pos);

constraint :: "Exact profile of parts along the x axis." implied_constraint(
forall(rc in Pos) (
    exact_fill(x, rc)
)
);

constraint :: "Exact profile of parts along the y axis." implied_constraint(
forall(rc in Pos) (
    exact_fill(y, rc)
)
);

Here, a utility function is added so that the right sum can be constructed for each column and for each row. The exact_fill function takes the positions of all the parts along either the x or y axis, and a specified row or column. Inside, a local array on_rc indexed by Parts of Boolean variables is constructed that indicates whether each part overlaps that row or column. Multiplying by the size of each part gives how much of the dimension is used, and that is required to be equal to the cardinality of the Pos set.

This addition is a huge improvement over the base model! A solution is found in less than 4 minutes instead of 3 and a half hours.

%%%mzn-stat: boolVariables=3960
%%%mzn-stat: failures=6155170
5 collapsed lines
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: propagations=1762225146
%%%mzn-stat: solveTime=225.119
%%%mzn-stat-end
Finished in 3m 45s.

This is starting to look like a viable model to use. Checking if the cumulative constraint might help now shows that it is still not a good addition, and it increased the search time to 4 minutes 33 seconds.

%%%mzn-stat: boolVariables=3960
%%%mzn-stat: failures=7544566
5 collapsed lines
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: propagations=2046103308
%%%mzn-stat: solveTime=272.594
%%%mzn-stat-end
Finished in 4m 33s.

Restrictions of placements close to the edge#

From the work that Mats Carlsson and Nicolas Beldiceanu did creating the geost constraint, there are several additional deductions that can be made based on placements of boxes. The core insight in this case is that since the board should be filled completely, then for every area created there must be parts that can fill it. Consider the below packing where a part has been placed on the board close to the edge.

The red area next to the border has a width of 2 and a height of 6. It can only be packed with parts that are at most size 2, and a total area of $2\cdot 6=12$ needs to be available. However, for parts up to size 2, this is not possible since there is one $1\times 1$ square and two $2\times 2$ squares, for a total area of 9. Trying to fill up the area between the size 6 part and the border would look like this.

Given the above reasoning, it is clear that any part of size 6 must either be placed next to a border, or at a distance of more than 2 from a border. In general, for a given size $n$, the sum of the areas of the smaller parts (up to size $n-1$) is the square of the triangular number for $n-1$. This reasoning can be generalized and implemented with the following MiniZinc code.

% The amount of available area from parts up to given size
function int: available_area(int: size) =
  let {
      % t is the triangular number of size
      int: t = (size * (size + 1)) div 2;
  } in
  t * t;

constraint :: "Edge-placement limits" implied_constraint(
forall(size in N where size > 1) (
  let {
      % Find the smallest distance from the edge that is possible to place.
      int: min_distance_from_edge = min({d | d in 1..size
                where d * size > available_area(d)}),
      % Placing in these positions is not packable for a full packing
      set of int: forbidden_placements =
          % Positions at low placement indices
          2..(1+min_distance_from_edge)
           union
          % positions at high placement indices
          max(Pos)-size-min_distance_from_edge..<max(Pos)-size,
      set of Pos: allowed_placements = Pos diff forbidden_placements
  } in
  forall(p in Parts where sizes[p] = size) (
       x[p] in allowed_placements
        /\
       y[p] in allowed_placements
   )
));

For each size of part, there is a custom calculation of the allowed_placements for that part. Since the parts and the board are squares, the same set can be used for both x and y placements. The calculation of min_distance_from_edge uses the idea that if the part is d steps away from the edge, then the available area of parts up to size d must be greater than the size length times the value d for it to be valid. Using this, the set of forbidden_placements is computed close to the edges, and the allowed_placements are the complement of that with respect to Pos. This is a conservative approximation of the packability: if this requirement is not satisfied, then there is no packing that would work.

Adding this constraint reduces the time significantly again. Running three times the time varies between 45 and 105 seconds due to the stochastic nature of the solving process. The median has the following statistics.

2 collapsed lines
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: boolVariables=3812
%%%mzn-stat: failures=2251053
%%%mzn-stat: propagations=558765030
3 collapsed lines
%%%mzn-stat: solveTime=73.3695
%%%mzn-stat: nSolutions=1
%%%mzn-stat-end
Finished in 1m 13s.

In the original geost work, this type of reasoning is not just limited to placements close to an edge, but for all different types of induced areas during the search. This is much stronger reasoning, but would not be readily expressible as fixed constraints. It requires careful implementation as a propagator in a system. SICStus Prolog has the original and probably most advanced implementation of geost with a large domain-specific language to express placements.

Symmetry breaking for same size parts#

Symmetry breaking is often crucial in many problems. Here, the focus is a symmetry that is introduced by the modeling: parts of the same size should be indistinguishable. The three parts of size $3\times 3$ are the same, but since they have different identifiers they are different to the solvers. A common way to break symmetries is to introduce an ordering among the alternatives.

constraint :: "Equal size squares symmetry breaking" symmetry_breaking_constraint(
forall (size in N) (
   let {
       set of Parts: PartsWithSize = {p | p in Parts where sizes[p] == size},
       set of int: P = 1..card(PartsWithSize),
       array[1..2, P] of var Pos: placements = array2d(1..2, P, [
           [x[p], y[p]][x_or_y]
         | x_or_y in 1..2, p in PartsWithSize
       ])
   } in
   lex_chain_less(placements)
)
);

For each size, the set of parts with that size are collected. Then, a matrix of placements is constructed where each column represents the x and y coordinates of a part for that size.⁷ The lex_chain_less constraint is used to order these tuples using lexicographic ordering.

Adding the symmetry reduces the solving time significantly again. In 10 runs, it is between 0.8 and 3.6 seconds, with an average of 1.9 seconds. The median has the following statistics.

2 collapsed lines
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: boolVariables=3700
%%%mzn-stat: failures=15018
%%%mzn-stat: propagations=8614720
3 collapsed lines
%%%mzn-stat: solveTime=1.4264
%%%mzn-stat: nSolutions=1
%%%mzn-stat-end
Finished in 1s 830msec.

How about the board symmetry?#

As mentioned above, the board has 8 symmetries (four rotations times flipping), and it is common to break them in many puzzle cases. Matt Parker argues in the video that for the purposes of this puzzle, they should be kept in. Also, it can be quite tricky to combine symmetry breaking techniques. For any way to order the symmetries of the board, that ordering would have to work jointly with the ordering of the parts.

The Full Program#

For testing, you can download the full MiniZinc model. Remember to set OR-Tools CP-SAT to use at least as many threads as you have cores, and to also check the free search box.

Solving larger instances#

In all the above examples, size 8 has been the instance solved. Using the model developed, let’s try larger sizes and see the performance for that.

Size 9#

In Matt Parker’s video that inspired this post, size 9 was the instance that was discussed. This is because size 9 has a side-length of 45, and thus the area of the board is $45^2=2025$, which is the year the video was published.

Remember, even though the step from 8 to 9 sounds small, the number of parts grows from 36 to 45. In a couple of tests, it took between 61 and 86 seconds to solve size 9.

2 collapsed lines
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: boolVariables=6060
%%%mzn-stat: failures=651221
%%%mzn-stat: propagations=323892330
3 collapsed lines
%%%mzn-stat: solveTime=61.1534
%%%mzn-stat: nSolutions=1
%%%mzn-stat-end
Finished in 1m 1s.

Size 10#

At size 10, there are 55 parts to pack on a board of 3025 squares, increasing the difficulty even more. Here OR-Tools CP-SAT is starting to struggle a bit more, and in two runs took about 13 and a half minutes. Here are the statistics for one of them.

2 collapsed lines
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: boolVariables=9319
%%%mzn-stat: failures=6516108
%%%mzn-stat: propagations=3208777732
3 collapsed lines
%%%mzn-stat: solveTime=804.504
%%%mzn-stat: nSolutions=1
%%%mzn-stat-end
Finished in 13m 25s.

As can be seen below, the two solutions found are quite different from each other.

Size 11#

Turning it up to eleven, it took OR-Tools CP-SAT a bit more than 51 minutes to solve the problem. With 66 parts and an area of 4356, it is significantly larger than size 10.

2 collapsed lines
%%%mzn-stat: objective=0
%%%mzn-stat: objectiveBound=0
%%%mzn-stat: boolVariables=13850
%%%mzn-stat: failures=15611863
%%%mzn-stat: propagations=10240280820
3 collapsed lines
%%%mzn-stat: solveTime=3078.61
%%%mzn-stat: nSolutions=1
%%%mzn-stat-end
Finished in 51m 19s.

Size 12#

Finding a solution of size 12 turned out to be too hard for the model. Running OR-Tools CP-SAT for 12 hours gave no result.

Other solvers#

In the above tests, only the OR-Tools CP-SAT solver has been used. This is both because initial experiments showed it was probably the best solver for this and because it has been dominant in the MiniZinc Challenge for more than a decade. A benefit of MiniZinc is that many different solvers can be tested, so let’s look at some alternatives.

Huub#

The new Huub solver was quite impressive in this year’s MiniZinc Challenge coming in third after OR-Tools CP-SAT and Chuffed in the Open category. Huub uses an external SAT solver, and runs single threaded. Running the model for size 8 with free search for ten rounds solves it in between 7.8 and 7.9 seconds, which is remarkably stable.

%%%mzn-stat: solveTime=7.474344917
%%%mzn-stat: failures=103390
%%%mzn-stat: peakDepth=4796
%%%mzn-stat: propagations=20878861
%%%mzn-stat: restarts=145
3 collapsed lines
%%%mzn-stat: oracleDecisions=123909
%%%mzn-stat: userDecisions=0
%%%mzn-stat-end
Finished in 7s 839msec.

This looked very promising, but increasing to size 9 Huub timed out after 12 hours.

Pumpkin 🎃#

Pumpkin is also an LCG solver like Huub, but it is more focused on proof logging. It is single-threaded like Huub, and uses a custom internal SAT solver. Here, solving size 8 took around 2 minutes (2 test runs).

%%%mzn-stat: nodes=838498
%%%mzn-stat: failures=427421
%%%mzn-stat: restarts=1706
%%%mzn-stat: variables=12219
%%%mzn-stat: propagators=14931
%%%mzn-stat: propagations=422962879
%%%mzn-stat: peakDepth=570
4 collapsed lines
%%%mzn-stat: nogoods=427421
%%%mzn-stat: backjumps=307815
%%%mzn-stat: solveTime=147.569079042
%%%mzn-stat-end
Finished in 2m 28s.

While size 8 was significantly slower for Pumpkin than for Huub, Pumpkin could actually solve size 9 in around 10 minutes.

%%%mzn-stat: nodes=3080585
%%%mzn-stat: failures=1547051
%%%mzn-stat: restarts=5243
%%%mzn-stat: variables=19208
%%%mzn-stat: propagators=23411
%%%mzn-stat: propagations=1673243823
%%%mzn-stat: peakDepth=925
4 collapsed lines
%%%mzn-stat: nogoods=1547051
%%%mzn-stat: backjumps=1108974
%%%mzn-stat: solveTime=642.870503792
%%%mzn-stat-end
Finished in 10m 43s.

Running size 10 with Pumpkin failed with an unspecified error after around 5 hours.

Chuffed, Gecode, and HiGHS#

None of these solvers were really useful for this problem. Chuffed is often a very good solver with really great automatic search heuristics, but sometimes it doesn’t work that well. Here, it took just over two hours to find a solution to the base size 8 packing. Chuffed is single-threaded, same as Huub and Pumpkin.

%%%mzn-stat: nodes=123635519
%%%mzn-stat: failures=60596049
%%%mzn-stat: restarts=73648
%%%mzn-stat: variables=34838
%%%mzn-stat: intVars=2734
%%%mzn-stat: boolVariables=32102
%%%mzn-stat: propagators=5521
%%%mzn-stat: propagations=96201866795
%%%mzn-stat: peakDepth=272
10 collapsed lines
%%%mzn-stat: nogoods=60596049
%%%mzn-stat: backjumps=59399516
%%%mzn-stat: peakMem=0.00
%%%mzn-stat: time=7432.788
%%%mzn-stat: initTime=0.078
%%%mzn-stat: solveTime=7432.710
%%%mzn-stat: baseMem=0.00
%%%mzn-stat: trailMem=0.12
%%%mzn-stat: randomSeed=-499155368
%%%mzn-stat-end
Finished in 2h 3m 53s.

Gecode is a competent classical constraint programming solver, and as such it doesn’t really have any effective automatic search heuristics. This is clearly visible for this problem, where it fails to solve the problem in 12 hours.

3 collapsed lines
%%%mzn-stat: initTime=0.0371863
%%%mzn-stat: solveTime=43199.8
%%%mzn-stat: solutions=0
%%%mzn-stat: variables=9550
%%%mzn-stat: propagators=9245
%%%mzn-stat: propagations=2213651468397
%%%mzn-stat: nodes=5060871988
%%%mzn-stat: failures=2530435922
%%%mzn-stat: restarts=0
%%%mzn-stat: peakDepth=108
%%%mzn-stat-end
Finished in 12h.

Since Gecode can really benefit from a search heuristic, I tried adding one. This heuristic uses the well-known left-bottom placement strategy, prioritizing placement of larger parts before placing smaller parts. This did not help.

% The position of a Part is essentially the index of the square.
array[Parts] of var int: position = [
    x[p] * card(Pos) + y[p]
  | p in Parts
];

% Search by placing the part with the smallest position/index at that position,
% breaking ties by input order (where larger parts are earlier).
solve :: int_search(position, smallest, indomain_min)
satisfy;

Finally, HiGHS is a modern open source MIP solver. Unfortunately, it also fails to solve this problem in 12 hours.

SICStus#

As mentioned above, the original development of the geost constraint was done in the SICStus Prolog solver. However, the MiniZinc model here does not translate to the geost constraint, nor is there support for using the specialized settings for the geost constraint.

Running the base MiniZinc model takes more than 4 hours to solve size 8.

2 collapsed lines
%%%mzn-stat: initTime=0.075
%%%mzn-stat: solveTime=15488.9
%%%mzn-stat: propagations=32188610389
3 collapsed lines
%%%mzn-stat: entailments=17947678933
%%%mzn-stat: prunings=58276738702
%%%mzn-stat: backtracks=381834851
%%%mzn-stat: restarts=0
%%%mzn-stat: solutions=1
%%%mzn-stat: optimalities=0
%%%mzn-stat: propagators=6651
%%%mzn-stat: variables=16508
%%%mzn-stat-end
Finished in 4h 24m 11s.

However, in the SICStus distribution, there is a partridge packing example with suitable geost arguments and a custom search predicate. Here we get the chance to compare a generic model with one that is customized for a solver, using that particular solver’s special features.

SICStus 4.10.1 (arm64-darwin-21.0.0): Sat Jun 28 12:23:49 CEST 2025
Licensed to Mikael Zayenz Lagerkvist
| ?- compile(user).
% compiling user...
| call_time(G,T) :-
   statistics(runtime,[T0|_]),
   G,
   statistics(runtime,[T1|_]),
   T is T1 - T0.
| ^D
% compiled user in module user, 2 msec 768 bytes
yes
| ?- ['lib/sicstus-4.10.1/library/clpfd/examples/partridge.pl'].
56 collapsed lines
% compiling /Users/zayenz/solvers/sicstus/lib/sicstus-4.10.1/library/clpfd/examples/partridge.pl...
%  module partridge imported into user
%  loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/lists.po...
%  module lists imported into partridge
%   loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/types.po...
%   module types imported into lists
%   loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/types.po in module types, 1 msec 6416 bytes
%  loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/lists.po in module lists, 3 msec 204320 bytes
%  loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/trees.po...
%  module trees imported into partridge
%  loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/trees.po in module trees, 1 msec 16336 bytes
%  loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/clpfd.po...
%  module clpfd imported into partridge
%   loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/atts.po...
%   module attributes imported into clpfd
%    module types imported into attributes
%   loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/atts.po in module attributes, 1 msec 32704 bytes
%   loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/fvar.po...
%    loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/ordsets.po...
%    module ordsets imported into fvar
%    loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/ordsets.po in module ordsets, 1 msec 50416 bytes
%    module attributes imported into fvar
%    module attributes imported into fvar
%   loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/fvar.po in module fvar, 1 msec 65376 bytes
%   loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/avl.po...
%   module avl imported into clpfd
%   loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/avl.po in module avl, 1 msec 68848 bytes
%   module lists imported into clpfd
%   module ordsets imported into clpfd
%   module trees imported into clpfd
%   module types imported into clpfd
%   loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/terms.po...
%   module terms imported into clpfd
%    module types imported into terms
%    module avl imported into terms
%   loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/terms.po in module terms, 1 msec 52656 bytes
%   loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/timeout.po...
%   module timeout imported into clpfd
%   loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/timeout.po in module timeout, 0 msec 1536 bytes
%   loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/ugraphs.po...
%   module ugraphs imported into clpfd
%    module ordsets imported into ugraphs
%    module lists imported into ugraphs
%    module avl imported into ugraphs
%    loading /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/random.po...
%    module random imported into ugraphs
%     module types imported into random
%     loading foreign resource /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/arm64-darwin-21.0.0/random.bundle in module random
%    loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/random.po in module random, 1 msec 31008 bytes
%    module types imported into ugraphs
%   loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/ugraphs.po in module ugraphs, 2 msec 104000 bytes
%   loading foreign resource /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/arm64-darwin-21.0.0/clpfd.bundle in module clpfd
%   module attributes imported into clpfd
%  loaded /Users/zayenz/solvers/sicstus/bin/sp-4.10.1/sicstus-4.10.1/library/clpfd.po in module clpfd, 19 msec 2004432 bytes
% compiled /Users/zayenz/solvers/sicstus/lib/sicstus-4.10.1/library/clpfd/examples/partridge.pl in module partridge, 29 msec 2250112 bytes
yes
| ?- call_time(partridge(8), T8).
placement space = 36x36
rectangles r(X,W,Y,H) = [r(1,8,1,8),r(1,8,13,8),r(1,8,21,8),r(1,8,29,8),r(9,8,1,8),r(9,8,9,8),r(9,8,17,8),r(17,8,1,8),r(9,7,30,7),r(16,7,30,7),r(17,7,18,7),r(23,7,30,7),r(30,7,16,7),r(30,7,23,7),r(30,7,30,7),r(24,6,18,6),r(24,6,24,6),r(25,6,1,6),r(25,6,7,6),r(31,6,1,6),r(31,6,7,6),r(9,5,25,5),r(14,5,25,5),r(17,5,13,5),r(19,5,25,5),r(22,5,13,5),r(1,4,9,4),r(5,4,9,4),r(17,4,9,4),r(21,4,9,4),r(27,3,15,3),r(31,3,13,3),r(34,3,13,3),r(27,2,13,2),r(29,2,13,2),r(30,1,15,1)]
T8 = 522 ?
yes
| ?- call_time(partridge(9), T9).
placement space = 45x45
rectangles r(X,W,Y,H) = [r(1,9,1,9),r(1,9,10,9),r(1,9,19,9),r(1,9,28,9),r(1,9,37,9),r(10,9,1,9),r(10,9,10,9),r(10,9,19,9),r(29,9,1,9),r(10,8,38,8),r(18,8,38,8),r(24,8,30,8),r(26,8,38,8),r(38,8,1,8),r(38,8,16,8),r(38,8,24,8),r(38,8,32,8),r(10,7,31,7),r(17,7,31,7),r(24,7,16,7),r(24,7,23,7),r(31,7,16,7),r(31,7,23,7),r(39,7,9,7),r(19,6,6,6),r(27,6,10,6),r(32,6,30,6),r(33,6,10,6),r(34,6,40,6),r(40,6,40,6),r(19,5,1,5),r(19,5,16,5),r(19,5,21,5),r(19,5,26,5),r(24,5,1,5),r(19,4,12,4),r(23,4,12,4),r(25,4,6,4),r(34,4,36,4),r(10,3,28,3),r(13,3,28,3),r(16,3,28,3),r(25,2,10,2),r(32,2,36,2),r(38,1,9,1)]
T9 = 60575 ?
yes
| ?- call_time(partridge(10), T10).
placement space = 55x55
rectangles r(X,W,Y,H) = [r(1,10,1,10),r(1,10,11,10),r(1,10,21,10),r(1,10,31,10),r(1,10,41,10),r(11,10,1,10),r(11,10,11,10),r(11,10,21,10),r(11,10,31,10),r(11,10,41,10),r(21,9,1,9),r(21,9,10,9),r(21,9,19,9),r(27,9,31,9),r(29,9,47,9),r(30,9,1,9),r(38,9,47,9),r(39,9,1,9),r(47,9,47,9),r(21,8,40,8),r(21,8,48,8),r(37,8,10,8),r(40,8,27,8),r(48,8,1,8),r(48,8,23,8),r(48,8,31,8),r(48,8,39,8),r(29,7,40,7),r(30,7,10,7),r(30,7,17,7),r(30,7,24,7),r(37,7,18,7),r(49,7,9,7),r(49,7,16,7),r(21,6,28,6),r(21,6,34,6),r(36,6,35,6),r(36,6,41,6),r(42,6,35,6),r(42,6,41,6),r(1,5,51,5),r(6,5,51,5),r(11,5,51,5),r(16,5,51,5),r(44,5,18,5),r(36,4,31,4),r(44,4,23,4),r(45,4,10,4),r(45,4,14,4),r(27,3,28,3),r(37,3,25,3),r(37,3,28,3),r(40,2,25,2),r(42,2,25,2),r(48,1,9,1)]
T10 = 1377485 ?
yes
| ?- call_time(partridge(11), T11).
placement space = 66x66
rectangles r(X,W,Y,H) = [r(1,11,1,11),r(1,11,12,11),r(1,11,23,11),r(1,11,34,11),r(1,11,45,11),r(1,11,56,11),r(12,11,1,11),r(12,11,12,11),r(12,11,23,11),r(12,11,45,11),r(12,11,56,11),r(23,10,1,10),r(23,10,11,10),r(23,10,49,10),r(33,10,1,10),r(33,10,11,10),r(39,10,25,10),r(43,10,9,10),r(57,10,32,10),r(57,10,42,10),r(57,10,52,10),r(23,9,21,9),r(39,9,40,9),r(39,9,49,9),r(39,9,58,9),r(48,9,40,9),r(48,9,49,9),r(48,9,58,9),r(49,9,23,9),r(58,9,23,9),r(12,8,34,8),r(23,8,59,8),r(25,8,41,8),r(31,8,59,8),r(43,8,1,8),r(49,8,32,8),r(51,8,1,8),r(59,8,1,8),r(20,7,34,7),r(32,7,21,7),r(32,7,28,7),r(53,7,9,7),r(53,7,16,7),r(60,7,9,7),r(60,7,16,7),r(27,6,35,6),r(33,6,35,6),r(33,6,41,6),r(33,6,47,6),r(33,6,53,6),r(43,6,19,6),r(27,5,30,5),r(39,5,35,5),r(44,5,35,5),r(57,5,62,5),r(62,5,62,5),r(21,4,41,4),r(23,4,30,4),r(39,4,21,4),r(49,4,19,4),r(12,3,42,3),r(15,3,42,3),r(18,3,42,3),r(23,2,45,2),r(23,2,47,2),r(20,1,41,1)]
T11 = 269799 ?
yes
| ?- call_time(partridge(12), T12).
placement space = 78x78
rectangles r(X,W,Y,H) = [r(1,12,1,12),r(1,12,13,12),r(1,12,25,12),r(1,12,37,12),r(1,12,49,12),r(1,12,61,12),r(13,12,1,12),r(13,12,19,12),r(13,12,31,12),r(13,12,43,12),r(13,12,55,12),r(13,12,67,12),r(25,11,1,11),r(25,11,24,11),r(25,11,35,11),r(25,11,46,11),r(25,11,57,11),r(25,11,68,11),r(36,11,1,11),r(44,11,26,11),r(44,11,47,11),r(47,11,1,11),r(58,11,1,11),r(34,10,12,10),r(45,10,37,10),r(59,10,50,10),r(59,10,60,10),r(69,10,1,10),r(69,10,11,10),r(69,10,30,10),r(69,10,40,10),r(69,10,50,10),r(69,10,60,10),r(25,9,12,9),r(36,9,38,9),r(51,9,12,9),r(52,9,70,9),r(60,9,12,9),r(61,9,21,9),r(61,9,70,9),r(70,9,21,9),r(70,9,70,9),r(36,8,22,8),r(36,8,30,8),r(36,8,47,8),r(36,8,55,8),r(36,8,63,8),r(36,8,71,8),r(44,8,63,8),r(44,8,71,8),r(44,7,12,7),r(44,7,19,7),r(52,7,63,7),r(55,7,36,7),r(55,7,43,7),r(62,7,36,7),r(62,7,43,7),r(1,6,73,6),r(7,6,73,6),r(13,6,13,6),r(19,6,13,6),r(55,6,26,6),r(63,6,30,6),r(44,5,58,5),r(49,5,58,5),r(51,5,21,5),r(54,5,58,5),r(56,5,21,5),r(55,4,32,4),r(55,4,50,4),r(55,4,54,4),r(59,4,32,4),r(25,3,21,3),r(28,3,21,3),r(31,3,21,3),r(34,2,22,2),r(61,2,30,2),r(44,1,37,1)]
T12 = 4276951 ?
yes
| ?-

Solving size 8 is really quick at around half a second (the timing is reported in milliseconds). Note also that SICStus is a single-threaded system. Size 9 took about a minute, size 10 took around 23 minutes, size 11 took 4 and a half minutes, and size 12 took 1 hour 11 minutes. It is expected that a larger instance can sometimes be faster (11 vs 10) when searching for a satisfying solution. Another things that is also worth noting that SICStus uses less than 20 MiB of memory when searching for a solution for size 12, while OR-Tools CP-SAT uses over 3 GiB of memory.

Here is the size 12 partridge packing that SICStus found. Since size 12 is the reason the Partridge Packing Problem got its name, it feels good to find a solution for this size as well.

At size 13, SICStus also starts to struggle with the search, with no solution produced in 12 hours.

Summary#

Solving the Partridge Packing Problem using MiniZinc is an interesting challenge. The base model performs poorly, and the usual trick (adding cumulative constraints) for improving a packing problem was not that useful. However, with some custom implied constraints and symmetry breaking, it was possible to get solutions for size 8 and 9 quite quickly.

As is common for CP problems modeled in MiniZinc, OR-Tools CP-SAT dominates in performance. However, it was interesting to see that the relatively new solvers Huub and Pumpkin are both promising.⁸ Moving from MiniZinc to the custom SICStus Partridge program showed the benefits of using a system with smart propagators and a custom search strategy.

There are better ways to solve this packing problem, giving faster solutions in a more scalable way. Still, it is a good example of how to incrementally develop a MiniZinc model and how to add strengthening constraints. A benefit of using a high-level modeling language for this type of problem is that it can be adapted to new constraints and changes in requirements. In many industrial problems, it is quite common for requirements to change frequently.

In the end though, the most important part was that it was fun to experiment with.