https://zayenz.se/blog/post/partridge-packing/ Solving the Partridge Packing Problem using MiniZinc Toggle menu Menu Navigation Home Research Blog On This Page * Solving the Partridge Packing Problem using MiniZinc * The Problem * The Base MiniZinc Model + Representing the Problem + The Viewpoint + Basic Packing Constraints + Search and output + Performance of base model * Improved model + Cumulative profile + Exact fill + Restrictions of placements close to the edge + Symmetry breaking for same size parts + How about the board symmetry? * The Full Program * Solving larger instances + Size 9 + Size 10 + Size 11 + Size 12 * Other solvers + Huub + Pumpkin + Chuffed, Gecode, and HiGHS + SICStus * Summary * Footnotes - Back to all blog posts - Newer Older - Solving the Partridge Packing Problem using MiniZinc DRAFT 2025-11-26 * 23 min read * Cite * Constraint programming MiniZinc models puzzles Part of the minizinc collection. A collection of posts about constraint programming using the MiniZinc modeling language, including tutorials, case studies, and benchmarking results. See all 4 posts in the collection 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.^1 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 nnn, the goal is to pack: * 1 square of size 1x11 \times 11x1 * 2 squares of size 2x22 \times 22x2 * ... * nnn squares of size nxnn \times nnxn into a square of size n(n+1)2xn(n+1)2\frac{n(n+1)}{2} \times \frac{n (n+1)}{2}2n(n+1) x2n(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 [?]i=1ni[?]i2=n(n+1)2[?]n(n+1)2\sum_{i=1}^{n} i \cdot i^2 = \frac{n(n+1)} {2} \cdot \frac{n(n+1)}{2}[?]i=1n i[?]i2=2n(n+1) [?]2n(n+1) 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 nnn 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.^3 122333444455555666666777777788888888 These parts can be packed in a square of size 36x3636\times 3636x36, where 363636 comes from 8x92=36\frac{8 \times 9}{2} = 3628x9 =36, and here is one such solution. 888888887777777666666555554444333221 This visualization shows how all the squares pack together perfectly to fill the 36x36 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 1x11\times 11x1 and 2 parts of size 2x22\times 22x2 that should be packed in a 3x33\times 33x3 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. 122 122 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 nnn, 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.^4 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. MiniZinc IDE instance dialog 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.^5 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. Terminal window 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. Visualization of packing 888888887777777666666555554444333221 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.^6 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. 02457902467911Capacity 812345466678TimeResource usage 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! Terminal window %%%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. Terminal window %%%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. Terminal window %%%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. 6area 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[?]6=122\cdot 6=122[?]6=12 needs to be available. However, for parts up to size 2, this is not possible since there is one 1x11\ times 11x1 square and two 2x22\times 22x2 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. 6221 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 nnn, the sum of the areas of the smaller parts (up to size n-1n-1n-1) is the square of the triangular number for n-1n-1n-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..12h - Huub 8s >12h - - - - Pumpkin 2m 28s 10m 43s >5h - - - Chuffed 2h 4m - - - - - Gecode >12h - - - - - HiGHS >12h - - - - - SICStus 4h 24m - - - - - SICStus Partridge 1s 1m 1s 23m 4m 30s 1h 11m >12h 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. Footnotes# 1. Personally, I think the name nooverlap is better, and that is the name used in Gecode. - 2. The formula n(n+1)2\frac{n(n+1)}{2}2n(n+1) is called the triangular number and represents the sum from 111 to nnn. - 3. For size n=8n=8n=8, we need to pack: + 1 square of size 1x11 \times 11x1 (area: 13=11^3=113=1) + 2 squares of size 2x22 \times 22x2 (area: 23=82^3=823=8) + 3 squares of size 3x33 \times 33x3 (area: 33=273^3=2733=27) + 4 squares of size 4x44 \times 44x4 (area: 43=644^3=6443=64) + 5 squares of size 5x55 \times 55x5 (area: 53=1255^3=12553=125 ) + 6 squares of size 6x66 \times 66x6 (area: 63=2166^3=21663=216 ) + 7 squares of size 7x77 \times 77x7 (area: 73=3437^3=34373=343 ) + 8 squares of size 8x88 \times 88x8 (area: 83=5128^3=51283=512 ) The total area is then 1+8+27+64+125+216+343+512=12961 + 8 + 27 + 64 + 125 + 216 + 343 + 512 = 12961+8+27+64+125+216+343+512=1296. The larger square has side lengths 8x92=36\frac{8 \times 9}{2} = 3628x9 =36, so the area is 36x36=129636 \times 36 = 129636x36= 1296, which matches. In Matt Parker's video, the size 9 is used since that gives a total area of 2025. - 4. Peter Stuckey has a good talk called There are no integers in discrete optimisation , where he argues that it is very important to use strong typing for domains in combinatorial optimization models. I generally agree, and try to use it when possible in MiniZinc models. Unfortunately, there are still some sharp corners in the experience, but it is improving. - 5. Given that it is so common to express features on aspects of a problem, I think it would be great if MiniZinc supported accessing features of records inside arrays. MiniZinc issue #970 suggests this as a feature. - 6. Note that the classical viewpoint for 8-queens in essence names the queens, but it also restricts it so that the first queen is always on the first row, the second queen on the second row, etc. This breaks the introduced symmetry directly. - 7. The construction [ [x[p], y[p]][x_or_y] | x_or_y in 1..2, p in PartsWithSize] is a bit unfortunate in my opinion. Having the variable x_or_y to index a temporary array is an artifact of MiniZinc only supporting one element in each iteration of the comprehension. In this issue there is a proposal for allowing multiple elements in each iteration, in which case the much more natural [ x[p], y[p] | p in PartsWithSize] could have been used instead. - 8. I hope that both Huub and Pumpkin will add parallelism in the future, as I think that would make the systems even more interesting to use. The simplest to add is portfolio parallelism, but the recent paper at CPAIOR 2025 from the OR-Tools team and Peter Stuckey on how to make parallel search work well with LCG solving is also very interesting. - Navigation * Home * Blog * Research * About On This Page * Solving the Partridge Packing Problem using MiniZinc DRAFT * The Problem DRAFT * The Base MiniZinc Model DRAFT + Representing the Problem DRAFT + The Viewpoint DRAFT + Basic Packing Constraints DRAFT + Search and output DRAFT + Performance of base model DRAFT * Improved model DRAFT + Cumulative profile DRAFT + Exact fill DRAFT + Restrictions of placements close to the edge DRAFT + Symmetry breaking for same size parts DRAFT + How about the board symmetry? DRAFT * The Full Program DRAFT * Solving larger instances DRAFT + Size 9 DRAFT + Size 10 DRAFT + Size 11 DRAFT + Size 12 DRAFT * Other solvers DRAFT + Huub DRAFT + Pumpkin DRAFT + Chuffed, Gecode, and HiGHS DRAFT + SICStus DRAFT * Summary DRAFT * Footnotes DRAFT Contact * fromweb@zayenz.se * GitHub * LinkedIn