Consider a grid with 40,000 cells. We have customers located on this grid (this is data). We want to place stores on the grid such that each customer has at least one store within reach. Finally, we want to minimize the number of stores we need to use.
Furthermore:
so in a situation like mine I could easily end up with a linear system with millions of rows and columns, which with integer/binary variables will take about the age of the universe to solve.
Well, let's see...
Data
The first thing we do is generating some data. We introduce a grid of \(200 \times 200\) with \(40,000\) cells. Then we randomly place \(M=500\) customers on this grid.
---- 22 PARAMETER cloc customer locations x y cust1 35 75 cust2 169 84 cust3 111 18 cust4 61 163 cust5 59 102 cust6 45 147 cust7 70 165 cust8 172 83 cust9 14 185 cust10 101 79 . . . cust495 83 186 cust496 174 159 cust497 196 148 cust498 115 136 cust499 63 101 cust500 92 87
Reach
Secondly, we create a large, sparse boolean parameter \(\mathit{reach}_{c,i,j}\) indicating if cell \((i,j)\) is within reach of customer \(c\). I used the rule: if the Manhattan distance between \(c\) and \((i,j)\) is 10 or less, we are within reach: \[ |i-\mathit{cloc}_{c,x}|+|j-\mathit{cloc}_{c,y}| \le 10\] This generates a large parameter. A fragment can look like:
---- 29 PARAMETER reach location (i,j) is within reach of customer j1 j2 j3 j4 j5 j6 j7 j8 j9 cust3 .i110 1 cust3 .i111 1 1 cust3 .i112 1 cust39 .i116 1 cust39 .i117 1 1 1 cust39 .i118 1 1 1 1 1 cust39 .i119 1 1 1 1 1 1 cust39 .i120 1 1 1 1 1 1 1 cust39 .i121 1 1 1 1 1 1 1 1
There are 106,376 elements in this data structure.
MIP model
With this, we can formulate the set covering model:
| Set covering model |
|---|
| \[\begin{align}\min\>& \color{darkred}{\mathit{numStores}} \\ & \color{darkred}{\mathit{numStores}} = \sum_{i,j} \color{darkred}{\mathit{placeStore}}_{i,j}\\ & \sum_{i,j} \color{darkblue}{\mathit{reach}}_{c,i,j} \cdot \color{darkred}{\mathit{placeStore}}_{i,j} \ge 1 & \forall c \\ &\color{darkred}{\mathit{placeStore}}_{i,j} \in \{0,1\} \end{align}\] |
This is a large, but easy MIP model. The size is large:
MODEL STATISTICS BLOCKS OF EQUATIONS 2 SINGLE EQUATIONS 501 BLOCKS OF VARIABLES 2 SINGLE VARIABLES 40,001 NON ZERO ELEMENTS 146,377 DISCRETE VARIABLES 40,000
But it solves very fast:
Found incumbent of value 495.000000 after 0.00 sec. (1.89 ticks) Tried aggregator 8 times. MIP Presolve eliminated 306 rows and 39764 columns. MIP Presolve modified 6 coefficients. Aggregator did 29 substitutions. Reduced MIP has 166 rows, 208 columns, and 658 nonzeros. Reduced MIP has 208 binaries, 0 generals, 0 SOSs, and 0 indicators. Presolve time = 0.36 sec. (72.59 ticks) Probing time = 0.00 sec. (0.05 ticks) Tried aggregator 1 time. Detecting symmetries... MIP Presolve eliminated 1 rows and 5 columns. Reduced MIP has 165 rows, 203 columns, and 644 nonzeros. Reduced MIP has 203 binaries, 0 generals, 0 SOSs, and 0 indicators. Presolve time = 0.00 sec. (0.42 ticks) Probing time = 0.00 sec. (0.05 ticks) Clique table members: 35. MIP emphasis: balance optimality and feasibility. MIP search method: dynamic search. Parallel mode: deterministic, using up to 8 threads. Parallel mode: deterministic, using up to 2 threads for parallel tasks at root LP. Tried aggregator 1 time. No LP presolve or aggregator reductions. Presolve time = 0.00 sec. (0.09 ticks) Iteration log . . . Iteration: 1 Dual objective = 57.000000 Iteration: 62 Dual objective = 93.000000 Iteration: 124 Dual objective = 109.000000 Iteration: 186 Dual objective = 111.500000 Root relaxation solution time = 0.06 sec. (0.94 ticks) Nodes Cuts/ Node Left Objective IInf Best Integer Best Bound ItCnt Gap * 0+ 0 495.0000 56.0000 88.69% Found incumbent of value 495.000000 after 0.44 sec. (77.77 ticks) * 0+ 0 118.0000 56.0000 52.54% Found incumbent of value 118.000000 after 0.44 sec. (78.15 ticks) 0 0 111.7778 78 118.0000 111.7778 203 5.27% * 0+ 0 116.0000 111.7778 3.64% Found incumbent of value 116.000000 after 0.45 sec. (79.32 ticks) 0 0 112.5000 78 116.0000 Cuts: 10 209 3.02% * 0+ 0 113.0000 112.5000 0.44% Found incumbent of value 113.000000 after 0.52 sec. (85.02 ticks) 0 0 cutoff 113.0000 112.5000 209 0.44% Elapsed time = 0.55 sec. (85.03 ticks, tree = 0.01 MB, solutions = 3) Zero-half cuts applied: 5 Lift and project cuts applied: 2 Gomory fractional cuts applied: 3 Root node processing (before b&c): Real time = 0.55 sec. (86.90 ticks) Parallel b&c, 8 threads: Real time = 0.00 sec. (0.00 ticks) Sync time (average) = 0.00 sec. Wait time (average) = 0.00 sec. ------------ Total (root+branch&cut) = 0.55 sec. (86.90 ticks) MIP status(101): integer optimal solution
We can see that the presolver is especially successful here: the model size is reduced by an extremely large amount. This will be less pronounced when the number of customers increases.
Solution
The solution looks like:
---- 51 VARIABLE numStores.L = 113 number of stores ---- 51 VARIABLE placeStore.L store locations j1 j6 j7 j8 j9 j15 j16 j17 j18 i4 1 i18 1 i40 1 i70 1 i79 1 i80 1 i107 1 i118 1 i136 1 i157 1 i167 1 i193 1 + j21 j23 j26 j28 j29 j31 j32 j36 j38 i10 1 i28 1 i54 1 i72 1 i96 1 i113 1 i147 1 i158 1 i179 1 i184 1 i198 1 + j39 j44 j45 j46 j49 j50 j56 j58 j59 i5 1 i18 1 i39 1 i62 1 i85 1 i102 1 i104 1 i133 1 i166 1 i195 1 + j62 j66 j67 j68 j69 j73 j74 j76 j80 i11 1 i16 1 i36 1 i61 1 i76 1 i105 1 i112 1 i117 1 i128 1 i146 1 i190 1 + j82 j84 j85 j88 j90 j92 j95 j96 j97 i17 1 i26 1 i35 1 i48 1 i68 1 i79 1 i97 1 i136 1 i156 1 i170 1 i183 1 i191 1 + j98 j102 j107 j111 j112 j114 j115 j116 j118 i4 1 i22 1 i36 1 i56 1 i63 1 i68 1 i88 1 i100 1 i101 1 i111 1 i129 1 i140 1 + j119 j121 j126 j127 j132 j133 j134 j136 j139 i11 1 i30 1 i53 1 i72 1 i111 1 i129 1 i144 1 i159 1 i183 1 i191 1 + j140 j147 j149 j150 j152 j153 j154 j156 j158 i14 1 i35 1 i48 1 i83 1 i98 1 i117 1 i158 1 i174 1 i194 1 + j161 j162 j163 j164 j166 j170 j172 j174 j175 i5 1 i32 1 i42 1 i61 1 i69 1 i103 1 i143 1 i145 1 i158 1 i192 1 i198 1 + j176 j178 j179 j180 j182 j183 j184 j188 j191 i6 1 i13 1 i23 1 i47 1 i61 1 i81 1 i93 1 i103 1 i125 1 i182 1 i193 1 + j192 j193 j196 i73 1 i120 1 i138 1 i167 1
We have really debunked that this is an intractable problem for a good MIP solver.
Note that the age of the universe is 13.7 billion years or 4.3e17 seconds. So we have achieved a speed-up of about 1e17. This is a record for me.
It is noted that this model does not really answer the question: what are the best locations for my stores? It only answers: how many stores do we need? In particular, it does not try to place stores as close as possible to customers. A subsequent optimization model can find good (or optimal) locations for the stores. Reference [2] discusses a similar issue for the continuous case, where facilities can be placed in any \((x,y)\) location instead of on a grid \((i,j)\).
Conclusion
We can solve this problem for this particular problem quite efficiently using a standard MIP solver. Quite often I see that these solution approaches are dismissed, due to a misunderstanding of complexity-theoretic arguments. Complexity theory says very little about the performance of a particular MIP model using a particular solver on a given data set. It is better not to reject a MIP approach even before trying.
References
- Minimize number of shops while reaching all customers, https://stackoverflow.com/questions/62224042/minimize-number-of-shops-while-reaching-all-customers
- Solving a facility location problem as an MIQP, https://yetanothermathprogrammingconsultant.blogspot.com/2018/01/solving-facility-location-problem-as.html (this is essentially a continuous version of the current problem)
No comments:
Post a Comment