This is the eightieth 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 Islanders game.

We have different types of buildings. Each building has a size and a range. When building is constructed, it automatically adds points for its base gain and for all buildings in its range constructed earlier (which may be negative). We want to maximize the income.

Let’s start with model:

1 2 3 4 5 6 |
class Building{ public string Name {get; set;} public int Range {get;set;} public int Size {get; set;} public int BaseScore {get;set;} } |

Each building has a name (for debugging only), range and size, and a base score (points added when the building is constructed).

Next, for each building we’ll create an ILP model:

1 2 3 4 5 6 7 8 |
class Placement{ public Building Building {get;set;} // Center of building public IVariable X {get;set;} public IVariable Y {get;set;} public IVariable BuildOrder {get;set;} public int Identifier {get;set;} } |

This has two variables for position, and third variable for build order. It also has an identifier for debugging.

Now, the code:

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 |
// Model var cityCenter = new Building(){ Name = "City Center", Range = 5, Size = 1, BaseScore = 15 }; var house = new Building(){ Name = "House", Range = 6, Size = 2, BaseScore = 1 }; var mansion = new Building(){ Name = "Mansion", Range = 3, Size = 2, BaseScore = 1 }; var fountain = new Building(){ Name = "Fountain", Range = 8, Size = 2, BaseScore = 0 }; var costs = new Dictionary<Tuple<Building, Building>, int>(){ { Tuple.Create(cityCenter, cityCenter), -40}, { Tuple.Create(cityCenter, house), 0}, { Tuple.Create(cityCenter, mansion), 0}, { Tuple.Create(cityCenter, fountain), 0}, { Tuple.Create(house, cityCenter), 6}, { Tuple.Create(house, house), 1}, { Tuple.Create(house, mansion), 0}, { Tuple.Create(house, fountain), 2}, { Tuple.Create(mansion, cityCenter), 8}, { Tuple.Create(mansion, house), 0}, { Tuple.Create(mansion, mansion), 1}, { Tuple.Create(mansion, fountain), 3}, { Tuple.Create(fountain, cityCenter), 7}, { Tuple.Create(fountain, house), 2}, { Tuple.Create(fountain, mansion), 3}, { Tuple.Create(fountain, fountain), -15} }; Building[] buildings = new Dictionary<Building, int>(){ { cityCenter, 1}, { house, 4}, { mansion, 0}, { fountain, 0} }.SelectMany(kv => Enumerable.Range(0, kv.Value).Select(i => kv.Key)).ToArray(); var boardSize = 20; // ILP Placement[] placements = buildings.Select((b, i) => new Placement(){ Building = b, X = solver.CreateAnonymous(Domain.PositiveOrZeroInteger), Y = solver.CreateAnonymous(Domain.PositiveOrZeroInteger), BuildOrder = solver.CreateAnonymous(Domain.PositiveOrZeroInteger), Identifier = i + 1 }).ToArray(); // Make sure buildings are on board foreach(var p in placements){ p.X.Set<LessOrEqual>(solver.FromConstant(boardSize - p.Building.Size)); p.Y.Set<LessOrEqual>(solver.FromConstant(boardSize - p.Building.Size)); } // Make sure buildings do not overlap foreach(var p in placements){ foreach(var p2 in placements){ if(p == p2)continue; var overlapX = p.X.Operation<Subtraction>(p2.X).Operation<AbsoluteValue>().Operation<Multiplication>(solver.FromConstant(2)).Operation<IsLessThan>(solver.FromConstant(p.Building.Size + p2.Building.Size)); var overlapY = p.Y.Operation<Subtraction>(p2.Y).Operation<AbsoluteValue>().Operation<Multiplication>(solver.FromConstant(2)).Operation<IsLessThan>(solver.FromConstant(p.Building.Size + p2.Building.Size)); var overlap = overlapX.Operation<Disjunction>(overlapY); overlap.Set<Equal>(solver.FromConstant(0)); } } // Make sure building order is different solver.Set<AllDifferent>(solver.FromConstant(0), placements.Select(p => p.BuildOrder).ToArray()); // Calculate cost IVariable baseInput = solver.Operation<Addition>(placements.Select(p => p.Building.BaseScore).Select(s => solver.FromConstant(s)).ToArray()); IVariable[] rangeInput = placements.SelectMany(p => { return placements.Where(p2 => p != p2).Select(p2 => { var isAfter = p.BuildOrder.Operation<IsGreaterOrEqual>(p2.BuildOrder); var isInRange = p.X.Operation<Subtraction>(p2.X).Operation<AbsoluteValue>().Operation<Addition>(p.Y.Operation<Subtraction>(p2.Y).Operation<AbsoluteValue>()).Operation<IsLessOrEqual>(solver.FromConstant(p.Building.Range)); return solver.Operation<Condition>( isAfter.Operation<Conjunction>(isInRange), solver.FromConstant(costs[Tuple.Create(p.Building, p2.Building)]), solver.FromConstant(0) ); }); }).ToArray(); var goal = baseInput.Operation<Addition>(rangeInput); solver.AddGoal("Goal", goal); solver.Solve(); Console.WriteLine("Result: " + solver.GetValue(goal)); var board = Enumerable.Range(0, boardSize).Select(i => new int[boardSize]).ToArray(); for(int i=0;i<boardSize;++i){ for(int j=0;j<boardSize;++j){ board[i][j] = 0; } } foreach(var p in placements){ int x = (int)Math.Round(solver.GetValue(p.X)); int y = (int)Math.Round(solver.GetValue(p.Y)); for(int i=0;i<p.Building.Size;++i){ for(int j=0;j<p.Building.Size;++j){ board[x+i][y+j] = p.Identifier; } } } for(int i=0;i<boardSize;++i){ for(int j=0;j<boardSize;++j){ if(board[i][j] != 0){ Console.Write(board[i][j]); }else{ Console.Write(" "); } } Console.WriteLine(); } |

In lines 1-25 we define buildings. Next, in lines 27-44 we specify points for constructing buildings in ranges. Finally, we ask to build one city center and for houses, just to keep model simple. We also specify that our board side length is 20.

We initialize variables in lines 55-62. Then, we start adding constraints.

First, we make sure that buildings are built on the board (lines 64-68).

Next, we want to make sure that buildings do not overlap. Building size is set to one means that it spans half a field each direction (so it forms a square of size one). To avoid floating point errors we multiply differences by two and not divide. This is in lines 70-83.

Next, we make sure that buildings are constructed in some order (line 86).

Finally, we calculate the cost in lines 88-101

Sample 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 |
Tried aggregator 2 times. MIP Presolve eliminated 545 rows and 105 columns. MIP Presolve modified 1807 coefficients. Aggregator did 456 substitutions. Reduced MIP has 1050 rows, 655 columns, and 3540 nonzeros. Reduced MIP has 560 binaries, 95 generals, 0 SOSs, and 0 indicators. Presolve time = 0.02 sec. (7.63 ticks) Probing fixed 160 vars, tightened 80 bounds. Probing changed sense of 40 constraints. Probing time = 0.00 sec. (6.62 ticks) Tried aggregator 2 times. MIP Presolve eliminated 320 rows and 160 columns. MIP Presolve modified 376 coefficients. Aggregator did 200 substitutions. Reduced MIP has 530 rows, 295 columns, and 1860 nonzeros. Reduced MIP has 200 binaries, 95 generals, 0 SOSs, and 0 indicators. Presolve time = 0.01 sec. (2.63 ticks) Probing time = 0.00 sec. (0.53 ticks) Tried aggregator 1 time. Reduced MIP has 530 rows, 295 columns, and 1860 nonzeros. Reduced MIP has 200 binaries, 95 generals, 0 SOSs, and 0 indicators. Presolve time = 0.00 sec. (0.85 ticks) Probing changed sense of 10 constraints. Probing time = 0.00 sec. (0.60 ticks) Clique table members: 310. MIP emphasis: balance optimality and feasibility. MIP search method: dynamic search. Parallel mode: deterministic, using up to 8 threads. Root relaxation solution time = 0.00 sec. (1.22 ticks) Nodes Cuts/ Node Left Objective IInf Best Integer Best Bound ItCnt Gap 0 0 55.0000 52 55.0000 59 0 0 49.0000 80 Cuts: 88 170 0 0 49.0000 64 Cuts: 104 217 0 0 49.0000 64 MIRcuts: 1 218 * 0+ 0 39.0000 49.0000 218 25.64% 0 2 49.0000 64 39.0000 49.0000 218 25.64% Elapsed time = 0.17 sec. (59.33 ticks, tree = 0.01 MB, solutions = 1) * 1155 462 integral 0 40.0000 49.0000 16771 22.50% * 4782+ 1764 43.0000 49.0000 73723 13.95% 5076 1836 48.8125 37 43.0000 49.0000 77941 13.95% * 5190+ 648 45.0000 49.0000 80331 8.89% 5623 126 47.9286 40 45.0000 49.0000 90495 8.89% Clique cuts applied: 10 Cover cuts applied: 22 Implied bound cuts applied: 10 Flow cuts applied: 3 Mixed integer rounding cuts applied: 63 Gomory fractional cuts applied: 12 Root node processing (before b&c): Real time = 0.16 sec. (59.06 ticks) Parallel b&c, 8 threads: Real time = 2.26 sec. (751.26 ticks) Sync time (average) = 0.33 sec. Wait time (average) = 0.00 sec. ------------ Total (root+branch&cut) = 2.42 sec. (810.32 ticks) Result: 45 55 55 33 33 1 22 22 44 44 |