Social Distancing as p-Dispersion Problem

Received July 9, 2020, accepted July 30, 2020, date of publication August 14, 2020, date of current version August 25, 2020. Digital Object Identifier 10.1109/ACCESS.2020.3016724 Social Distancing as p-Dispersion Problem JAKUB KUDELA Institute of Automation and Computer Science, Brno University of Technology, 601 90 Brno, Czech Republic e-mail: [email protected] This work was supported in part by the Ministry of Education, Youth and Sports of the Czech Republic (Computer Simulations for Effective Low-Emission Energy) under Project CZ.02.1.01/0.0/0.0/16_026/0008392, and in part by IGA under Grant BUT:FSI-S-20-6538. ABSTRACT 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 p-dispersion problem. We use recent algorithmic advancements, that were developed for the p-dispersion problem, and combine them with discretization schemes 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 on the other. INDEX TERMS Social distancing, p-dispersion problem, decremental clustering, COVID-19. 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 pof 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 ppoints 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. II. THE p-DISPERSION PROBLEM A. DEFINITION AND FORMULATIONS In the p-dispersion problem (pDP), we are given a set of npoints, a dissimilarity (or distance) matrix 149402 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ VOLUME 8, 2020 J. Kudela: Social Distancing as p-Dispersion Problem D= {D(i,j):1≤i,j≤n}satisfying D(i,j)≥0 for every 1≤i,j≤nand D(i,i)=0 for every 1 ≤i≤n, and an integer p≥2. The goal is to choose ppoints from the set of npoints 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 pas 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 psuch 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 Dis 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 E(d)= {(i,j)∈E:D(i,j)<d} ⊆ E. 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 D0<D1<· · · <Dkmax be the different non-zero values in D. The associated index sets are K= {1,2,...,kmax}and K0= {0} ∪ K. This formulation uses two types of binary variables: The binary location variable xiindicates if the point i∈Iis selected. For k∈K, the binary variable zkindicates if the location decisions (the particular selection of ppoints) satisfy a minimum distance of at least Dk. The pure binary program is the following: max D0+X k∈K (Dk−Dk−1)zk(1) s.t. X i∈I xi=p(2) zk≤zk−1,k∈K,k>1 (3) xi+xj+zk≤2,k∈K,(i,j)∈E(Dk)\E(Dk−1) (4) xi∈ {0,1},i∈I(5) zk∈ {0,1},k∈K(6) 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 nnodes, denoted by Cis a family {Ci:i=1,...,m}such that Ci∩Cj= ∅ for every 1≤i<j≤mand ∪{Ci:1,...,m} = I. A clustering Cis said to be sufficiently refined if, for every set Ci∈ C,D(Ci):= max{D(u,v):u,v∈Ci,u<v}<z∗, 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 Cbe a sufficiently refined clustering of the nodes of size m. Let DCbe a m×mdissimilarity matrix where DC(i,j)=max{D(u,v):u∈Ci,v∈Cj}. The optimal value ζ∗of the problem pDP(DC,p)provides an upper bound of problem pDP(D,p). 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 Uis computed as the largest dissimilarity between any two points. Using the lower bound L, we build an initial sufficiently refined clustering Cand a reduced dissimilarity matrix DC, using the procedure initialClustering(D,p,L), VOLUME 8, 2020 149403 J. Kudela: Social Distancing as p-Dispersion Problem 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: ki←point in Ciclosest to its center 8: end for 9: d←min(D(ki,kj):1≤i<j≤p) 10: if L<dthen 11: streak ←streak +1 12: else 13: L←d 14: streak ←0 15: end if 16: until streak =s 17: return L Algorithm 2 initialClustering(D,p,L) 1: D,p,L←inputs 2: m←p 3: C←k-means(D,p) 4: compute DC 5: while max(DC(i,i):1≤i≤m)>Ldo 6: i∗←arg max(DC(i,i):1≤i≤m) 7: m←m+1 8: C1 ∗,C2 ∗←k-means(Ci∗,2) 9: Ci∗←C1 ∗ 10: Cm←C2 ∗ 11: recompute DC 12: end while 13: return C,DC 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 Sand W(with S⊆W), where Srepresents the set of optimal nonsingleton clusters (S= {Ci: |Ci| ≥ 2,i=1,...,m}), and Wis the complete optimal solution to the restricted pDP (w.r.t. DC). Iteratively, the sets Sand Ware used to refine the current clustering, resulting in a refined clustering Cand dissimilarity matrix DC, using the procedure splitAndAdd(S,W,C,DC), which is described in Algorithm 3. The resulting reduced pDP is then solved, yielding an upper bound Uon the full problem, and its optimal solution is used to update the sets S,W. The solution procedure solvePDP(DC,p,U,W) has two parts – a 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 Algorithm 3 splitAndAdd(S,W,C,DC) 1: S,W,C,DC←inputs 2: if S= ∅ then 3: C←C,DC←DC(i.e., do nothing) 4: else 5: m←size DC 6: (s∗,w∗)←arg min(DC(s,w):s∈S,w∈W) 7: if w∗∈Sthen 8: i∗←arg max(DC(u,u):u∈ {s∗,w∗}) 9: else 10: i∗←s∗ 11: end if 12: C1 ∗,C2 ∗←k-means(Ci∗,2) 13: Ci∗←C1 ∗ 14: Cm+1←C2 ∗ 15: recompute DC 16: end if 17: return C,DC Algorithm 4 solvePDP(DC,p,U,W) 1: DC,p,U,W←inputs 2: for k=[1 :p+1] do 3: Uh(k)←min(DC(i,j):1≤i,j≤p+1,i6= k,j6= k) 4: end for 5: if max(Uh(k):1≤k≤p+1) =Uthen 6: k∗←arg max(Uh(k):1≤k≤p+1) 7: W← {W(i):1≤i≤p+1,i6= k∗} 8: 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 ppoints out of the p+1 points obtained from the splitAndAdd procedure. If the value of this solution equals the upper bound Ufrom 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 Lof 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 149404 VOLUME 8, 2020 J. Kudela: Social Distancing as p-Dispersion Problem FIGURE 1. Initial mesh (left) and mesh after one round of refinement (right), first example. Algorithm 5 decrementalClustering(D,p,L) 1: D,p,L←inputs 2: L0←heuristicPDP(D,p,10) 3: L←max(L,L0),U←max(D(i,j):1≤i<j≤n) 4: C,DC←initialClustering(D,p,L) 5: S← ∅,W← ∅ 6: repeat 7: C,DC←splitAndAdd(S,W,C,DC) 8: U,W←solvePDP(DC,p,U,W) 9: S← {w∈W: |Cw| ≥ 2} 10: until S= ∅ or L=U 11: return U,X← {Cw:w∈W} FIGURE 2. Initial mesh, second example. 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 pand Dis 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 1and VOLUME 8, 2020 149405 J. Kudela: Social Distancing as p-Dispersion Problem FIGURE 3. Optimal placements, p=5, r=1. Initial mesh with n=102, zr=1,273.88 (left). Final mesh with n=8,745, zr=1,338.09 (right). FIGURE 4. Optimal placements, p=10, r=1. Initial mesh with n=102, zr=749.72 (left). Final mesh with n=8,745, zr=805.25 (right). FIGURE 5. Optimal placements, p=15, r=1. Initial mesh with n=102, zr=552.63 (left). Final mesh with n=8,745, zr=622.75 (right). 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 Twas 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: nthe number of points, 1the square root of the area of the largest triangle in the mesh (a useful measure of the granularity of the grid), rthe rounding factor, kmax the number of distinct elements in D0,z0∗ the optimal objective value of the problem with D0,zrthe ‘‘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 Uis reported in square brackets. If the computation of the instance was terminated prematurely, because during 149406 VOLUME 8, 2020 J. Kudela: Social Distancing as p-Dispersion Problem FIGURE 6. Optimal placements, p=30, r=2. Initial mesh with n=102, zr=325.77 (left). Final mesh with n=1,891, zr=395.42 (right). FIGURE 7. Optimal placements, p=5, r=1. Initial mesh with n=150, zr=1,139.49 (left). Final mesh with n=9,313, zr=1,183.01 (right). FIGURE 8. Optimal placements, p=10, r=1. Initial mesh with n=150, zr=694.62 (left). Final mesh with n=9,313, zr=747.72 (right). the computation, the upper bound Uwas equal to the lower bound Lfrom 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 zr) are shown in Figures 3-6for the first example and in Figures 7-10 for the second one. The first general observations is that in order to decrease 1by half, the number of points nneeds 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 VOLUME 8, 2020 149407 J. Kudela: Social Distancing as p-Dispersion Problem FIGURE 9. Optimal placements, p=15, r=1. Initial mesh with n=150, zr=518.40 (left). Final mesh with n=2,273, zr=555.93 (right). FIGURE 10. Optimal placements, p=30, r=1. Initial mesh with n=150, zr=323.69 (left). Final mesh with n=575, zr=352.96 (right). FIGURE 11. The effect of the mesh refinement on the value of zr. First example on the left, second example on the right. 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 zr. 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 zrand the computational time t, with finer meshes having higher objective value zr, but because of the increase in the number of variables, take progressively 149408 VOLUME 8, 2020 J. Kudela: Social Distancing as p-Dispersion Problem TABLE 1. Results of the computation, first example, max(D)=3,240.8. longer to compute. The improvement of the objective value zrbased 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 palso 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 pincreases. 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 nand unique values kmax 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, VOLUME 8, 2020 149409 J. Kudela: Social Distancing as p-Dispersion Problem TABLE 2. Results of the computation, second example, max(D)=2,180. Algorithm 6 refinePDP(D,p,L,r,T) 1: D,p,L,r,t←inputs 2: D0←round(D,r) 3: U,X←decrementalClustering(D0,p,L) 4: L←U 5: repeat 6: D,p←refineMesh(D,p) 7: D0←round(D,r) 8: U,X←decrementalClustering(D0,p,L) 9: L←U 10: until runTime >T 11: return U,X the optimal placements in Figures 7-10 have a straightforward ‘‘symmetric’’ counterpart with the same objective value. Lastly, the hexagonal pattern, that seems to emerge in Figures 6and 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. 149410 VOLUME 8, 2020