In this post, I consider the following question:
Given \(m\) consumer locations, what is the minimum number of facilities (warehouses) needed, \(n\), so each consumer is within \(K\) kilometers of a warehouse.
A second, related question is:
Where should we place these \(n\) warehouses?
I looked at this earlier in [1]. I want to add some new experiments:
- Do something about symmetry: breaking symmetry by simple ordering constraints. These can be highly effective (but can also disappoint).
- See if I can formulate the "min sum distance" model as a convex SOCP (second-order cone programming) problem.
- Compare with a discrete location problem: place facilities on a grid. This makes the model very large but linear. How does it compare to the quadratic models that place facilities anywhere?
Generate consumer locations.
I just use random locations drawn from a uniform distribution: \[\color{darkblue}p_{i,c} \sim U(0,\color{darkblue}S)\] where \(c \in \{x,y\}\) are coordinates and \(\color{darkblue}S=100\). Of course, this is not very realistic.
Q: What would be a good way to generate more realistic locations? Typically, locations are clustered a bit. Maybe: uniform locations with some random count of consumers at each location.
For my experiments, I used a dataset \(m=50\) consumer locations. This is both small enough to make the models solve quickly and large enough to show differences in performance.
---- 37 PARAMETER p points x y p1 17.175 84.327 p2 55.038 30.114 p3 29.221 22.405 p4 34.983 85.627 p5 6.711 50.021 p6 99.812 57.873 p7 99.113 76.225 p8 13.069 63.972 p9 15.952 25.008 p10 66.893 43.536 p11 35.970 35.144 p12 13.149 15.010 p13 58.911 83.089 p14 23.082 66.573 p15 77.586 30.366 p16 11.049 50.238 p17 16.017 87.246 p18 26.511 28.581 p19 59.396 72.272 p20 62.825 46.380 p21 41.331 11.770 p22 31.421 4.655 p23 33.855 18.210 p24 64.573 56.075 p25 76.996 29.781 p26 66.111 75.582 p27 62.745 28.386 p28 8.642 10.251 p29 64.125 54.531 p30 3.152 79.236 p31 7.277 17.566 p32 52.563 75.021 p33 17.812 3.414 p34 58.513 62.123 p35 38.936 35.871 p36 24.303 24.642 p37 13.050 93.345 p38 37.994 78.340 p39 30.003 12.548 p40 74.887 6.923 p41 20.202 0.507 p42 26.961 49.985 p43 15.129 17.417 p44 33.064 31.691 p45 32.209 96.398 p46 99.360 36.990 p47 37.289 77.198 p48 39.668 91.310 p49 11.958 73.548 p50 5.542 57.630
Model 1: Minimize the number of facilities
In the following, \(i\) indicates the consumer locations, and \(j\) refers to warehouses (facilities). The number of possible warehouses we can use is encoded in the set \(j \in J\). It should be larger than the number we expect, but not crazy large as this makes the problem much larger. If you choose it too small, the problem will be infeasible.
| MIQCP model 1: min num facilities | |
|---|---|
| Data | \[\begin{align}&i &&\text{consumers, customers}\\ & j &&\text{facilities, warehouses (overestimate number)} \\ & c=\{x,y\} && \text{use 2D coordinates}\\ & \color{darkblue}p_{i,c} && \text{locations of customers} \\ & \color{darkblue}K && \text{maximum distance} \\ & \color{darkblue} S = 100 && \text{side of square we are considering} \\ & \color{darkblue}M_1= 2\color{darkblue}S^2 - \color{darkblue}K^2&& \text{big-M value}\end{align}\] |
| Variables | \[\begin{align} &\color{darkred}w_{j,c} \in [0,\color{darkblue}S] && \text{location of facility}\\ &\color{darkred}a_{i,j} \in \{0,1\} && \text{assignment of consumer to facility}\\ & \color{darkred}u_j \in \{0,1\} && \text{facility is used} \\ & \color{darkred} n && \text{number of facilities needed}\end{align}\] |
| Equations | \[\begin{align}\min\>&\color{darkred}n = \sum_j \color{darkred} u_j \\ & \sum_c \left( \color{darkblue}p_{i,c} - \color{darkred}w_{j,c}\right)^2 \le \color{darkblue}K^2 + \color{darkblue}M_1 (1- \color{darkred}a_{i,j}) && \forall i,j && \text{distance}\\ & \sum_j \color{darkred}a_{i,j} = 1 && \forall i && \text{assignment}\\ & \color{darkred}a_{i,j} \le \color{darkred}u_j && \forall i,j && \text{closed facility} \\ & \color{darkred}u_j \ge \color{darkred}u_{j+1} && && \text{ordering (optional)} \\ & \color{darkred}w_{j,x} \le \color{darkred}w_{j+1,x}&& && \text{ordering (optional)} \end{align} \] |
We restricted the search area to \([0,\color{darkblue}S]^2\). Assuming facilities will be placed within the convex hull of customer points, we could as well require: \[ \min_i \color{darkblue}p_{i,c}\le \color{darkred}w_{j,c} \le \max_i \color{darkblue}p_{i,c}\] For simplicity, I chose \(\color{darkred}w_{j,c} \in [0,\color{darkblue}S]\).
The distance equation says: for consumer \(i\) there is (at least) one facility \(j\) that is closer than \(K\) km. We do this by assigning a warehouse to each consumer and imposing the distance constraint for that assigned pair. A general nonlinear distance constraint could look like: \[\color{darkred}a_{i,j}=1 \Rightarrow \sqrt{\sum_c \left( \color{darkblue}p_{i,c} - \color{darkred}w_{j,c}\right)^2} \le \color{darkblue} K\] This could be implemented as a quadratic indicator constraint: \[\color{darkred}a_{i,j}=1 \Rightarrow \sum_c \left( \color{darkblue}p_{i,c} - \color{darkred}w_{j,c}\right)^2 \le \color{darkblue} K^2\] Interestingly, solvers that support indicator constraints, seem to allow only linear constraints on the right. (Why? I don't know) Here we have a quadratic constraint, so we have to reformulate this as a big-M constraint \[ \sum_c \left( \color{darkblue}p_{i,c} - \color{darkred}w_{j,c}\right)^2 \le \color{darkblue}K^2 +\color{darkblue}M_1 (1-\color{darkred}a_{i,j})\] A conservative bound gives us the following value \(\color{darkblue}M_1 = 2\color{darkblue}S^2 - \color{darkblue}K^2\).
This is a convex quadratic constraint, so we can use well-known solvers for this.
We have two ordering constraints. The first one says: the used warehouses are the first ones. I.e. \(u_j\) will look like \([1,1,1,0,0,0,0]\). This actually helped me in the implementation of the subsequent models. The first \(n\) should be reused in those models. But it also imposes an ordering that can help against a lot of symmetry. To further reduce symmetry we also order the facilities by their \(x\)-coordinate.
The addition of these symmetry-breaking constraints is not a guaranteed success. We make the search space much smaller, so this is a big advantage. But we also may make finding a feasible solution much more difficult. Before doing some experimentation, it is very hard to predict how these trade-offs will play out.
The solution can look like:
We see we need 3 facilities, but the locations and the assignments are not optimal. We will use a second model to address this.
Model 2: Minimize the sum of the squared distances
| MIQCP model 2: min sum squared distances | |
|---|---|
| Data | \[\begin{align}& \color{blue}n && \text{number of facilities needed}\\ &i &&\text{consumers, customers}\\ & j = \{1,\dots,\color{darkblue}n\}&&\text{facilities, warehouses} \\ & c=\{x,y\} && \text{use 2D coordinates}\\ & \color{darkblue}p_{i,c} && \text{locations of customers} \\ & \color{darkblue}K && \text{maximum distance} \\ & \color{darkblue}M_2 = 2\color{darkblue}S^2&& \text{big-M value} \end{align}\] |
| Variables | \[\begin{align} &\color{darkred}\Delta_{i,j}\in [0,\color{darkblue}K^2]&& \text{squared distance between assigned facility and customer (or zero)} \\ &\color{darkred}w_{j,c} \in [0,\color{darkblue}S] && \text{location of facility}\\ &\color{darkred}a_{i,j} \in \{0,1\} && \text{assignment of consumer to facility} \end{align}\] |
| Equations | \[\begin{align}\min\>&\sum_{i,j} \color{darkred} \Delta_{i,j} \\ & \color{darkred}\Delta_{i,j} \ge \sum_c \left( \color{darkblue}p_{i,c} - \color{darkred}w_{j,c}\right)^2 - \color{darkblue}M_2 (1- \color{darkred}a_{i,j}) && \forall i,j && \text{squared distance}\\ & \sum_j \color{darkred}a_{i,j} = 1 && \forall i && \text{assignment}\\ & \color{darkred}w_{j,x} \le \color{darkred}w_{j+1,x}&& && \text{ordering (optional)} \end{align} \] |
The distance constraint implements the implication: \[\color{darkred}a_{i,j} = 1 \Rightarrow \color{darkred}\Delta_{i,j} \ge \sum_c \left( \color{darkblue}p_{i,c} - \color{darkred}w_{j,c}\right)^2\] Again, indicator constraints are only for linear constraints, so we need a big-M formulation. A conservative value for the big-M constant is \(2\color{darkblue}M_2 = \color{darkblue}S^2\). This will allow \(\color{darkred}\Delta_{i,j}=0\) when \(\color{darkred}a_{i,j} = 0\).
The result can be:
This picture is much better. We see that the location of warehouses and the assignments make much more sense.
We can look at this as a multi-objective problem where we have the objectives (1) minimize the number of facilities and (2) minimize the sum of (squared) distances. Here we used a lexicographic approach: solve for the first objective and add this as a "constraint" to the second objective problem. In this case, the constraint turns out just a question of dimensioning the second problem.
Model 3: Minimize the sum of the distances
It looks like converting model 2 from minimizing the sum of squared distances to the sum of distances should be an easy task. However, this is not exactly the case. A first attempt leads to an MINLP model:
| MINLP model 3a: min sum distances | |
|---|---|
| Equations | \[\begin{align}\min\>&\sum_{i,j} \sqrt{\color{darkred} \Delta_{i,j}} \\ & \color{darkred}\Delta_{i,j} \ge \sum_c \left( \color{darkblue}p_{i,c} - \color{darkred}w_{j,c}\right)^2 - \color{darkblue}M_2 (1- \color{darkred}a_{i,j}) && \forall i,j && \text{squared distance}\\ & \sum_j \color{darkred}a_{i,j} = 1 && \forall i && \text{assignment}\\ & \color{darkred}w_{j,x} \le \color{darkred}w_{j+1,x}&& && \text{ordering (optional)} \end{align} \] |
The data and the variables are the same as for model 2. The objective makes this problem a general MINLP model which is more difficult to solve and requires different solvers.
| Non-convex MIQCP model 3b: min sum distances | |
|---|---|
| Variables | \[\begin{align} &\color{darkred}d_{i,j}\in [0,\color{darkblue}K]&& \text{distance between assigned facility and customer (or zero)} \\ &\color{darkred}w_{j,c} \in [0,\color{darkblue}S] && \text{location of facility}\\ &\color{darkred}a_{i,j} \in \{0,1\} && \text{assignment of consumer to facility} \end{align}\] |
| Equations | \[\begin{align}\min\>&\sum_{i,j} \color{darkred} d_{i,j} \\ & \color{darkred}d^2_{i,j} \ge \sum_c \left( \color{darkblue}p_{i,c} - \color{darkred}w_{j,c}\right)^2 - \color{darkblue}M_2 (1- \color{darkred}a_{i,j}) && \forall i,j && \text{squared distance}\\ & \sum_j \color{darkred}a_{i,j} = 1 && \forall i && \text{assignment}\\ & \color{darkred}w_{j,x} \le \color{darkred}w_{j+1,x}&& && \text{ordering (optional)} \end{align} \] |
| MISOCP Model 3c: min sum distances | |
|---|---|
| Variables | \[\begin{align} &\color{darkred}z_{i,j}\in [0,\color{darkblue}K]&& \text{distance between assigned facility and customer (or zero)} \\ &\color{darkred}d_{i,j}\in [0,\sqrt{\color{darkblue}M_2}]&& \text{distance between any facility and customer}\\ & \color{darkred} y_{i,j,c} && \text{intermediate free variable} \\ &\color{darkred}w_{j,c} \in [0,\color{darkblue}S] && \text{location of facility}\\ &\color{darkred}a_{i,j} \in \{0,1\} && \text{assignment of consumer to facility} \end{align}\] |
| Equations | \[\begin{align}\min\>&\sum_{i,j} \color{darkred} z_{i,j} \\ & \color{darkred}z_{i,j} \ge \color{darkred}d_{i,j} - \sqrt{\color{darkblue} M_2 }(1-\color{darkred}a_{i,j}) && \forall i,j && \text{distance is zero when not assigned}\\ & \color{darkred}d^2_{i,j} \ge \sum_c \color{darkred}y^2_{i,j,c} && \forall i,j && \text{SOCP constraint} \\ & \color{darkred}y_{i,j,c} = \color{darkblue}p_{i,c} - \color{darkred}w_{j,c} && \forall i,j,c && \text{intermediate}\\ & \sum_j \color{darkred}a_{i,j} = 1 && \forall i && \text{assignment}\\ & \color{darkred}w_{j,x} \le \color{darkred}w_{j+1,x}&& && \text{ordering (optional)} \end{align} \] |
We needed to add a lot of extra variables and constraints to make this model work for an SOCP solver.
Model 4: grid search
| MIP Model 4: locations on a grid | |
|---|---|
| Data | \[\begin{align}&g&&\text{grid points}\\ &\color{darkblue}{\mathit{dist}}_{i,g}&&\text{distances}\\ &\color{darkblue}{\mathit{ok}}_{i,g}&&\text{admissible assignments} \end{align}\] |
| Variables | \[\begin{align}&\color{darkred}{\mathit{place}}_g \in \{0,1\}&&\text{whether to place a facility at $g$}\\ &\color{darkred}{\mathit{ag}}_{i,g} \in \{0,1\}&&\text{assign customer $i$ to $g$} \end{align}\] |
| Model 4a | \[\begin{align}\min\>&\color{darkred}n\\&\sum_{g|\color{darkblue}{\mathit{ok}}(i,g)} \color{darkred}{\mathit{ag}}_{i,g} = 1 &&\forall i &&\text{assignment} \\ & \color{darkred}{\mathit{ag}}_{i,g} \le \color{darkred}{\mathit{place}}_g &&\forall i,g|\color{darkblue}{\mathit{ok}}(i,g) &&\text{use location}\\ &\color{darkred}n = \sum_g \color{darkred}{\mathit{place}}_g && && \text{number of facilities}\end{align}\] |
| Model 4b | \[\begin{align} \min & \sum_{i,g|\color{darkblue}{\mathit{ok}}(i,g)} \color{darkblue}{\mathit{dist}}_{i,g} \cdot \color{darkred}{\mathit{ag}}_{i,g} && && \text{sum distances} \\ & \sum_{g|\color{darkblue}{\mathit{ok}}(i,g)} \color{darkred}{\mathit{ag}}_{i,g} = 1 &&\forall i &&\text{assignment} \\ & \color{darkred}{\mathit{ag}}_{i,g} \le \color{darkred}{\mathit{place}}_g &&\forall i,g|\color{darkblue}{\mathit{ok}}(i,g) &&\text{use location}\\ & \sum_g \color{darkred}{\mathit{place}}_g = \color{darkblue}n^* && && \text{fix number of facilities} \end{align} \] |
The pictures look very much the same.
Performance
| Model | Type | Objective | Sum dist^2 | Sum dist | Time | Nodes | Iterations |
|---|---|---|---|---|---|---|---|
| 1 | MIQCP | min num fac | 42,127 | 1,380 | 20.27 | 15,487 | 490,331 |
| 1 + order | MIQCP | min num fac | 39,288 | 1,313 | 6.375 | 1,230 | 28,284 |
| 2 | MIQCP | min sum dist^2 | 21,105 | 937.8 | 49.23 | 255,858 | 976,174 |
| 2 + order | MIQCP | min sum dist^2 | 21,105 | 937.8 | 15.31 | 87,319 | 328,474 |
| 3 | MISOCP | min sum dist | 21,186 | 935.1 | 2,981 | 5,896,361 | 66,465,832 |
| 3 + order | MISOCP | min sum dist | 21,186 | 935.1 | 431.9 | 933,002 | 14,612,297 |
| 4a | MIP | min num fac | 38,134 | 1,273 | 31.17 | 0 | 21 |
| 4b | MIP | min sum dist | 21,313 | 944.7 | 15.94 | 0 | 17 |
A lot to digest here. But what we see is:
- Adding the ordering constraints to help with symmetry does indeed help. Sometimes a lot.
- The MIP model is super easy to solve. Of course, the problem is reduced to more or less an assignment problem. This is something solvers can do very well. In this case, the model can be solved completely during preprocessing. But note that it is very big.
- The second-order cone model is much more difficult to solve than the MIQCP model that minimizes the sum of the squared distances. My intuition: too many epigraph-type indirections in the model so the pricing signals from the objective are less obvious.
Conclusions
- Simple symmetry-breaking constraints that enforce an ordering on the decision variables can make a large difference.
- Indicator constraints are for linear constraints only. So we had to implement big-M formulations.
- Formulating SOCPs can require some real effort. Tools (modeling systems, solvers) can do a better job in applying appropriate transformations automatically.
- Experimentation with different formulations can pay off. Here we see that min sum distance^2 is preferable above min sum distance: much easier to solve and almost the same optimal solution. Not obvious (to me) before implementing and trying both formulations.
- It may be better to use a grid of candidate locations. This way distances can be precomputed, and we have a linear problem. However, the model can become very big if the grid is too fine.
- Just a simple question: what is the minimum number of warehouses needed so each customer is within a given distance? Not so simple to answer. There are quite a few angles to explore.
References
- Solving a facility location problem as a MIQCP, https://yetanothermathprogrammingconsultant.blogspot.com/2018/01/solving-facility-location-problem-as.html
- Second-order cone programming, https://en.wikipedia.org/wiki/Second-order_cone_programming
I sometimes wonder why solvers don't provide built-in support to compare the performance of different problem formulations. Often, the large number of possible formulations makes it difficult to compare by hand.
ReplyDelete