3
$\begingroup$

I need to solve an optimization problem for a given set of polyominoes, for example the five Tetrominoes known from Tetris. The goal is to place each one of the figures on the 2d grid, so the area they enclose is maximal.

Staying with our Tetromino example, the optimal solution would be:

Maximal area enclosed with Tetrominoes

Bruteforcing the problem becomes unfeaseable very quickly, especially since my check with bfs search alone that an area is even enclosed at all must be performed countless times. I tried to implement systems to recognize when it won't be possible to close the circle with the remaining polyominoes, but didn't achieve a significant speed-up with them.

I also researched into strategies for the related polyomino tiling problem to see if I can adapt them, but had no success so far. The most promising approach I found was to translate it into Mixed-Integer Linear Programming like here or in this variation. But to apply this approach to my problem I need to both add the restriction of the enclosed area and then formulate its size as a linear objective function, whereby especially the encoding of the latter may not be possible at all.

Are there other algorithms to find the global maximum for this type of optimization problem that I've overlooked? Or are there data structures or strategies for branch elimination to drastically speed up the search of all possibilities?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

Here is one possible approach. Introduce boolean variables $x_\ell,y_\ell,z_{\ell,p,r}$, where $\ell$ ranges over locations in the grid, $p$ ranges over polyominoes, and $r$ over rotation amounts (0, 90, 180, 270), with the following meaning:

  • $x_\ell$ is true means that $\ell$ is part of the interior region that is enclosed (not including the boundary/perimeter)

  • $y_\ell$ is true means that $\ell$ is part of the boundary/perimeter (it is covered by a polyomino)

  • $z_{\ell,p,r}$ means that polyomino $p$ is rotated by $r$ degrees and then placed at location $\ell$

Now you can write down some boolean constraints regarding the relationships between these. For instance:

  • $x_\ell \land \neg y_{\ell'} \implies x_{\ell'}$, for all pairs of adjacent locations $\ell,\ell'$

  • $z_{\ell,p,r} \implies y_{\ell+t}$ for all offsets $t$ that are part of $p$ after being rotated by $r$ degrees

  • $y_{\ell} \implies \bigvee_{p,r,t} z_{\ell-t,p,r}$ where the logical-OR ranges over all $p,r$ and all offsets $t$ that are part of $p$ after being rotated by $r$ degrees

  • no two polyominoes overlap, which can be expressed in terms of $\neg (z_{\ell,p,r} \land z_{\ell',p',r'})$ where $p,p',r,r'$ range over all possibilities and then $\ell,\ell'$ have an appropriate relationship given $p,p',r,r'$

  • For each $p$, exactly one of $z_{\ell,p,r}$ is true (ranging over all $\ell,r$) - this is a one-out-of-$n$ constraint, which can be encoded in standard ways

  • $\neg x_\ell$, for each location $\ell$ along the four outermost edges of the grid of cells - this forces them to be on on the outside of the shape

  • $\neg y_\ell$, for each location $\ell$ along the four outermost edges of the grid of cells - this forces them to not be covered by any polyomino

Finally the goal is to find a satisfying assignment that maximizes the number of $x_\ell$ that are true. This could be solved using ILP (each variable can be treated as a 0-or-1 integer variable, each boolean constraint can be expressed as one or more linear inequalities, and the objective function is $\sum_\ell x_\ell$) or using SAT (use binary search on the number of cells the region encloses, add a pseudo-boolean constraint that requires at least that many of the $x_\ell$ are true). I don't know whether this will be efficient enough to work well in practice, but you could give it a try. That would let you benefit from all of the heuristics and algorithms that have been developed over the years for ILP and/or SAT.

Improve this answer
$\endgroup$
19
  • $\begingroup$ This is not enough. Your constraints do not stop the solver from setting $x_l$ to true everywhere (apart from the origin). Same goes for $y_l$ - the solver can set values that the boundary is everywhere. You can't use material implication just like that. You need to use equivalence (or sums). Or you need to add constraint that a field is either internal or external. $\endgroup$ Commented Aug 27, 2022 at 7:12
  • $\begingroup$ @user1543037, That doesn't seem right to me. The last constraint ensures that $x_0$ is false (at the origin). Then, the first constraint (about adjacent locations) implies that $x_\ell$ is false at the locations that are immediately adjacent to the origin, and so on, until you hit the boundary -- so the combination of the last constraint and the first constraint imply that $x_\ell$ is false everywhere outside the enclosed region. Am I missing something? $\endgroup$ Commented Aug 27, 2022 at 7:15
  • $\begingroup$ When $x_0$ is zero, then the left side of the implication is zero, so right side can be anything. $\endgroup$ Commented Aug 27, 2022 at 7:17
  • $\begingroup$ Same issue applies to the boundary. You add a constraint which makes sure that there is boundary formed on polyominoes. But there is no constraint making sure that there is no boundary outside of polyominoes. $\endgroup$ Commented Aug 27, 2022 at 7:18
  • $\begingroup$ @user1543037, I don't think that's right. Could you take another look? If $\ell$ is adjacent to the origin and the origin is not part of the boundary, consider the case $\ell'=0$; then the constraint $x_\ell \land \neg y_{\ell'} \implies x_{\ell'}$ (i.e., $x_\ell \land \neg y_0 \implies x_0$) combined with $x_0=\text{False}$ and $y_0=\text{False}$ implies $x_\ell=\text{False}$. This continues onward, implying that every cell in the exterior must have $x_\ell$ be false. Do you see something I'm missing? $\endgroup$ Commented Aug 27, 2022 at 7:23
2
$\begingroup$

In case anyone is interested, this is the solution implemented with https://github.com/afish/MilpManager

It's inspired by https://cs.stackexchange.com/a/153823/79175 (including the discussions we had in the comments). Source code:

var integerWidth = 15; var epsilon = 0.001; var storeDebugExpressions = false; var cacheConstants = false; var storeDebugConstraints = false; var solver = new CplexMilpSolver(new CplexMilpSolverSettings { IntegerWidth = integerWidth, Epsilon = epsilon, StoreDebugExpressions = storeDebugExpressions, CacheConstants = cacheConstants, StoreDebugConstraints = storeDebugConstraints }); solver.Cplex.SetParam(ILOG.CPLEX.Cplex.Param.MIP.Tolerances.Integrality, 0.000000000001); var pieces = new []{ new []{ "XX", "X ", "X " }, new[]{ " XX", "XX " }, new[]{ "XX ", " XX" }, new[]{ "X ", "XX", "X " }, new[]{ "XX", "XX" }, new[]{ "XXXX" }, new[]{ "X ", "X ", "XX" } }; var totalLength = pieces.Select(p => (int)Math.Max(p.Length, p[0].Length)).Sum(); var maxLength = pieces.Max(p => (int)Math.Max(p.Length, p[0].Length)); var width = 11; // Or totalLength + 2 to make sure it's big enough for any shapes var height = 9; // Or totalLength + 2 to make sure it's big enough for any shapes var piecesCount = pieces.Length; // Right, down, left, up var directions = 4; var heads = Enumerable.Range(0, height).Select(h => Enumerable.Range(0, width).Select(w => Enumerable.Range(0, piecesCount).Select(p => Enumerable.Range(0, directions).Select(d => solver.CreateAnonymous(Domain.BinaryInteger) ).ToArray() ).ToArray() ).ToArray() ).ToArray(); var isInterior = Enumerable.Range(0, height).Select(h => Enumerable.Range(0, width).Select(w => solver.CreateAnonymous(Domain.BinaryInteger) ).ToArray() ).ToArray(); var isBoundary = Enumerable.Range(0, height).Select(h => Enumerable.Range(0, width).Select(w => solver.CreateAnonymous(Domain.BinaryInteger) ).ToArray() ).ToArray(); Func<string[], int, string[]> rotatePiece = (piece, direction) => { if(direction == 0){ return piece; }else if(direction == 2){ return piece.Select(l => string.Join("", l.Reverse())).Reverse().ToArray(); }else if(direction == 1){ return Enumerable.Range(0, piece[0].Length) .Select(c => string.Join("", Enumerable.Range(0, piece.Length).Select(w => piece[piece.Length - 1 - w][c]))) .ToArray(); }else { return Enumerable.Range(0, piece[0].Length) .Select(c => string.Join("", Enumerable.Range(0, piece.Length).Select(w => piece[w][piece[0].Length - 1 - c]))) .ToArray(); } }; // Put pieces on the board for(int p=0;p<piecesCount;++p){ solver.Set<Equal>( solver.FromConstant(1), solver.Operation<Addition>( Enumerable.Range(0, height).SelectMany(h => Enumerable.Range(0, width).SelectMany(w => Enumerable.Range(0, directions).Select(d => heads[h][w][p][d] ) ) ).ToArray() ) ); } // Make sure pieces don't fall off the board for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ for(int p=0;p<piecesCount;++p){ for(int d = 0; d < directions;++d){ var piece = rotatePiece(pieces[p], d); if(d == 0){ if(w + piece[0].Length - 1 >= width){ heads[h][w][p][d].Set<Equal>(solver.FromConstant(0)); } }else if(d == 1){ if(h + piece.Length - 1 >= height){ heads[h][w][p][d].Set<Equal>(solver.FromConstant(0)); } }else if(d == 2){ if(w - piece[0].Length + 1 < 0){ heads[h][w][p][d].Set<Equal>(solver.FromConstant(0)); } }else if(d == 3){ if(h - piece.Length + 1 < 0){ heads[h][w][p][d].Set<Equal>(solver.FromConstant(0)); } } } } } } // Set impacts IList<Tuple<int, int, int, int>>[][] allImpacts = Enumerable.Range(0, height).Select(h => Enumerable.Range(0, width).Select(w => new List<Tuple<int, int, int, int>>() ).ToArray() ).ToArray(); for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ for(int p=0;p<piecesCount;++p){ for(int d = 0; d < directions;++d){ var piece = rotatePiece(pieces[p], d); var pieceWidth = piece[0].Length; var pieceHeight = piece.Length; for(int y=0;y<pieceHeight;++y){ for(int x=0;x<pieceWidth;++x){ var isPieceX = piece[y][x] == 'X'; if(!isPieceX) continue; var finalY = d == 0 || d == 1 ? h + y : y - pieceHeight + 1 + y; var finalX = d == 0 || d == 1 ? w + x : x - pieceWidth + 1 + x; if(finalY >= 0 && finalY < height && finalX >=0 && finalX < width){ allImpacts[finalY][finalX].Add(Tuple.Create(h, w, p, d)); } } } } } } } // Set the boundary for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ isBoundary[h][w].Set<Equal>( solver.Operation<Addition>( allImpacts[h][w].Select(t => heads[t.Item1][t.Item2][t.Item3][t.Item4]).ToArray() ) ); } } // Make sure pieces do not overlap for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ solver.Set<LessOrEqual>( solver.Operation<Addition>( allImpacts[h][w].Select(t => heads[t.Item1][t.Item2][t.Item3][t.Item4]).ToArray() ), solver.FromConstant(1) ); } } // Extend interior for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ for(int y=-1;y<=1;++y){ for(int x=-1;x<=1;++x){ if(y == 0 && x == 0){ continue; } if(h + y < 0 || h + y >= height){ continue; } if(w + x < 0 || w + x >= width){ continue; } solver.Set<Equal>( solver.FromConstant(1), solver.Operation<MaterialImplication>( solver.Operation<Conjunction>( isInterior[h][w], isBoundary[h + y][w + x].Operation<BinaryNegation>() ), isInterior[h + y][w + x] ) ); } } } } // Make sure there is some exterior for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ if(h == 0 || h == height - 1 || w == 0 || w == width - 1){ solver.Set<Equal>( solver.FromConstant(0), isInterior[h][w] ); solver.Set<Equal>( solver.FromConstant(0), isBoundary[h][w] ); } } } // Make sure a field cannot be both exterior and boundary for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ solver.Set<LessOrEqual>( solver.Operation<Addition>( isInterior[h][w], isBoundary[h][w] ), solver.FromConstant(1) ); } } var goal = solver.Operation<Addition>(isInterior.SelectMany(x => x).ToArray()); solver.AddGoal("goal", goal); solver.Solve(); Console.WriteLine(solver.GetValue(goal)); var board = Enumerable.Range(0, height).Select(h => new String(' ', width).ToCharArray()).ToArray(); for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ for(int p=0;p<piecesCount;++p){ for(int d = 0; d < directions;++d){ var piece = rotatePiece(pieces[p], d); var pieceWidth = piece[0].Length; var pieceHeight = piece.Length; if(solver.GetValue(heads[h][w][p][d]) > 0.5){ for(int y=0;y<pieceHeight;++y){ for(int x=0;x<pieceWidth;++x){ var finalY = d == 0 || d == 1 ? h + y : y - pieceHeight + 1 + y; var finalX = d == 0 || d == 1 ? w + x : x - pieceWidth + 1 + x; if(finalY >= 0 && finalY < height && finalX >=0 && finalX < width){ board[finalY][finalX] = piece[y][x]; }else { Console.WriteLine("Error in " + h + " " + w + " " + p + " " + d + " " + y + " " + x); } } } } } } } } foreach(var l in board){ Console.WriteLine(string.Join("", l)); } Console.WriteLine(); for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ Console.Write(Math.Round(solver.GetValue(isInterior[h][w]))); } Console.WriteLine(); } Console.WriteLine(); for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ Console.Write(Math.Round(solver.GetValue(isBoundary[h][w]))); } Console.WriteLine(); } Console.WriteLine();

It should generalize to any polyominoes you have (just replace the array at the beginning). This can be easily translated to SAT if you need, so you have plenty of solvers out there to play with it.

CPLEX output for the sample above:

Tried aggregator 2 times. MIP Presolve eliminated 4659 rows and 4747 columns. MIP Presolve modified 688 coefficients. Aggregator did 1432 substitutions. Reduced MIP has 809 rows, 847 columns, and 4647 nonzeros. Reduced MIP has 847 binaries, 0 generals, 0 SOSs, and 0 indicators. Presolve time = 0.03 sec. (10.92 ticks) Found incumbent of value 0.000000 after 0.06 sec. (33.44 ticks) Probing time = 0.00 sec. (2.70 ticks) Tried aggregator 1 time. Reduced MIP has 809 rows, 847 columns, and 4647 nonzeros. Reduced MIP has 847 binaries, 0 generals, 0 SOSs, and 0 indicators. Presolve time = 0.02 sec. (2.75 ticks) Probing time = 0.02 sec. (2.69 ticks) Clique table members: 1918. MIP emphasis: balance optimality and feasibility. MIP search method: dynamic search. Parallel mode: deterministic, using up to 8 threads. Root relaxation solution time = 0.02 sec. (17.35 ticks) Nodes Cuts/ Node Left Objective IInf Best Integer Best Bound ItCnt Gap * 0+ 0 0.0000 35.0000 565 --- 0 0 29.7678 243 0.0000 29.7678 565 --- * 0+ 0 2.0000 29.7678 657 --- 0 0 29.2337 242 2.0000 Cuts: 36 657 --- * 0+ 0 3.0000 29.2337 657 874.46% 0 0 28.7133 242 3.0000 Cuts: 99 816 857.11% 0 0 28.5944 248 3.0000 Cuts: 26 881 853.15% 0 0 28.3575 243 3.0000 Cuts: 21 958 845.25% 0 0 28.1942 256 3.0000 Cuts: 27 1036 839.81% 0 0 28.1076 244 3.0000 Cuts: 25 1120 836.92% * 0+ 0 4.0000 28.1076 1120 602.69% 0 0 28.1076 245 4.0000 Cuts: 12 1130 602.69% 0 0 28.0851 255 4.0000 Cuts: 9 1150 597.76% 0 0 27.9794 239 4.0000 Cuts: 28 1229 597.76% 0 0 27.8878 280 4.0000 Cuts: 25 1330 597.20% 0 0 27.8093 275 4.0000 Cuts: 23 1402 595.23% 0 0 27.6921 275 4.0000 Cuts: 40 1520 592.30% 0 0 27.6921 276 4.0000 Cuts: 7 1525 592.30% * 0+ 0 20.0000 27.6921 1525 38.46% 0 2 27.6921 276 20.0000 27.6060 1525 38.03% Elapsed time = 0.97 sec. (335.52 ticks, tree = 0.01 MB, solutions = 5) * 5+ 5 21.0000 27.5336 1944 31.11% * 56 40 integral 0 22.0000 27.4817 5646 24.92% * 63+ 38 25.0000 27.4817 5959 9.93% Clique cuts applied: 15 Cover cuts applied: 19 Implied bound cuts applied: 63 Mixed integer rounding cuts applied: 35 Zero-half cuts applied: 53 Gomory fractional cuts applied: 1 Root node processing (before b&c): Real time = 0.97 sec. (335.26 ticks) Parallel b&c, 8 threads: Real time = 0.52 sec. (250.13 ticks) Sync time (average) = 0.26 sec. Wait time (average) = 0.00 sec. ------------ Total (root+branch&cut) = 1.49 sec. (585.39 ticks) 25 XXXXXX XX X XX X X XX X X XX XX XXXXXXX 00000000000 00000000000 00001111000 00011111000 00111111000 00111111000 00011110000 00000000000 00000000000 00000000000 00011111100 00110000100 01100000100 01000000110 01000000100 01100001100 00111111100 00000000000
Improve this answer
$\endgroup$
3
  • $\begingroup$ This is awesome, thank you very much for the implementation :) $\endgroup$ Commented Aug 27, 2022 at 15:38
  • 1
    $\begingroup$ Can you share the main ideas? We're not a coding site, and we prefer answers with ideas, explanations, pseudocode, etc., over a big block of code. Thank you! $\endgroup$ Commented Aug 27, 2022 at 23:09
  • $\begingroup$ @D.W. As mentioned, it's the implementation of your post cs.stackexchange.com/a/153823/79175 There are no new "main ideas", it's nearly a transliteration. $\endgroup$ Commented Aug 28, 2022 at 6:35

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.