Social Distancing as p-Dispersion Problem

The spread of COVID-19 and similar viruses poses new challenges for our society. There is a strong incentive towards safety measures that help to mitigate the outbreaks. Many countries have imposed social distancing measures that require a minimum distance between people in given places, such as schools, restaurants, shops, etc. This in turn creates complications for these places, as their function is to serve as many people as they were originally designed for. In this article, we pose the problem of using the available space in a given place, such that the social distancing measures are satisfied, as a <inline-formula> <tex-math notation="LaTeX">$p$ </tex-math></inline-formula>-dispersion problem. We use recent algorithmic advancements, that were developed for the <inline-formula> <tex-math notation="LaTeX">$p$ </tex-math></inline-formula>-dispersion problem, and combine them with discretization schemes to find computationally attainable solutions to the <inline-formula> <tex-math notation="LaTeX">$p$ </tex-math></inline-formula>-dispersion problem and investigate the trade-off between the level of discretization and computational efforts on one side, and the value of the optimal solution on the other.


I. INTRODUCTION
The outbreak of the COVID-19 had an enormous impact on the world at large. To mitigate the spread of the virus, various technologies, such as Internet of Things, Unmanned Aerial Vehicles, blockchain, Artificial Intelligence, and 5G are already in use [1]. In this article, we take a look at the problem of positioning people in a given area, such as in a restaurant, school, office, etc., in order to minimize the spread of viruses such as COVID-19. After the initial lockdown, many countries imposed a set of social distancing measures that should help to slow down the spread of the virus. These measures impose a minimum distance between people in a given area. This means that spaces that could previously serve a large number people need to be adjusted for these new measures. As it seems unlikely that we will see the construction of new places that will be designed to abide by these (hopefully temporary) measures, it is only natural to try to find the best use of the ''facilities'' that are already available. However, as the social distancing measures do not have to be stable and can change over time, we will pose the problem of using the available space to its full extent in the following way: Given a fixed number p of people, fit them into a predefined space in such a way, that the minimum distance between any two persons is maximized. Afterwards, by varying p, we can get the optimal (largest) distance that the people can The associate editor coordinating the review of this manuscript and approving it for publication was Wei Liu. be separated by, and, given a particular social distancing rule, we can determine the maximum number (and placement) of people that will fit into the predefined space. The problem of selecting p points in order to maximize the minimum distance between any pair is called the p-dispersion problem [2] and it is one of the classical combinatorial optimization problems. Although easy to formulate, effective and provably optimal methods for solving this problem are quite a recent development. Most notably, the state-of-the-art methods are based on the formulation developed in 2017 in [3] and the most successful method is the decremental clustering scheme published in 2020 in [4]. It is these advancements that made it possible to solve instances of the size sufficient for our purpose. In this article, we devise a discretization scheme that is build on top of the decremental clustering to find computationally attainable solutions to the p-dispersion problem and investigate the trade-off between the level of discretization and computational efforts on one side, and the value of the optimal solution (the minimum distance between any two points) on the other. The investigation is carried out on two numerical examples, the first one is a place with a ''general'' shape, the second one with an ''auditorium-like'' shape.

A. DEFINITION AND FORMULATIONS
In the p-dispersion problem (pDP), we are given a set of n points, a dissimilarity (or distance) matrix D = {D(i, j) : 1 ≤ i, j ≤ n} satisfying D(i, j) ≥ 0 for every 1 ≤ i, j ≤ n and D(i, i) = 0 for every 1 ≤ i ≤ n, and an integer p ≥ 2. The goal is to choose p points from the set of n points in such a way, that the minimum pairwise dissimilarity (the distance between any two points) within the selected points is maximized. The pDP is an NP hard problem [5]. We denote this problem for given input parameters D and p as pDP(D, p). One of the standard applications of the pDP is the location of nuclear power plants, where one is interested in minimizing the risk of losing multiple plants in the event that only one plant is subjected to an enemy attack. To achieve this, the desired selection of plants is that in which the interplant distances are as large as possible [3]. A more peaceful applications of the pDP can be found in location analysis of services, e.g., schools, hospitals, electoral districts, or waste collection plants. A comprehensive survey of the location applications of pDP can be found in [6] and [7]. Another application of the pDP is found in multiobjective optimization -if the Pareto frontier of a problem contains multiple solutions, one can solve a pDP to find p such solutions with most distinct features [8]. In the same paper, an application in portfolio optimization is presented -given a set of potential investment opportunities, one wishes to choose a subset that reduces the closeness in terms of features between the different investment options, which reduces the risk associated with the portfolio.
Within the methodological contributions to the solution of the pDP, several articles have dealt with the problem of solving the pDP to proven optimality. A mixed-integer quadratic formulation was introduced by [9], which can be partially solved by a series of relaxations and reformulationlinearization. A mixed-integer linear formulation of the problem using the ''big M'' constraints was defined in [6]. This formulation can be retroactively thought of as a linearization of the mixed integer quadratic model, that was developed 20 years afterward. Although the linear model is more compact than the quadratic one, it provides much weaker upper bounds.
Our attention will focus on a formulation introduced in [3], which is a novel pure binary compact formulation. Using this formulation, the authors reported substantial computational advancements when compared with the other formulations. Without loss of generality, we can assume that the dissimilarity matrix D is symmetric. Let (I , E) be the complete graph in which points I = {1, . . . , n} are the vertices and E = {(i, j) ∈ I × I : i < j} are the edges. Given any distance d, we define subsets of edges as The compact pure binary formulation exploits the fact that the optimal distance is identical to at least one of the entries in the dissimilarity matrix. Let D 0 < D 1 < · · · < D k max be the different non-zero values in D. The associated index sets are K = {1, 2, . . . , k max } and K 0 = {0} ∪ K . This formulation uses two types of binary variables: The binary location variable x i indicates if the point i ∈ I is selected. For k ∈ K , the binary variable z k indicates if the location decisions (the particular selection of p points) satisfy a minimum distance of at least D k . The pure binary program is the following: The formulation (1)-(6) can be further strengthen using clique-like inequalities and computation can be sped-up by exploiting valid lower and upper bounds [3].

B. DECREMENTAL CLUSTERING METHOD
The decremental clustering method introduced in [4] for the pDP utilizes the formulation (1)- (6). The usage of clustering techniques for finding feasible solutions for combinatorial problems is hardly new. For example, in vehicle routing and scheduling, the ''cluster-first, route-second'' (see [10], [11]) and ''route-first, cluster-second'' (see [12], [13]) paradigms were used to ease the computational burden of the hard combinatorial problem. What sets the decremental clustering method apart is that it provides guarantees for optimality. Decremental clustering was also proposed for the solution of the vertex p-center problem in [14] and [15]. We present the decremental clustering method with the same notation and vocabulary as it was developed in [4]. A clustering of the n nodes, denoted by C is a family where z * is the optimal value of pDP(D, p). The correctness of the decremental clustering method is supported by the following result (proved in [4]).
Lemma 1: Let C be a sufficiently refined clustering of the nodes of size m. Let D C be a m×m dissimilarity matrix where The decremental clustering method works as follows. A lower bound L ≤ z * is computed using a k-means algorithm [16], in a procedure named heuristicPDP(D, p, s) whose pseudocode is described in Algorithm 1. Since the k-means clustering is a stochastic method, it is repeated multiple times as long as the value of the lower bound keeps increasing (the authors of [4] stop after s = 10 iterations without being able to improve its value). An initial upper bound U is computed as the largest dissimilarity between any two points. Using the lower bound L, we build an initial sufficiently refined clustering C and a reduced dissimilarity matrix D C , using the procedure initialClustering(D, p, L), Algorithm 1 heuristicPDP(D, p, s) 1: D, p, s ← inputs 2: L ← 0 3: streak ← 0 4: repeat 5: C ← k-means(D, p) 6: for i = [1 : p] do 7: k i ← point in C i closest to its center 8: end for 9: 10: if L < d then 11: streak ← streak + 1 12: else 13: L ← d 14: streak ← 0 15: recompute D C 12: end while 13: return C, D C whose pseudocode is described in Algorithm 2. The main idea of the procedure is to find the clusters with the largest inter-node distances and split them into two, until the internode distance in all clusters is less than L (which implies that it is less than z * ). After these initial steps, the method uses two auxiliary sets S and W (with S ⊆ W ), where S represents the set of optimal nonsingleton clusters (S = {C i : |C i | ≥ 2, i = 1, . . . , m}), and W is the complete optimal solution to the restricted pDP (w.r.t. D C ). Iteratively, the sets S and W are used to refine the current clustering, resulting in a refined clustering C and dissimilarity matrix D C , using the procedure splitAndAdd(S, W , C, D C ), which is described in Algorithm 3. The resulting reduced pDP is then solved, yielding an upper bound U on the full problem, and its optimal solution is used to update the sets S, W . The solution procedure solvePDP(D C , p, U , W ) has two partsa heuristic ''preprocessing'' method and an exact solver. The heuristic procedure is based on the observation that in a large number of iterations, the optimal value of the pDP problem does not decrease from one iteration to the next, which is 13: U ← U 9: else 10: U , W ← solve (1)-(6) 11: end if 12: return U , W a common feature of decremental relaxation schemes [17]. Therefore, before executing the exact solver, the heuristic procedure checks the best possible selection of p points out of the p + 1 points obtained from the splitAndAdd procedure. If the value of this solution equals the upper bound U from the previous iteration, the associated solution is then optimal, and there is no need to execute the exact solver. In our implementation, the exact solver comprise of solving the model (1)-(6) using the modelling package JuMP [18] in Julia [19], and the GUROBI solver [20]. The pseudocode of the solvePDP procedure is described in Algorithm 4. The whole decremental clustering algorithm is described in Algorithm 5. It also incorporates a possible knowledge on a lower bound L of the optimal value of (1)-(6), which will be explained in the forthcoming section.

III. DISCRETIZATIONS
The continuous variant of the pDP is extremely difficult to solve and the techniques for approaching it are usually bound to convex feasible spaces [21]. In order to apply the pDP framework developed earlier for general spaces, such as classrooms, restaurants, beaches, etc., the feasible space of the problem -the possible locations of the points -is discretized. This discretization is carried out by triangulation (or mesh generation) of the two dimensional feasible area, with the vertices of the triangles being the possible feasible points [22]. Naturally, a question arises about the relationship between the granularity of the triangulation and the objective value of the optimal solution of the associated pDP -the finer the mesh, the better the solution (with higher smallest distance between any two selected points). We address this issue by devising a mesh refinement scheme and solving a series of pDPs on a progressively finer mesh. The mesh refinement works as follows: First, the area of each triangle in the triangulation is computed. The triangles with larger than average area are then split into 4 triangles by adding points is the middle of their sides. The process is illustrated in Figure 1. A side effect of the refinement scheme is that by solving the pDP on a coarser mesh, we obtain a guaranteed lower bound on the optimal value of the pDP on a finer mesh.
Another ''discretization'' can be devised in the space of possible values of the dissimilarity matrix D. As the problem gets more difficult with more unique values in D, with each unique value having a separate binary variable in the model (1)-(6), we can greatly reduce the number of these added variables by rounding the elements of D. The trade-off between the optimal objective value and computational efforts of these two discretization schemes are investigated in the following section. The pseudocode of the method used in the computational experiments, called refinePDP(D, p, L, r, T ), is described in Algorithm 6. It supposes that an initial mesh with p and D is available. The other inputs are a lower bound L (if there is no prior knowledge about a possible lower bound, then L = 0), a rounding factor r (with r = −1 being rounding to the nearest tenths, r = 1 to nearest integers, r = 2 to nearest tens, etc., and r = 0 no rounding) and a time limit for the computation T .

IV. COMPUTATIONAL EXPERIMENTS
We investigate the effect of discretization on two examples. The first one is a general shape depicted in Figure 1 and VOLUME 8, 2020   the second one an auditorium-like shape shown in Figure 2.
In both examples, we examine the computational efforts to solve the pDP problem for progressively finer mesh granularity and for different values of p ∈ {5, 10, 15, 30} and r ∈ {0, 1, 2}. In both examples, the time limit T was set to 24 hours. For the first example, the maximal distance between any two points (''the problem diameter'') was max(D) = 3240.8, for the second example, it was max(D) = 2180. For each problem instance, we report: n the number of points, the square root of the area of the largest triangle in the mesh (a useful measure of the granularity of the grid), r the rounding factor, k max the number of distinct elements in D , z * the optimal objective value of the problem with D , z r the ''real'' objective value (without rounding), and t the time it took to find the optimum. If the computations were not finished (the time limit was reached), the best upper bound U is reported in square brackets. If the computation of the instance was terminated prematurely, because during the computation, the upper bound U was equal to the lower bound L from the solution of coarser discretization (i.e., the finer discretization did not improve on the optimal value of the solution) the instance is marked with a '*'. The instances that were not computed, because the time limit was already reached, are marked by a '-'. The computations were carried out on an ordinary computer with 3.2 GHz i5-4460 CPU and 16 GB RAM.
The numerical results of the computations are reported in Table 1 for the first example and Table 2 for the second one. The optimal placements (the ones with the best value of z r ) are shown in Figures 3-6 for the first example and in Figures 7-10 for the second one. The first general observations is that in order to decrease by half, the number of points n needs to be roughly quadrupled (which follows from the way the mesh gets refined). The second general observation can be made about the impact of the rounding factor r: Apart from a single instance, there was no difference in the ''real'' objective value between instances with r = 0 and r = 1, while there was a huge difference   between the computational time t. From these computational experiments, there is no doubt that the rounding procedure presents a substantial benefit, as the instances with rounding were computed around an order of magnitude faster than the instances without rounding, and some large instances could not be computed within the time limit without the use of rounding.
The difference between r = 1 and r = 2 is much more nuanced. In 83 % of the instances, the computations with r = 2 were faster. On the other hand, using r = 2 instead of r = 1 results on average in 0.43 % worse value of z r . Premature termination is also more prevalent in the r = 2 case. Of the 48 successfully computed instances it occurred 14 times for r = 2, compared to 8 times for r = 1.
The coarseness of the mesh naturally plays a crucial role in both the objective value z r and the computational time t, with finer meshes having higher objective value z r , but because of the increase in the number of variables, take progressively longer to compute. The improvement of the objective value z r based on the mesh refinement is captured in Figure 11, which shows the percentage improvements caused by the increases in the number of mesh points n. Additionally, the value of p also has an extensive impact on the computational time. Similar to the findings in [4], we also find that the decremental clustering scheme works very well for smaller values of p, but the computations become progressively more costly as p increases. The value of p = 30 seems to be close to the limit of applicability of the method. On the one hand, it means that in the context of social distancing it is well applicable for use in situations, when the available space does not allow for more than 30 persons, such as in classrooms, restaurants, or offices. On the other hand, it still can be used to compute valid upper bounds on the objective value even for larger problems, which can be explored by various heuristics.
There is also a significant difference in the ''difficulty'' between the two examples. Although the number of points n and unique values k max were similar for both examples in the individual mesh refinement steps, the computational times differ quite a lot, mainly for larger values of p. This can be attributed to the ''dual degeneracy'' [4] occurring when a larger number of clusters can be rearranged from one iteration to the next to find solutions of the same cost. The second example is more symmetric than the first one, meaning that it has a larger number of the possible rearrangements. Naturally, Lastly, the hexagonal pattern, that seems to emerge in Figures 6 and 10 (both with p = 30) is no accident -the hexagonal pattern is optimal for many location problems (including the p-dispersion problem) with numerous facilities covering a large area [23].

V. CONCLUSION
In this article we have studied the problem of locating persons in a given area, that should abide to social distancing measures such as those arising in the time of COVID-19 and similar viruses. We have argued that the p−dispersion problem can be used to efficiently model these situations. We devised a discretization scheme that was build on top of the decremental clustering method to get computationally attainable solutions, which worked very well in the computational study on the two artificial examples, especially for smaller values of p.
We have investigated the effect of rounding of the dissimilarity matrix D on the computational effort and conclude that it is an indispensable part of the discretization scheme that has virtually no disadvantages in terms of the quality of the obtained solution. We have also seen the substantial increase in computational efforts for higher values of p, which can be contributed to the ''dual degeneracy'' of the clustering scheme.
For future research, fast heuristics that run parallel to the decremental clustering scheme might further improve on the computational time and the size of the problems that are solvable by the presented method. Also, a mesh refinement scheme based on the current optimal placement (instead of the triangle size as presented here) could lead to improvements in the objective value.