This is the sixtieth seventh part of the ILP series. For your convenience you can find other parts in the table of contents in Part 1 – Boolean algebra
Today we are going to solve Parcelacja.
We have a rectangular board without top left corner. We need to put two 3×2 rectangles and 5 triminoes to fill the board. Additional constraint is: no two triminoes can make a 3×2 rectangle. Rotations are allowed.
Solution:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 |
var board = new string[]{ " 000000", "0000000", "0000000", "0000000" }; var shapes = new Dictionary<string[][], int>(){ { new string[][] { new string[] { "000", "000" }, new string[] { "00", "00", "00" } }, 2}, { new string[][] { new string[] { "00", "0 " }, new string[] { "00", " 0" }, new string[] { " 0", "00" }, new string[] { "0 ", "00" } }, 5} }; var width = board[0].Length; var height = board.Length; var solver = new CplexMilpSolver(new CplexMilpSolverSettings { IntegerWidth = 12, Epsilon = 0.1, StoreDebugExpressions = false, CacheConstants = false, StoreDebugConstraints = false }); Console.WriteLine("Create variables"); var fields = Enumerable.Range(0, height).Select(r => Enumerable.Range(0, width).Select(c => board[r][c] != ' ' ? new System.Dynamic.ExpandoObject() : null).ToArray()).ToArray(); var allNonNullFields = fields.SelectMany(r => r).Where(f => f != null).ToArray(); Console.WriteLine("Set ranges"); for(int row = 0; row < height;++row){ for(int column = 0; column < width; ++ column){ dynamic self = fields[row][column]; if(self == null)continue; self.PieceId = solver.CreateAnonymous(Domain.PositiveOrZeroInteger); self.PieceType = solver.CreateAnonymous(Domain.PositiveOrZeroInteger); self.AllowedShapes = new List<int>(); solver.Set<LessOrEqual>(((IVariable)self.PieceType), solver.FromConstant(shapes.Keys.Sum(v => v.Length))); solver.Set<Equal>(solver.FromConstant(1), solver.Operation<MaterialImplication>(solver.Operation<IsNotEqual>((IVariable)self.PieceType, solver.FromConstant(0)), solver.Operation<IsEqual>((IVariable)self.PieceId, solver.FromConstant(row * width + column + 1)))); } } Console.WriteLine("Make shape rules"); for(int row = 0; row < height;++row){ for(int column = 0; column < width; ++ column){ string[][] sketches = shapes.Keys.SelectMany(v => v).ToArray(); for(int sketchId=0;sketchId < sketches.Length;++sketchId){ bool isAccessible = true; System.Dynamic.ExpandoObject designatedField = null; List<IVariable> parts = new List<IVariable>(); string[] sketch = sketches[sketchId]; for(int y=0;y<sketch.Length;++y){ for(int x=0;x<sketch[0].Length;++x){ if(sketch[y][x] != ' '){ if(row + y >= height || column + x >= width){ isAccessible = false; break; } dynamic self = fields[row+y][column+x]; if(self == null){ isAccessible = false; break; } if(designatedField == null){ designatedField = self; } parts.Add(self.PieceId); } } } if(isAccessible){ ((dynamic)designatedField).AllowedShapes.Add(sketchId); foreach(var field in parts){ solver.Set<Equal>( solver.FromConstant(1), solver.Operation<MaterialImplication>( solver.Operation<IsEqual>(solver.FromConstant(sketchId + 1), (IVariable)((dynamic)designatedField).PieceType), solver.Operation<IsEqual>(parts[0], field) ) ); } } } } } Console.WriteLine("Prohibit invalid shapes"); for(int row = 0; row < height;++row){ for(int column = 0; column < width; ++ column){ if(fields[row][column] == null)continue; dynamic self = fields[row][column]; List<int> allowedShapes = self.AllowedShapes; string[][] sketches = shapes.Keys.SelectMany(v => v).ToArray(); for(int sketchId=0;sketchId < sketches.Length;++sketchId){ if(!allowedShapes.Contains(sketchId)){ solver.Set<Equal>(solver.FromConstant(1), solver.Operation<IsNotEqual>((IVariable)self.PieceType, solver.FromConstant(sketchId + 1))); } } } } Console.WriteLine("Make shape counts"); int lowerBound = 1; foreach(var sketches in shapes.Keys){ int upperBound = lowerBound + sketches.Length - 1; int expectedCount = shapes[sketches]; solver.Set<Equal>( solver.Operation<Addition>( allNonNullFields.Select(f => ((IVariable)((dynamic)f).PieceType)).Select(f => solver.Operation<Conjunction>( solver.Operation<IsGreaterOrEqual>(f, solver.FromConstant(lowerBound)), solver.Operation<IsLessOrEqual>(f, solver.FromConstant(upperBound)) )).ToArray() ), solver.FromConstant(expectedCount) ); lowerBound = upperBound + 1; } Console.WriteLine("Ban trimino creating horizontal rectangle"); for(int row = 0; row < height;++row){ for(int column = 0; column < width; ++ column){ if(fields[row][column] == null)continue; if(column + 2 >= width)continue; if(fields[row][column+2] == null)continue; dynamic self = fields[row][column]; dynamic next = fields[row][column + 2]; solver.Set<Equal>(solver.FromConstant(0), solver.Operation<Conjunction>(solver.Operation<IsEqual>(solver.FromConstant(6), (IVariable)self.PieceType), solver.Operation<IsEqual>(solver.FromConstant(4), (IVariable)next.PieceType))); solver.Set<Equal>(solver.FromConstant(0), solver.Operation<Conjunction>(solver.Operation<IsEqual>(solver.FromConstant(3), (IVariable)self.PieceType), solver.Operation<IsEqual>(solver.FromConstant(5), (IVariable)next.PieceType))); } } Console.WriteLine("Ban trimino creating vertical rectangle"); for(int row = 0; row < height;++row){ for(int column = 0; column < width; ++ column){ if(fields[row][column] == null)continue; if(row + 1 >= height)continue; if(column + 1 >= width)continue; if(fields[row+1][column+1] == null)continue; dynamic self = fields[row][column]; dynamic next = fields[row + 1][column + 1]; solver.Set<Equal>(solver.FromConstant(0), solver.Operation<Conjunction>(solver.Operation<IsEqual>(solver.FromConstant(3), (IVariable)self.PieceType), solver.Operation<IsEqual>(solver.FromConstant(5), (IVariable)next.PieceType))); } } Console.WriteLine("Ban trimino creating vertical rectangle"); for(int row = 0; row < height;++row){ for(int column = 0; column < width; ++ column){ if(fields[row][column] == null)continue; if(row + 1 >= height)continue; if(fields[row+1][column] == null)continue; dynamic self = fields[row][column]; dynamic next = fields[row + 1][column]; solver.Set<Equal>(solver.FromConstant(0), solver.Operation<Conjunction>(solver.Operation<IsEqual>(solver.FromConstant(4), (IVariable)self.PieceType), solver.Operation<IsEqual>(solver.FromConstant(6), (IVariable)next.PieceType))); } } solver.AddGoal("goal", solver.FromConstant(0)); solver.Cplex.SetParam(ILOG.CPLEX.Cplex.Param.MIP.Pool.Intensity, 4); solver.Cplex.SetParam(ILOG.CPLEX.Cplex.Param.MIP.Limits.Populate, 1000); solver.Cplex.SetParam(ILOG.CPLEX.Cplex.Param.MIP.Pool.AbsGap, 0); solver.Cplex.Populate(); Console.WriteLine(solver.Cplex.SolnPoolNsolns); for(int i=0;i<solver.Cplex.SolnPoolNsolns;++i){ for(int row = 0; row < height;++row){ for(int column = 0; column < width; ++ column){ Console.Write(fields[row][column] != null ? Math.Round(solver.Cplex.GetValue(((ICplexVariable)(((dynamic)fields[row][column]).PieceId)).Var, i)).ToString("00") + "," : " ,"); } Console.WriteLine(); } Console.WriteLine(); for(int row = 0; row < height;++row){ for(int column = 0; column < width; ++ column){ Console.Write(fields[row][column] != null ? Math.Round(solver.Cplex.GetValue(((ICplexVariable)(((dynamic)fields[row][column]).PieceType)).Var, i)).ToString("00") + "," : " ,"); } Console.WriteLine(); } Console.WriteLine(); Console.WriteLine(); Console.WriteLine(); } |
Let’s analyze the code. First, in lines 1-6 we define the board. Spaces represent missing fields, zeroes show allowed positions.
Next, in lines 8-38 we define shapes with rotations. Each two dimensional array of strings represent one “shape” with all rotations. We also specify number of expected pieces of given kind.
Next, in lines 40-50 we prepare solver and variables.
In lines 52-70 we define our model. Each field is represented by expando. Each dynamic object has 3 fields: PieceId for numbering pieces (and making sure they are continuous etc). PieceType for deciding which shape and rotation we use (we number them from 1). AllowedShapes is to ban invalid positions to start a shape. In line 66 we make sure we only use valid number of piece type. In line 68 we make sure that if a field starts a shape then its PieceId is equal to the position on the board (so we uniquely identify each field).
Next, we define shapes. In lines 72-122 we do the magic. We iterate over each field. Next, for each field we iterate over each shape. We check every part of the shape and check if it’s non-empty (line 86). Next, we check if we didn’t jump out of the board (line 87), if there is a field over there (line 92). For each shape we define a “designatedField” being first non-empty part of the shape. So in line 98 we check if we already found a designated field and otherwise we set it. Next, in line 103 we store field for this shape’s part.
Next, in lines 107-119 we are sure that current shape can fit on the board and that we have a designated part. So we remember that designated part can hold this shape (line 108) and then in lines 110-117 we encode the constraint that all parts for the shape must have PieceId equal to the designated part if we decide to use the shape here.
Next, in lines 124-140 we prohibit invalid shapes (invalid as they don’t fit on the board in given position). We iterate over all fields, check if given field can be a designated point for given shape, and if not then we explicitly mark it as banned.
In lines 142-159 we count shapes and make sure we meet constraints (2 rectangles and 5 triminoes).
Finally, in lines 161-204 we ban triminoes from forming rectangles. In lines 161-175 we ban horizontal rectangle, then in line 177-190 we ban one type of vertical rectangle and then the another in lines 192-204.
Next, we add the goal in line 206. Next, in lines 208-209 we populate the solution pool (to find all solutions). Finally, we print them out. First part prints ids (so you can see how pieces were placed) and then we print the designated parts determining the type.
That’s it. The output:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 |
Populate: phase I Tried aggregator 2 times. MIP Presolve eliminated 3090 rows and 1572 columns. MIP Presolve modified 6123 coefficients. Aggregator did 2528 substitutions. Reduced MIP has 4862 rows, 2299 columns, and 12196 nonzeros. Reduced MIP has 2252 binaries, 47 generals, 0 SOSs, and 0 indicators. Presolve time = 0.03 sec. (15.50 ticks) Probing fixed 10 vars, tightened 0 bounds. Probing changed sense of 23 constraints. Probing time = 0.13 sec. (22.20 ticks) Tried aggregator 1 time. MIP Presolve eliminated 1956 rows and 965 columns. MIP Presolve modified 43 coefficients. Reduced MIP has 2906 rows, 1334 columns, and 8128 nonzeros. Reduced MIP has 1290 binaries, 44 generals, 0 SOSs, and 0 indicators. Presolve time = 0.02 sec. (8.51 ticks) Probing time = 0.00 sec. (3.78 ticks) Tried aggregator 1 time. Reduced MIP has 2906 rows, 1334 columns, and 8128 nonzeros. Reduced MIP has 1290 binaries, 44 generals, 0 SOSs, and 0 indicators. Presolve time = 0.02 sec. (3.79 ticks) Probing time = 0.01 sec. (3.80 ticks) Clique table members: 5449. 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. (16.53 ticks) Nodes Cuts/ Node Left Objective IInf Best Integer Best Bound ItCnt Gap 0 0 0.0000 463 0.0000 526 0 0 0.0000 298 Cuts: 528 845 0 0 0.0000 459 Cuts: 169 1045 0 0 0.0000 286 Cuts: 22 1125 0 0 0.0000 410 Cuts: 227 1353 0 2 0.0000 299 0.0000 1353 Elapsed time = 2.08 sec. (514.47 ticks, tree = 0.01 MB, solutions = 0) 310 93 0.0000 354 0.0000 19346 763 107 0.0000 414 0.0000 41850 1115 133 0.0000 416 0.0000 63794 1514 164 0.0000 449 0.0000 90196 1972 181 infeasible 0.0000 118810 2331 218 0.0000 560 0.0000 143347 2842 206 0.0000 459 0.0000 167442 3324 247 0.0000 265 0.0000 198386 3761 197 0.0000 488 0.0000 224293 5182 210 infeasible 0.0000 325163 Elapsed time = 17.41 sec. (3684.80 ticks, tree = 0.06 MB, solutions = 0) 7229 227 0.0000 453 0.0000 430764 7259 241 0.0000 432 0.0000 434041 7667 74 0.0000 435 0.0000 493078 7906 76 0.0000 452 0.0000 545085 8194 67 0.0000 455 0.0000 604941 8469 46 infeasible 0.0000 661968 * 8731+ 49 0.0000 0.0000 717433 0.00% 8753 42 0.0000 422 0.0000 0.0000 722464 0.00% GUB cover cuts applied: 73 Clique cuts applied: 54 Cover cuts applied: 222 Implied bound cuts applied: 60 Flow cuts applied: 157 Mixed integer rounding cuts applied: 71 Zero-half cuts applied: 23 Gomory fractional cuts applied: 5 Root node processing (before b&c): Real time = 2.11 sec. (512.94 ticks) Parallel b&c, 8 threads: Real time = 45.81 sec. (11210.76 ticks) Sync time (average) = 6.52 sec. Wait time (average) = 0.05 sec. ------------ Total (root+branch&cut) = 47.92 sec. (11723.70 ticks) Populate: phase II Nodes Cuts/ Node Left Objective IInf Best Integer Best Bound ItCnt Gap 8789 40 0.0000 542 0.0000 0.0000 730782 0.00% Elapsed time = 48.16 sec. (11754.16 ticks, tree = 0.06 MB, solutions = 1) 8852 47 0.0000 335 0.0000 0.0000 744239 0.00% 8946 45 0.0000 489 0.0000 0.0000 762836 0.00% 9026 44 infeasible 0.0000 0.0000 778848 0.00% 9119 22 0.0000 513 0.0000 0.0000 796941 0.00% 9241 25 0.0000 291 0.0000 0.0000 812743 0.00% 9334 10 infeasible 0.0000 0.0000 823224 0.00% 9362 2 0.0000 481 0.0000 0.0000 828305 0.00% GUB cover cuts applied: 73 Clique cuts applied: 54 Cover cuts applied: 222 Implied bound cuts applied: 60 Flow cuts applied: 157 Mixed integer rounding cuts applied: 71 Zero-half cuts applied: 23 Gomory fractional cuts applied: 5 Root node processing (before b&c): Real time = 0.06 sec. (0.34 ticks) Parallel b&c, 8 threads: Real time = 6.11 sec. (1829.20 ticks) Sync time (average) = 1.61 sec. Wait time (average) = 0.00 sec. ------------ Total (root+branch&cut) = 6.17 sec. (1829.54 ticks) 3 ,02,03,03,03,06,06, 02,02,03,03,03,13,06, 15,15,15,18,13,13,21, 15,15,15,18,18,21,21, ,05,01,00,00,04,00, 00,00,00,00,00,05,00, 01,00,00,06,00,00,05, 00,00,00,00,00,00,00, ,02,03,03,05,05,05, 02,02,10,03,05,05,05, 15,10,10,18,19,19,19, 15,15,18,18,19,19,19, ,05,04,00,01,00,00, 00,00,05,00,00,00,00, 06,00,00,05,01,00,00, 00,00,00,00,00,00,00, ,02,02,04,04,06,06, 08,08,02,04,12,12,06, 08,08,17,17,17,12,21, 08,08,17,17,17,21,21, ,04,00,03,00,04,00, 02,00,00,00,04,00,00, 00,00,01,00,00,00,05, 00,00,00,00,00,00,00, |
You can see it found 3 solutions, presented below: