## Friday, December 20, 2019

### A facility location problem

The following problem is from [1]:

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:

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}$

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:

1. 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.
2. 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:

1. $$x,y$$ location of each facility (continuous variables)
2. 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

1. 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
2. Solving a facility location problem as an MIQCP, https://yetanothermathprogrammingconsultant.blogspot.com/2018/01/solving-facility-location-problem-as.html
3. 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.