I am trying to solve an uncapacitated facility location problem in which I would like to find the location to open at most \(k\) facilities, and assign \(n\) clients to those facilities, each with cost \(c_i\) (proportionate to the distance of the client from the facility it's assigned to) for assignment with a budget \(B\). So I would like to find the locations of the facilities, and assign clients to those facilities such that the maximum number of clients are serviced within the budget \(B\). I haven't found a variant like this in which there is an overall budget; I've only seen versions where we want to minimize the overall cost while assigning all the clients to facilities.
Let's give this a shot.
Example Data
A good way to familiarize yourself with a problem, is to invent some example data. Here is some random test data I used for the models below:
---- 35 SET i clients client1 , client2 , client3 , client4 , client5 , client6 , client7 , client8 , client9 client10, client11, client12, client13, client14, client15, client16, client17, client18 client19, client20, client21, client22, client23, client24, client25 ---- 35 SET j facilities facility1, facility2, facility3, facility4 ---- 35 SET k xy-coordinates x, y ---- 35 PARAMETER U = 100.000 box size for generating random data PARAMETER maxdist = 141.421 maximum distance (big-M) PARAMETER budget = 500.000 available budget ---- 35 PARAMETER client location and cost data x y cost client1 17.175 84.327 6.611 client2 55.038 30.114 7.558 client3 29.221 22.405 6.274 client4 34.983 85.627 2.839 client5 6.711 50.021 0.864 client6 99.812 57.873 1.025 client7 99.113 76.225 6.413 client8 13.069 63.972 5.453 client9 15.952 25.008 0.315 client10 66.893 43.536 7.924 client11 35.970 35.144 0.728 client12 13.149 15.010 1.757 client13 58.911 83.089 5.256 client14 23.082 66.573 7.502 client15 77.586 30.366 1.781 client16 11.049 50.238 0.341 client17 16.017 87.246 5.851 client18 26.511 28.581 6.212 client19 59.396 72.272 3.894 client20 62.825 46.380 3.587 client21 41.331 11.770 2.430 client22 31.421 4.655 2.464 client23 33.855 18.210 1.305 client24 64.573 56.075 9.334 client25 76.996 29.781 3.799
The clients are uniformly distributed over the \([0,100]\times[0,100]\) square.
 |
| Client locations |
1. High-level MINLP Formulation
The first step is to translate the "word-problem" into a high-level Mixed-Integer Nonlinear Programming (MINLP) model. The goal is not to achieve the best performance possible, but rather to have a "reference model": a precise mathematical description of the problem.
We use the following naming:
We use the following naming:
| Names | |
|---|---|
| Indices | |
| \(i\) | Clients |
| \(j\) | Facilities |
| \(k\) | x or y coordinate |
| Data | |
| \(\color{darkblue}{\mathit{Cost}}_i\) | Unit cost |
| \(\color{darkblue}{\mathit{Client}}_{i,k}\) | Client locations |
| \(\color{darkblue}{\mathit{Budget}}\) | Available budget |
| \(\color{darkblue}U\) | Size of box containing client locations |
| Variables | |
| \(\color{darkred}{\mathit{Assign}}_{i,j} \in \{0,1\} \) | Assignment of client to facility |
| \(\color{darkred}{\mathit{ClientUsed}}_{i} \in \{0,1\} \) | Client is assigned to a facility |
| \(\color{darkred}{\mathit{Facility}}_{j,k} \in [0,\color{darkblue}U]\) | Facility locations |
| \(\color{darkred}{\mathit{Dist}}_{i,j} \in [0,\sqrt{2}\color{darkblue}U] \) | Distance between client and facility Zero if not assigned |
With this we can start developing our model:
| MINLP Problem |
|---|
| \[\begin{align}\max & \sum_i \color{darkred}{\mathit{ClientUsed}}_i\\ & \color{darkred}{\mathit{ClientUsed}}_i = \sum_j \color{darkred}{\mathit{Assign}}_{i,j} \\ &\color{darkred}{\mathit{Dist}}_{i,j} = \color{darkred}{\mathit{Assign}}_{i,j} \cdot \sqrt{\sum_k \left(\color{darkblue}{\mathit{Client}}_{i,k}-\color{darkred}{\mathit{Facility}}_{j,k} \right)^2} \\ & \sum_{i,j} \color{darkblue}{\mathit{Cost}}_i \cdot \color{darkred}{\mathit{Dist}}_{i,j} \le \color{darkblue}{\mathit{Budget}}\\ &\color{darkred}{\mathit{ClientUsed}}_{i}, \color{darkred}{\mathit{Assign}}_{i,j} \in \{0,1\} \\ & \color{darkred}{\mathit{Facility}}_{j,k} \in [0,\color{darkblue}U]\\ & \color{darkred}{\mathit{Dist}}_{i,j} \in [0,\sqrt{2}\color{darkblue}U]\end{align}\] |
We don't need to add an assignment constraint like\[\sum_j {\mathit{Assign}}_{i,j} \le 1 \>\>\forall i \] as this is already captured in the constraint \[\mathit{ClientUsed}_i = \sum_j {\mathit{Assign}}_{i,j}\]
The distance constraint can be interpreted as \[{\mathit{Dist}}_{i,j} = \begin{cases} \sqrt{\displaystyle\sum_k \left({\mathit{Client}}_{i,k}-{\mathit{Facility}}_{j,k} \right)^2} & \text{if $\mathit{Assign}_{i,j}=1$} \\ 0 & \text{if $\mathit{Assign}_{i,j}=0$}\end{cases}\]
The distance constraint can be interpreted as \[{\mathit{Dist}}_{i,j} = \begin{cases} \sqrt{\displaystyle\sum_k \left({\mathit{Client}}_{i,k}-{\mathit{Facility}}_{j,k} \right)^2} & \text{if $\mathit{Assign}_{i,j}=1$} \\ 0 & \text{if $\mathit{Assign}_{i,j}=0$}\end{cases}\]
We can even use a global MINLP solver to get optimal solutions for smaller data sets (or "good" solutions for slightly larger problems). This will give us more confidence that the model is correctly representing the problem.
2. Mixed Integer Quadratically Constrained Model
The second step is to get rid of some of the non-linearities so that a quadratic model remains. We need to work on the distance constraint. By making it a big-M constraint, we can make everything quadratic. We get rid of the square root and the multiplication by the Assign variable.
| MIQCP Problem |
|---|
| \[\begin{align}\max & \sum_i \color{darkred}{\mathit{ClientUsed}}_i\\ & \color{darkred}{\mathit{ClientUsed}}_i = \sum_j \color{darkred}{\mathit{Assign}}_{i,j} \\ &\color{darkred}{\mathit{Dist}}_{i,j}^2 \ge \sum_k \left(\color{darkblue}{\mathit{Client}}_{i,k}-\color{darkred}{\mathit{Facility}}_{j,k} \right)^2 - \color{darkblue}M (1-\color{darkred}{\mathit{Assign}}_{i,j}) \\ & \sum_{i,j} \color{darkblue}{\mathit{Cost}}_i \cdot \color{darkred}{\mathit{Dist}}_{i,j} \le \color{darkblue}{\mathit{Budget}} \\ &\color{darkred}{\mathit{ClientUsed}}_{i}, \color{darkred}{\mathit{Assign}}_{i,j} \in \{0,1\}\\ & \color{darkred}{\mathit{Facility}}_{j,k} \in [0,\color{darkblue}U]\\ & \color{darkred}{\mathit{Dist}}_{i,j} \in [0,\sqrt{2}\color{darkblue}U]\\&\color{darkblue} M = 2\color{darkblue}U^2 \end{align}\] |
In this model we have only a bound on the variable \(\mathit{Dist}_{i,j}\). That means it can be a bit larger than needed, but only as long as it obeys the Budget constraint. If we want to report on the actual distances, we need to recalculate them.
If we run this with Cplex we see:
CPLEX Error 5002: 'QCP_row_for_distance2(client1.facility1)' is not convex.
The new version Gurobi 9 can solve this as a nonconvex MIQCP. Of course we can also solve this with a global MINLP solver.
The problem would be easier if the costs do not increase proportionally with the distance but with the squared distance. That would yield a simpler and convex MIQCP problem. [2]
3. Mixed Integer Second Order Cone Formulation
The constraint \[y^2 \ge \sum_i x_i^2 \>\>\text{where $y\ge 0$}\] is a Second Order Cone constraint. This type of constraint is supported by Cplex. But only if this is specified (almost) exactly like this. I.e. no additional linear terms in the constraints. The mechanics are not complicated, just a bit tedious. We need to split our distance constraint in parts to shoehorn it into a SOCP framework:
| MISOCP Problem |
|---|
| \[\begin{align}\max & \sum_i \color{darkred}{\mathit{ClientUsed}}_i\\ & \color{darkred}{\mathit{ClientUsed}}_i = \sum_j \color{darkred}{\mathit{Assign}}_{i,j} \\ & \color{darkred}{\mathit{Dist}}_{i,j} \ge \color{darkred}\Delta_{i,j} - \color{darkblue}M (1-\color{darkred}{\mathit{Assign}}_{i,j}) \\ &\color{darkred}{\mathit{\Delta}}_{i,j}^2 \ge \sum_k \color{darkred} D_{i,j,k}^2\\ & \color{darkred} D_{i,j,k} = \color{darkblue}{\mathit{Client}}_{i,k}-\color{darkred}{\mathit{Facility}}_{j,k} \\ & \sum_{i,j} \color{darkblue}{\mathit{Cost}}_i \cdot\color{darkred}{\mathit{Dist}}_{i,j} \le \color{darkblue}{\mathit{Budget}} \\ &\color{darkred}{\mathit{ClientUsed}}_{i}, \color{darkred}{\mathit{Assign}}_{i,j} \in \{0,1\}\\ & \color{darkred}{\mathit{Facility}}_{j,k} \in [0,\color{darkblue}U]\\ & \color{darkred}{\mathit{Dist}}_{i,j}, \color{darkred}\Delta_{i,j} \in [0,\sqrt{2}\color{darkblue}U]\\ & \color{darkred}D_{i,j,k} \>\text{free} \\ & \color{darkblue} M = \sqrt{2}\color{darkblue}U\end{align}\] |
This is now a convex problem. It can be solved with MISOCP solvers like Cplex, Gurobi and Mosek.
This is just a reformulation of the previous MIQCP model. Solvers do not typically recognize that form as a second order cone problem. So often the task is left to the modeler to help the solver and produce a convex model.
There are two objections against this special way of writing SOCP constraints:
- We need to have fairly active knowledge of how SOCP constraints can look like. I.e. we need to remember a bag of reformulation tricks. If you don't, you may miss out on solver capabilities.
- Even if we know how to do this, the reformulation is making the model less clear and less intuitive. We need to add extra comments why additional, at first sight unnecessary, variables and constraints were added.
Some high level modeling systems can apply reformulations automatically [3]. I think that is a good idea in order to keep formulations as natural, and as close to the original problem as possible.
This is a small problem. But Cplex has no easy time in solving it. Interestingly, it found the optimal solution completely at the end.
6939104 1303772 20.7395 24 19.0000 20.8047 2.53e+08 9.50% 6971225 1302201 20.0000 20 19.0000 20.7915 2.54e+08 9.43% 7005795 1298106 infeasible 19.0000 20.7779 2.55e+08 9.36% 7039089 1292724 20.0000 24 19.0000 20.7573 2.56e+08 9.25% 7072254 1283943 infeasible 19.0000 20.7430 2.58e+08 9.17% 7107348 1281824 infeasible 19.0000 20.7274 2.58e+08 9.09% 7143017 1262286 infeasible 19.0000 20.7002 2.60e+08 8.95% 7179356 1248370 20.0000 29 19.0000 20.6153 2.61e+08 8.50% 7218255 1222703 infeasible 19.0000 20.2398 2.62e+08 6.53% *7227772 1220913 integral 0 20.0000 20.1541 2.62e+08 0.77% Found incumbent of value 19.999999 after 1932.55 sec. (890057.09 ticks) *7247324 926397 integral 0 20.0000 20.0058 2.63e+08 0.03% Cone: 48 Found incumbent of value 20.000000 after 1940.95 sec. (892272.37 ticks) 7252826 926392 infeasible 20.0000 20.0058 2.63e+08 0.03% Elapsed time = 1945.70 sec. (892962.40 ticks, tree = 766.57 MB, solutions = 28) 7280114 64647 infeasible 20.0000 20.0002 2.64e+08 0.00% 7303087 40831 cutoff 20.0000 20.0000 2.64e+08 0.00% 7323676 15673 infeasible 20.0000 20.0000 2.66e+08 0.00% Cover cuts applied: 1 Mixed integer rounding cuts applied: 20 Gomory fractional cuts applied: 10 Cone linearizations applied: 3328 Root node processing (before b&c): Real time = 0.25 sec. (50.77 ticks) Parallel b&c, 16 threads: Real time = 1995.13 sec. (905620.31 ticks) Sync time (average) = 265.77 sec. Wait time (average) = 0.00 sec. ------------ Total (root+branch&cut) = 1995.38 sec. (905671.08 ticks) MIQCP status(101): integer optimal solution
The solution looks like:
---- 156 VARIABLE assign.L client to facility facility1 facility2 facility3 facility4 client1 1.000 client3 1.000 client4 1.000 client5 1.000 client6 1.000 client8 1.000 client9 1.000 client10 1.000 client11 1.000 client12 1.000 client14 1.000 client15 1.000 client16 1.000 client17 1.000 client18 1.000 client20 1.000 client21 1.000 client22 1.000 client23 1.000 client25 1.000 ---- 156 VARIABLE facility.L coordinates x y facility1 17.146 84.398 facility2 23.055 66.584 facility3 29.228 22.420 facility4 66.892 43.536 ---- 156 VARIABLE clientUsed.L client is assigned to facility client1 1.000, client3 1.000, client4 1.000, client5 1.000 client6 1.000, client8 1.000, client9 1.000, client10 1.000 client11 1.000, client12 1.000, client14 1.000, client15 1.000 client16 1.000, client17 1.000, client18 1.000, client20 1.000 client21 1.000, client22 1.000, client23 1.000, client25 1.000 ---- 156 PARAMETER DistanceReport distances for assignments facility1 facility2 facility3 facility4 client1 0.084 client3 0.029 client4 22.499 client5 23.346 client6 35.973 client8 10.337 client9 13.634 client10 0.020 client11 14.488 client12 17.746 client14 0.047 client15 17.008 client16 20.391 client17 3.082 client18 6.750 client20 4.991 client21 16.154 client22 17.930 client23 6.298 client25 17.089 ---- 156 PARAMETER CostAllocation cost = unit cost * distance for assignment s facility1 facility2 facility3 facility4 client1 0.557 client3 0.180 client4 63.867 client5 20.177 client6 36.877 client8 56.367 client9 4.298 client10 0.155 client11 10.542 client12 31.173 client14 0.352 client15 30.295 client16 6.962 client17 18.036 client18 41.934 client20 17.903 client21 39.260 client22 44.184 client23 8.219 client25 64.929 ---- 156 PARAMETER CostReport total cost vs budget total 496.267, limit 500.000
Basically: we can cover 20 clients with the given budget and 4 facilities. We recalculated the cost, as the optimization model may overestimate the distance and the cost (it only worries about the cost to be below the budget). In the picture below we see the assigned clients (and their corresponding facility) as colored clusters. Remaining unassigned clients are drawn as black points.
 |
| Solution |
One thing we can see is that a small distance does not guarantee to be selected: it can be too costly because of high unit cost. This is visible in the light blue cluster: a far away point is included while a closer by point is not included. We should also note that this is not a min cost solution, just a solution that obeys the budget constraint. We are not minimizing the cost.
Symmetry
In the models above we can reduce symmetry by having the locations \(\mathit{facility}_{j,k}\) ordered by the \(x\)-coordinate. We can do this by adding the (linear) constraint \[\mathit{facility}_{j,x} \le \mathit{facility}_{j+1,x}\] This makes the solution easier to read, but more importantly will remove a lot of symmetry. This can reduce solution times.
4. Approximation: restrict locations to grid points
In the previous models, we calculated two quantities:
- \(x,y\) location of each facility (continuous variables)
- Assignment of clients to a facility (binary variable)
This turned out to be a difficult model to solve. The usual approach is: if a model is too difficult to solve, solve a different problem. If we set up a grid, and allow a facility only to be placed on the grid points, we get a much simpler model. We can pre-compute the distance and cost between each client location and each possible facility location. The end result is just a linear, but large MIP model.
The model looks like:
| MIP Problem |
|---|
| \[\begin{align}\max & \sum_i \color{darkred}{\mathit{ClientUsed}}_i\\ & \color{darkred}{\mathit{ClientUsed}}_i = \sum_p \color{darkred}{\mathit{Assign}}_{i,p} \\ & \sum_p \color{darkred}{\mathit{PointUsed}}_p = \color{darkblue}{\mathit{NumFacilities}}\\ & \color{darkred}{\mathit{Assign}}_{i,p} \le \color{darkblue}{\mathit{NumClients}} \cdot\color{darkred}{\mathit{PointUsed}}_p\\ & \sum_{i,p} \color{darkblue}{\mathit{Cost}}_{i,p} \cdot \color{darkred}{\mathit{Assign}}_{i,p} \le \color{darkblue}{\mathit{Budget}} \\ & \color{darkred}{\mathit{Assign}}_{i,p}, \color{darkred}{\mathit{PointUsed}}_p, \color{darkred}{\mathit{ClientUsed}}_i \in \{0,1\} \end{align}\] |
Here \(p\) are the grid points. As we have box of \([0,100]\times[0,100]\) for the client locations, I used \(101 \times 101 = 10201\) grid points. This is a fine grid, so we expect similar results as for the models where facilities can be located anywhere on the \(x\)-\(y\) plane. The parameter \(\mathit{Cost}_{i,p}\) is precalculated. The resulting model is very large (considering the original data set was small):
MODEL STATISTICS BLOCKS OF EQUATIONS 5 SINGLE EQUATIONS 255,053 BLOCKS OF VARIABLES 4 SINGLE VARIABLES 265,252 NON ZERO ELEMENTS 1,030,352 DISCRETE VARIABLES 265,251
However it solves very fast (less than a minute). The results are:
---- 121 VARIABLE assign.L client to facility gp1803 gp2953 gp6043 gp6811 client1 1.000 client3 1.000 client4 1.000 client5 1.000 client6 1.000 client9 1.000 client10 1.000 client11 1.000 client12 1.000 client13 1.000 client15 1.000 client16 1.000 client17 1.000 client18 1.000 client19 1.000 client20 1.000 client21 1.000 client22 1.000 client23 1.000 client25 1.000 ---- 121 VARIABLE clientUsed.L client is assigned to gridpoint client1 1.000, client3 1.000, client4 1.000, client5 1.000, client6 1.000, client9 1.000 client10 1.000, client11 1.000, client12 1.000, client13 1.000, client15 1.000, client16 1.000 client17 1.000, client18 1.000, client19 1.000, client20 1.000, client21 1.000, client22 1.000 client23 1.000, client25 1.000 ---- 121 VARIABLE numClients.L = 20.000 objective variable ---- 121 PARAMETER DistanceReport distances for assignments gp1803 gp2953 gp6043 gp6811 client1 0.696 client3 0.635 client4 17.994 client5 35.027 client6 36.025 client9 13.202 client10 0.546 client11 14.002 client12 17.751 client13 0.126 client15 16.483 client16 32.622 client17 2.452 client18 6.111 client19 10.735 client20 5.372 client21 16.678 client22 18.504 client23 6.820 client25 16.573 ---- 121 PARAMETER CostAllocation cost = unit cost * distance for assignments gp1803 gp2953 gp6043 gp6811 client1 4.599 client3 3.981 client4 51.078 client5 30.272 client6 36.931 client9 4.162 client10 4.328 client11 10.189 client12 31.181 client13 0.661 client15 29.360 client16 11.137 client17 14.346 client18 37.964 client19 41.799 client20 19.269 client21 40.534 client22 45.598 client23 8.901 client25 62.968 ---- 121 PARAMETER CostReport total cost vs budget total 489.260, limit 500.000
We see again that we can cover 20 clients with the given budget.
Obviously the problem becomes smaller and easier when we use a coarser grid. But we give up some precision. I tried with a \(19 \times 19 = 361\) grid, with coordinates along the axes: \(5,10,\dots,90,95\). This model had only 9,411 binary variables and solved instantaneously. But that optimal solution covered only 19 clients.
Of course this we can also use this approximation model as a way to get a good initial solution. We can solve the MISOCP model using the MIPSTART option. This can shave some time from a solve from scratch.
Conclusion
A variant of the facility location problem seeks to find how many clients we can serve for a given budget. The costs are expressed as \(c_i \times \mathbf{dist}(i,j)\), where \(c_i\) is the unit cost for client \(i\) and \(\mathbf{dist}(i,j)\) is the distance between client \(i\) and facility \(j\).
We formulated the model in four different ways:
- An MINLP model: this is the simplest model, and is a good start formalizing the problem,
- an MIQCP model: this model needs a non-convex solver, so limited progress compared to the previous model,
- an MISOCP model: this model can be solved with convex solvers, but at the expense of being more complex,
- and an approximation: facilities are only to be placed on grid points instead of on any point \((x,y)\). This yielded a large but simple MIP model.
The MIP model was by far the easiest to solve. This is not unusual: quadratic (and second order) models with discrete variables are often much more difficult to solve than linear ones.
References
- Uncapacitated facility location problem with at most k facilities and budget b, https://stackoverflow.com/questions/59327296/uncapacitated-facility-location-problem-with-at-most-k-facilities-and-budget-b
- Solving a facility location problem as an MIQCP, https://yetanothermathprogrammingconsultant.blogspot.com/2018/01/solving-facility-location-problem-as.html
- Jared Erickson, Robert Fourer, Detection and Transformation of Second-Order Cone Programming Problems in a General-Purpose Algebraic Modeling Language, 2019, http://www.optimization-online.org/DB_HTML/2019/05/7194.html
Appendix: MIPstart can disappoint
The MISOCP model is difficult. The final optimal solution of 20 covered clients is found very late in the search. Below is the plot of the MIP bounds: the best found integer solution and the best possible objective:
 |
| MISOCP model: MIP bounds |
We see the integer solutions (light curve) jumping by one (the objective is integer valued). The optimal solution of 20 is not found before spending around 1800 seconds.
In the next picture, I show what happens if we find a solution with an objective of 19 very quickly using the approximation model discussed in the post, and use a MIPstart to feed this into the MISOCP solver. Finding the objective of 19 took us 1000 seconds when solving the MISOCP model from scratch. See the picture above. So this should be a promising approach. Unfortunately, what we see is:
 |
| MISOCP model solved with MIPSTART |
The total solution time has increased. We see the best integer solution of 19 is accepted. However the time to get to 20 is longer than expected. MIPstarts often help in increasing performance, but, as shown here, not always.
