extfieldsofstudy:
Mathematics
Computer Science

text:
The intersection of location-allocation and clustering

Location-allocation and partitional spatial clustering both deal with spatial data, seemingly from different viewpoints. Partitional clustering analyses datasets by partitioning them into separate groups, while location-allocation places facilities in locations that best meet the needs of demand points. However, both partitional clustering and location-allocation can be formulated as optimization problems minimizing the distances of (demand) points from their associated centers (facilities). Further, both techniques consider certain extensions such as weighted data points, different distance metrics, capacity constraints, different membership types, outliers, and selecting the number of clusters or facilities. In this article, we highlight and review the similarities and differences of these techniques, compare them with model-placed clustering, and provide a unified theoretical framework covering both approaches. We look at a number of extensions common for both approaches, propose a new spatial clustering method, PACK, which combines and adjusts those extensions, and provide software tools (rpack) for conducting spatial analysis with the new method.

In this work, we study partitional (also referred to as non-hierarchical) clustering. Partitional clustering aims to partition the data set into nonoverlapping subsets such that each data point is in exactly one subset and each cluster can be represented by one point, referred to as a prototype (Jin and Han, 2010;Xiao and Yu, 2012). Partitional clustering is especially useful in applications which use the partitions for further analyses. As an example, direct marketing campaigns often start with segmenting customers into groups (see further examples in Banerjee and Ghosh, 2006).
Further, clustering often involves a spatial dimension with geographic information related to the studied phenomenon. Such a setting, referred to as spatial clustering, requires an appropriate and meaningful treatment of space, spatial relationships, and the attributes of locations (Grubesic et al., 2014). We focus in particular on partitional spatial clustering, where points of interest are partitioned into disjoint clusters. Partitional spatial clustering is used, for example, for the spatial analysis of Internet of Things (IoT) sensor data. IoT is a key emerging technology where the spatial dimension is often important (Lee and Lee, 2015). IoT sensors observe their local environment, creating massive amounts of data to be analyzed and acted upon. Partitioning the sensor data into local clusters can help in finding the local features of the observed phenomena, but also in distributing the computational burden of the data analytics, especially if local or edge-based computing capacity is available (Yi et al., 2015;Xu et al., 2017;Lovén et al., 2019).
Location-allocation, a concept of operations research, aims to solve a logistic optimization problem of locating facilities (e.g., hospitals) to best meet the needs of demand points (Hale and Moberg, 2003;Revelle et al., 2008). Location-allocation was first studied in the sixties (Cooper, 1963), but it has recently regained considerable interest. Faced with increasing amounts of spatial data, business intelligence and business analytics are turning into location intelligence and location analytics (Panian, 2012;Garber, 2013;Open Geospatial Consortium, 2012). The most common types of location-analytic tasks include finding the best places to locate facilities, for example stores and warehouses (Garber, 2013), given geospatial information on, for example, local population densities or their mobility patterns. Such tasks can essentially be interpreted as location-allocation problems.
Albeit the aims of location-allocation and partitional clustering may sound different, there is a strong connection between them. Both approaches can be formulated as optimization problems, minimizing the distances from the data points to some center points, be they cluster centers or facilities. Both can also be connected to model-based clustering, where the data points of each cluster are assumed to emerge from a parametric distribution family (see e.g. Raykov et al., 2016;Bishop, 2006). Further, both approaches have similar extensions. These include different types of distance metrics, restrictions on the locations of the centers, weighted data points, controls on the number of data points assigned to a given cen-ter, varying membership types, outliers, and ways to determine the number of facilities or clusters. Both approaches implement these extensions similarly (e.g., by including an optimization constraint), despite the differences in interpretation.
There is a vast number of articles proposing different approaches for conducting location-allocation or clustering analyses with these extensions. However, typically only one or two extensions are considered at a time. Still, in many applications, it is necessary to be able to consider several extensions simultaneously. Further, currently, there is no open source or commercial software readily available for location-allocation and clustering with all the aforementioned extensions and constraints.
The aim of this article is twofold. First, we bring partitional clustering and location-allocation, two seemingly very different approaches, closer together by reviewing and highlighting their similarities and building a unified theoretical framework which encompasses both. Further, for getting insights on the implicit assumptions the approaches make, we draw parallels between them and model-based clustering. Second, we present a novel partitioning approach that facilitates all the above mentioned extensions, and provide an extensive, open-source toolkit for conducting the proposed analyses 1 .
We refer to the proposed approach as PACK (PlAcement with Capacitated K-family 2 ). We have used earlier versions of PACK, with fewer features, in our previous articles for placing edge servers Lovén et al., 2020) and maternity hospitals (Huotari et al., 2020). However, compared with the earlier versions, we make several important enhancements. First, similarly as Lovén et al. (2020), we consider a scenario where some centers are fixed, i.e. have a known location. However, such fixed centers may result in a sub-optimal partition, and relocating some of the fixed centers may produce better results. In the location-allocation context, such relocation causes expenses. Thus, we extend the approach such that the fixed centers can be relocated if a given penalty is paid.
The revised PACK identifies outliers as well, that is, isolated points in the data that could distort the clustering structure or be considered too expensive to be offered the service of facility in location-allocation. Accordingly, we allow a point not to be a member of any center, if it appears to be an outlier. PACK also considers the selection of the number of centers.
Furthermore, in some clustering applications there may be a need to incorporate prior information about the locations of the cluster heads. As an example, one may have access to the clustering results of earlier data sets which could then be used to enhance the clustering of a new data set. However, there appears to be lack of methods that directly address this problem.
PACK enables the incorporation of such prior information as follows. In Lähderanta et al. (2020), some of the demand locations were assigned with a location priority, indicating they should have a center close to them. We show that in the clustering framework, this leads to a semi-Bayesian procedure where prior information related to the locations of the centers is incorporated via pseudo-data.
We thoroughly examine the effects of different features and their combinations using simulated data as well as two real world examples. The first real world example is related to the placement of edge computing servers in Shanghai, a topic currently under active research and the other is related to location intelligence, placing recycling centers in a town in Slovakia.
The rest of the article is organized as follows. In Section 2 we first overview the theory of partitional clustering and location-allocation, and their connection to model-based clustering, and further discuss extensions related to the locations of the centers (3.1), distance metrics (3.2), weighted data points (3.3), capacity (3.4), membership (3.5), outliers (3.6) and the number of centers (3.7). In Subsection 4.1 we describe the PACK method that considers all the extensions simultaneously, and in Subsection 4.2 we consider how the features of the proposed method can be interpreted as prior information on the location of the centers. In Subsection 4.3 we describe the proposed algorithm for solving the problem and discuss how the different extensions affect the algorithm. Finally, we demonstrate the effect of the extensions with both artificial and real data in Section 5 and discuss our findings in Section 6.

Framework
We start by formulating both location-allocation and partitional clustering as optimization problems, and then continue by considering their connection to model-based clustering.

Partitional clustering and location-allocation
Given a data set of n points, let x i = [x i1 , x i2 ] be the spatial coordinates of a point i. We denote by c j = [c j1 , c j2 ], j = 1, . . . , k the coordinates for the k facilities or cluster heads. For unifying the terminology, we will refer to c j as centers. The distance between a point i and a center j is denoted by d(x i , c j ). In addition, let us further denote by y ij the membership of a point i to the center j, where y ij equals one if the point i belongs to the cluster j and zero otherwise.
In both frameworks, the aim is to minimize the sum of the distances between points and the centers they are assigned to. Thus, we consider following objective function with the following two constraints: The constraint (2) dictates that a point i is either assigned to cluster j or not and constraint (3) ensures that each point is assigned to exactly one cluster.

Model-based clustering
We consider two approaches to model-based clustering, the mixture model approach and the classification maximum likelihood approach.

Mixture-models
In a mixture model (see e.g. McLachlan and Peel, 2000), data is assumed to have emerged from a mixture of parametric density family where f (x|θ j ) is the density function of a data point x given that it has emerged from cluster j and p(π j ) is the probability of observing a data point from cluster j. The number of clusters k is predetermined. This corresponds to a data generating process where for each x i , i = 1, . . . , n. Given a data set x 1 , . . . , x n , the aim is to estimate the parameter vectors θ = [θ 1 , . . . , θ k ] and π = [π 1 , . . . , π k ] by maximizing the logarithm of the likelihood function l(θ, π; x 1 , . . . , Each point x i is then assigned with a probability of belonging to each cluster j. The data set x 1 , . . . , x n can then be clustered by assigning each x i to the cluster with the highest likelihood. See Fraley and Raftery (2002) for further information.
Maximizing the log-likelihood of the mixture model (5) can be connected to the minimization problem (1) when specific assumptions are made regarding the distance d, the parameters θ j and proportions π j of the mixture densities.
For example, if the distance d is measured with a squared Euclidean distance metric, it can be shown that the minimization of (1) corresponds to estimating the parameters of an isotropic (i.e., spherical) Gaussian mixture model where each cluster is assumed to have the same proportion π j and the cluster variances tend to zero (see e.g. Bishop, 2006). This result is known as the small variance asymptotic derivation of K-means clustering (Jiang et al., 2012).

Classification maximum likelihood approach
Let us denote by P = [P 1 , . . . , P k ] the partition of the n points into k groups where P j contains the points assigned to group j. These groups correspond to the points assigned to a given facility in location-allocation and points assigned to one cluster in partitional clustering.
In the classification maximum likelihood approach (Celeux and Govaert, 1992), this unknown partition P is a parameter vector and the log-likelihood function to be maximized is Celeux and Govaert (1992) use also a more flexible version obtained by incorporating the proportions of clusters π j as follows l(P, θ, π; x 1 , . . . , where n j is the number of points in the group j. This corresponds to (6) if the proportions π j are assumed to be equal for each j = 1, . . . , n. Thus (6) implicitly assumes that the true clusters have equal sizes.
We continue to follow Celeux and Govaert (1992) and explore the connection between classification maximum likelihood and partitional clustering and location-allocation. Let us assume that f (·|θ j ) ∼ N (µ j , σ 2 I) for all j. Then the classification likelihood is l(P, µ, σ; x 1 , . . . , where a is the dimension of data (here a = 2), and c is a constant independent of P, µ and σ 2 .
If σ 2 is known, With unknown σ 2 , it is straightforward to show that (7) and by plugging it into (7) l(µ, P |x 1 , . . . , Hence, both with known and unknown σ 2 the classification maximum likelihood with Gaussian distribution corresponds to minimizing (1) with the squared Euclidean distance metric.
On the other hand, if for each cluster, x i1 and x i2 are independently generated from univariate Laplace distribution, i.e. f (·|θ j ) ∼ Laplace(µ j , σI) for all j, the classification likelihood function is With an analogous reasoning as with the squared Euclidean distance metric, maximizing (9) is equivalent to minimizing the sum of Manhattan distances k j=1 xi∈Pj Classification maximum likelihood thus corresponds to solving the clustering problem (1) with the Manhattan distance metric.
Similarly, minimizing (1) with the Euclidean distance metric corresponds to maximizing the classification maximum likelihood, assuming that the data of each cluster is generated from an isotropic multivariate Laplace distribution defined as a special case of the multivariate exponential power distribution (Naik and Plungpongpun, 2006;Arslan, 2010) as follows:

Extensions
In the following subsections, we consider certain extensions common to locationallocation and partitional clustering and highlight their connections to modelbased clustering. Such extensions include constraints on the locations of the centers, different distance metrics, weighted data points, capacity constraints on the clusters, different membership types, outliers, and selection of the number of centers.
3.1 Constraints on the possible locations of the centers So far, we have assumed that the centers can be located anywhere in the region of interest. In location-allocation, problems where the facilities can be placed anywhere in the area are referred to as continuous location-allocation problems or site generating problems (Brimberg et al., 2008). In clustering, algorithms with unrestricted placement of centers are referred to as virtual point prototype algorithms (Xiao and Yu, 2012). However, in both location-allocation and clustering, there are applications where the centers can be placed only on a predetermined set of points {p 1 , p 2 , ..., p m }. This can be achieved by including a constraint c j ∈ {p 1 , p 2 , ..., p m } ∀j to the optimization problem (1).
In location-allocation, this is a common constraint and the modified problem is referred to as the discrete location-allocation problem or the site selection problem. In clustering, there is less often a need to constrain the locations of the cluster centers, but if such a constraint is applied, the cluster centers are typically constrained to belong to the data points. Thus, the constrained algorithms are referred to as actual point-prototype-based clustering algorithms (Xiao and Yu, 2012).
In the following sections we will refer to the setup with constrained locations of centers as the discrete setting, and unconstrained locations as the continuous setting.
A third option for the cluster centers is that the locations of some centers may be fixed. In location-allocation, this is a typical setting, with new facilities being added to an existing network (see e.g. Rahman and Smith, 2000). For clustering, such a variant may correspond e.g. to a setting where, based on prior information, it is known that there are clusters with centers in specific locations.

Distance metrics
Section 2.2 described the connection between model-based clustering and the optimization problem (1) while applying the squared Euclidean, Manhattan and Euclidean distance metrics. However, the obtained methods are referred to differently depending on the framework and whether the continuous or discrete setting is applied. Here, we summarize the distance metrics along with the names of the algorithms in both settings and with different frameworks. We also consider other typical distance metrics considered in the frameworks. The metrics, along with the names of the methods, are listed in Table 1.
Squared Euclidean distance. As demonstrated in Section 2.2.2, minimizing the loss function (1) in the continuous setting with the squared Euclidean distance metric corresponds to maximizing the classification maximum likelihood with an isotropic Gaussian distributions with equal variance and equal proportions.
For the continuous setting in partitional clustering , the squared Euclidean distance metric corresponds to the famous K-means clustering (see a recent review by Zhao et al., 2018), whereas in the discrete setting, the metric corresponds to the K-medoid clustering (Kaufman and Rousseeuw, 1987).
The squared Euclidean distance metric has no physical interpretation and is not representative of travel, transportation or movement distance (Murray and Estivill-Castro, 1998). It is thus less popular in location-allocation. However, in some special occasions, travel costs increase quadratically with the Euclidean distance. As an example, Zhou and Li (2013) studied placing emergency service facilities, where the costs caused by fires can be expected to be quadratic with respect to the distance from the emergency center. Further, in the continuous setting, squared Euclidean distances may be used as approximations to Euclidean distances, as faster algorithms can be applied (Farahani and Hekmatfar, 2009).
Euclidean distance. Minimizing Euclidean distances in (1), on the other hand, corresponds to maximizing the classification maximum likelihood, assuming that the points have emerged from an isotropic Laplace distribution (2.2.2). In partitional clustering, the Euclidean distance metric is more rarely used, and the approaches are referred to as Euclidean K-median (see e.g. Kolliopoulos and Rao, 1999) or simply K-median (see e.g. Har-Peled and Mazumdar, 2004). In the discrete setting, this approach, too, is referred to as K-medoid clustering (Kaufman and Rousseeuw, 1987).
In location-allocation, the Euclidean distance metric is often used, as the aim is to directly minimize the physical distances between the demand points and the assigned facilities. The discrete setting is referred to as the p-median method (see reviews by Mladenović et al., 2007;Daskin and Maass, 2015) and the continuous setting as the multifacility Weber problem (see the review by Brimberg et al., 2008).
Manhattan distance. In model-based clustering, the Manhattan distance metric corresponds to assuming that the coordinates of each data point are generated from a univariate Laplace distribution (see also Cord et al., 2006). With the Manhattan distance metric, the method is referred to as K-median clustering in the continuous setting and, in the discrete setting, again as Kmedoid. In location-allocation this distance, also known as the city-block dis-tance, may be used as a predictor of physical distance for example in an urban environment (Brimberg et al., 2008).
Mahalanobis distance. The squared Euclidean distance metric implies an underlying assumption that the clusters are spherical, with similar shapes and sizes. With Euclidean or Manhattan distance metrics, remoteness from centers is not penalized as strongly. Still, if the aim is to detect the natural groupings in the data, both distances may produce sub-optimal clustering if the natural cluster sizes in the data are heterogeneous with, for example, elongated shapes. A generalization that alleviates this problem is to apply the squared Mahalanobis distance metric, which allows elongated shapes and different spatial sizes for each cluster (Cerioli, 2005;Sung and Poggio, 1998). This approach corresponds to assuming that the data of each cluster has emerged from an anisotropic multivariate Gaussian distribution. Sung and Poggio (1998) referred to the approach as Elliptical K-means.
The challenge with the Mahalanobis distance metric is the difficulty in the initial estimation of the related covariance matrices. Several approaches have been proposed (Xiang et al., 2008;Melnykov and Melnykov, 2014;Chokniwal and Singh, 2016). Further, Lapidot (2018) recently highlighted convergence problems with K-means algorithms based on the Mahalanobis distance metric, and included regularization parameters for each cluster to achieve convergence. The first proposed regularization penalizes large eigenvalues in the covariance matrices and the second one constrains the volumes of each cluster as constant. The Mahalanobis distance metric has also been combined with the K-median method by Ackermann et al. (2010), who specify a generalized K-median problem in the discrete setting, where an arbitrary dissimilarity measure is minimized and the Mahalanobis distance metric is used as an example.
Threshold distance. In location-allocation, typically Mahalanobis distance type generalization is not needed, as the aim is to minimize the travel distances regardless of the underlying structure of the data. However, different distance metrics have been proposed also in the location-allocation context. A loss function based on the Euclidean distance metric can be considered to minimize the average travel distance from the centers to the demand points. However, sometimes the aim may, instead, be to offer a facility within a given distance to as many demand points as possible. This corresponds to minimizing the percentage of demand where the distance to the nearest facility exceeds a given threshold. In location-allocation literature this is referred to as the maximal coverage location problem (MCLP) (Church and ReVelle, 1976;Church and Weaver, 1986). Church and Weaver (1986) have shown that the MCLP can be connected to the p-median problem by replacing the Euclidean distance metric with a thresholded distance metric where the distance is 0 if the Euclidean distance is under the given threshold and 1 if the threshold distance is exceeded.
-Discr p-median -Distance matrix based distance. Further, in location-allocation, Euclidean or Manhattan distance metrics are often only approximations of the actual travel distances or costs. In the discrete setting, the actual travel distances, or even the travel costs incurred, can be used in the objective function by pre-computing all the pairwise distances between the demand points and the facility candidates. Such approaches using pre-computed distances are also referred to as p-median problems (Mladenović et al., 2007), as the term p-median is not tied to the Euclidean distance metric.

Weights
For location-allocation, each point i is typically assigned a weight w i that corresponds to the amount of demand at that point such as the number of people living in an area. Weighted data points are rarer in clustering, but may similarly correspond to, e.g., the duplication of data (Borgwardt et al., 2017). For such weighted optimization problems, the objective function (1) Ackerman et al. (2012) analyze clustering methods with respect to how they respond to the weighted data points, and summarize that K-means, K-median and K-medoid are weight separable. In other words, by modifying the weights of points, these algorithms separate any points to different clusters. Indeed, the effect of the distortion of a point i in the objective function (1) is proportional to the weight w i . Thus, such weight-sensitive techniques have two potentially conflicting considerations: the weight of the points and the geometry of the data. The points with heavy weights tend to attract the center even if they appear to be isolated.
Due to the tendency of heavy points to attract cluster centers, the weight can be used also for indicating points that are more important when locating a center. Wheeler (2019), for example, create optimal police patrol areas by considering the number of calls weighted by the priority or duration of the call. Similarly, in case of clustering, Tseng (2007) interpret weights as an importance measure or, more generally, a prior measure which preferred or prohibited patterns of cluster selections.
Also in model-based clustering, weights can be interpreted as measures of importance. Gebru et al. (2016), for example, consider Gaussian mixture models with weighted data. They interpret the weight of a point i as its relevance, that is, the higher the weight of i, the stronger its impact on the clustering. They conclude that w i can be considered as point-wise precision. In more detail, given a cluster, the points with a heavy weight are assumed to be drawn from a distribution which has a smaller variance than points with a small weight and are thus presumably located closer to the cluster center.
Above, weight is considered as an attribute of a data point. In clustering, weights can more generally be assigned to the variables, i.e. the dimensions of the data points (see a recent survey by de Amorim, 2016). Thus, some variables are considered more important than others in clustering. Typically, such weights are estimated during the clustering process. As our focus is on spatial data, the weights are only assigned to the points.

Capacities of facilities and sizes of clusters
Clustering is often used for dissecting data into separate groupings for further processing or analysis. Such cases often aim at balanced clustering, where each cluster has a similar number of data points in each cluster, and too large or too small clusters are undesirable. Such balance may be desirable even if the "natural" clusters in the data were imbalanced (Banerjee and Ghosh, 2006). Section 2.2.2 showed that the minimization of (1) with the Euclidean, squared Euclidean and Manhattan distance metrics implies an assumption that the true clusters have equal proportions. This assumption was demonstrated also by Xiong et al. (2008) who formally illustrate that K-means tends to produce nearly equal-sized clusters even if the data has varying "true" cluster sizes. K-means and K-medians thus have a built-in tendency to produce balanced cluster sizes. However, neither guarantees balanced cluster sizes if the true sizes vary severely. Clustering methods that specifically aim at balancing the partitions are thus needed.
Clustering methods that aim at balance can be partitioned into two groups, balance constrained and balance driven approaches (Malinen and Fränti, 2014). In balance constrained clustering, the cluster size balance must be met, and the minimization of distances is only a secondary criterion. Balance driven clustering, on the other hand, is a compromise between balance and distance.
Several approaches have been proposed for both approaches. Balance constrained include Banerjee and Ghosh (2006), Malinen and Fränti (2014) and Elliott (2011). Other constrained approaches include Hu et al. (2018), Zhu et al. (2010) and Ganganath et al. (2014), where Hu et al. (2018) constrain the size of the clusters only from below, and Zhu et al. (2010) and Ganganath et al. (2014) allow predefined size constraints to vary between clusters. Li et al. (2018) and Liu et al. (2017) studied balance-driven clustering with a soft constraint that penalizes large deviations in cluster sizes. See also a recent survey on balanced data clustering algorithms by Gupta (2017).
In model-based clustering, the constraints in the cluster size can be set as prior distributions. Jitta and Klami (2018) proposed flexible priors where each possible cluster size can be assigned with a prior probability. On the other hand, Levine et al. (2012) used a Poisson distribution as a prior distribution for the number of points in a cluster.
In location-allocation, the facilities may have a restricted amount of supplies, or capacity, that can not be surpassed. Thus the facilities must be placed and the demand allocated such that these capacity limits are followed. Such placement problems are referred to as capacitated location-allocation problems and they include a hard upper capacity constraint, in the optimization problem (10).
In his seminal work, Cooper (1972) applied capacity constraints in the continuous setting, focusing on the Euclidean distance metric. Later, capacity constraints have been used with various distances both in the discrete (e.g. capacitated p-median problem, CPMP) and the continuous setting (e.g. capacitated multifacility Weber problem, CMWP).
Note that in (11) the constants a i may not be equal to the weights w i . For example, Negreiros and Palhano (2006) assign no weights in the objective function but use the demand only in the constraint to satisfy the capacity constraints. Further, also heterogeneous capacity constraints have been considered, as some applications have facilities with varying capacities (Baranwal and Salapaka, 2017;Liao and Guo, 2008).
Hard capacity constraints have also been proposed in the clustering framework by Mulvey and Beck (1984), who assumed that the weights w i = 1 and coefficients a i = 1 for all i and referred to the problem as the capacitated clustering problem (CCP). Later, Negreiros and Palhano (2006) used squared Euclidean distances in a continuous setting and referred to the problem as capacitated centered clustering problem (CCCP).
Location-allocation may also aim for balance, as facilities with underused capacity may be undesirable. For example, Negreiros and Palhano (2006) combine the aims of obeying the capacity constraints and balancing the data by setting a constraint that controls the number of points in each cluster, while Aras et al. (2007) assume that the sum of demand assigned to a facility must be equal to its capacity.
Another possibility for combining balance and capacity constraints is to apply both a lower and an upper limit for the used capacity (Borgwardt et al., 2017) by including a constraint In model-based clustering, such a hard capacity interval can be interpreted as a uniform prior distribution for used capacity (Klami and Jitta, 2016).

Membership
Four different membership types are considered in both location-allocation and partitional clustering, namely, hard, uncertain, overlapping and fractional.
Hard membership. The objective function (1) forces each point to be assigned to exactly one cluster. In clustering, this is referred to as hard or crisp clustering, and the membership of a point i to center j as hard membership or hard assignment. In location-allocation literature, this corresponds to an underlying assumption that each demand point is served by exactly one facility and thus referred to as single-sourcing.
Uncertain membership. Uncertain membership means that a point can belong to many clusters in some degree. In model-based clustering with the mixture model (4), each point is assigned the probabilities of belonging to each cluster. Thus, the memberships are of an uncertain type per se. In partitional clustering and location-allocation, an uncertain membership typically corresponds to fuzzy logic (Klawonn and Höppner, 2003;Canós et al., 2001), where each point is assigned to each cluster with a membership degree that varies between 0 and 1.
A fuzzy clustering version of the K-means is referred to as fuzzy c-means (FCM) (Dunn, 1973;Bezdek, 1981). The objective function is The parameter 1 ≤ m < ∞ determines the level of fuzziness in the clustering, the larger the m the fuzzier the result. While a typical choice is m = 2, limit m = 1 corresponds to hard clustering, and as m → ∞, the points have equal memberships for all clusters and the cluster centers converge to the center of the data (Klawonn and Höppner, 2003).
Both model-based and FCM (Bezdek, 1981) apply a constraint requiring the memberships of a point i to sum up to 1. Krishnapuram and Keller (1993) omit this constraint, and the resulting membership values may be interpreted as degrees of possibility of the points belonging to the clusters, i.e., the compatibilities of the points with the prototypes. Ben-Israel and  propose yet another alternative, referred to as probabilistic d-clustering, that is based on an assumption that the probability of a point belonging to a cluster is inversely proportional to its distance from the cluster center. D-clustering was subsequently extended for unequal cluster sizes using the Mahalanobis distance metric .
Overlapping membership. Overlapping membership means a point can be assigned to several centers. In their review article, N'Cir et al. (2015) divide partitional clustering approaches with overlapping membership into two groups: those with uncertain overlapping memberships, and those with hard overlapping memberships.
Fuzzy c-means can be extended to overlapping clustering by fixing a threshold such that all observations whose degrees of memberships exceed the threshold are assigned to the respective clusters. Similarly, in model-based clustering, overlapping membership could be implemented by assigning a threshold for the membership probability. However, there is a fundamental problem with such straightforward modifications. Since both algorithms work in two phases, they initially aim at optimizing an objective function under constraints that do not match with the expected overlapping structure of the clustering. Further, there is a problem of choosing the suitable threshold (Cleuziou, 2008;Banerjee et al., 2005).
There are approaches that specifically combine uncertain membership clustering with the intuition of having overlapping clusters (see a review of N'Cir et al., 2015). As an example, in Evidential c-means (Masson and Denoeux, 2008) all the possible combinations of clusters are evaluated and a mass of belief is allocated within each possible combination. In the model-based framework, Banerjee et al. (2005) combine overlapping and uncertain memberships using a data generating function that omits the multinomial constraint.
There are several approaches for hard overlapping membership. For instance, Cleuziou (2008) propose an overlapping K-means algorithm where each point is assigned to at least one cluster, and the aim is to minimize the squared distances between the points and the averages of the centers they are assigned to. Other partitional clustering approaches with overlapping membership are reviewed by Baadel et al. (2016) andN'Cir et al. (2015).
In the location-allocation context, the view point for overlapping membership is different. Instead of identifying the points that should be assigned to several centers, either all or a predetermined subset of the points are assigned to several facilities, one of which is considered as a primary facility and the others as back-up facilities. This results, for example, in a more reliable system, where the assigned back-up facilities can be utilized in case the primary facility is not in use or is engaged with other demands. If back-up facilities are assigned for all points, such models are referred to as Q-coverage (or K-coverage) location models, and the constraint (3) where Q is the number of assigned facilities (Karatas et al., 2016(Karatas et al., , 2017. Fractional membership. In fractional (also known as partial) clustering, a point can be shared by several clusters (Borgwardt et al., 2017). Thus, the membership of a point i to cluster j is given by replacing the constraint (2) with a constraint y ij ∈ [0, 1] ∀i, j.
In location-allocation, fractional membership is referred to as multi-sourcing, and interpreted such that the demand of a point can be catered by several facilities. If no capacity constraints are applied, each point is assigned to its closest center and thus, like with fuzzy c-means with m = 1, fractional membership corresponds to hard membership. However, if capacity limits are applied, fractional membership allows a more flexible allocation of points into centers than the capacitated hard membership. In the extreme case, with hard membership, it may even be that the data points can not be partitioned according to given capacity limits (Borgwardt et al., 2017).
In model-based clustering, Heller et al. (2008) utilize fractional membership in the Bayesian framework, assuming that the memberships of points follow a Dirichlet distribution instead of the multinomial distribution as in the mixture models.
Thus, with both uncertain and fractional memberships, a point can belong to several clusters with different degrees, with the degrees summing up to one. Still, these membership types are different on a conceptual level. Indeed, uncertainty on which center a point is assigned to is different from the certainty of a point belonging partially to several clusters (Heller et al., 2008).
The conceptual difference between overlapping and fractional clustering becomes evident with capacity constraints. In fractional clustering, fractions of the weight of a point are assigned to clusters, whereas in overlapping clustering the whole weight would be assigned to several points, thus requiring more capacity in the system.

Outliers
With the aforementioned membership types, each point is assigned to a cluster, albeit the assignment may be uncertain, overlapping or fractional. However, if a point does not seem to belong to any cluster, a forced assignment may distort information and the interpretation of the cluster it is assigned to (Tseng, 2007). Similarly, in location-allocation, it may be economically essential to be allowed to exclude remote locations from service (Charikar et al., 2001). Thus, in addition to allowing a point to be a member of several centers, it may be necessary to allow a point not to be a member of any center. In both locationallocation and partitional clustering, such points that are left unassigned are referred to as outliers. Further, in partitional clustering, a partition where every point is assigned is referred to as complete clustering, and a partition where some points may be left unassigned is referred to as partial clustering (Steinbach et al., 2019).
The sensitiveness of the partition to the outliers depends on the method applied. For example, point prototype based clustering methods, where cluster centers are restricted to the data points, are more robust to outliers than the virtual point prototype based algorithms where the centers can be located anywhere (Xiao and Yu, 2012). Further, the K-median method is known to be less sensitive to outliers than the K-means method (García-Escudero et al., 2010), and can thus be considered as a robust alternative to K-means (Raykov et al., 2016).
Unlike the squared Euclidean distance metric, an Euclidean distance metric may mitigate the effect of outliers. Still, if outlier points are sufficiently far away from cluster centers, the outliers do distort the result even with the Euclidean distance metric (cf. García-Escudero et al., 2010). Thus, an explicit treatment of outliers is needed.
Outliers could be removed prior to the clustering by using some outlier detection method. However, the clustering structure may be highly dependent on the outliers, and on the other hand, the determination of the outliers depends on the clustering structure. Identifying the outliers should thus not be a separate step, but a part of the clustering process (Liu et al., 2019).
Outliers can be identified by determining the number of outliers prior to analysis, and iteratively assigning the most remote points as outliers. Such an approach has been proposed both in location-allocation (Charikar et al., 2001) and clustering (Chawla and Gionis, 2013). Whang et al. (2015) combine outlier analysis with overlapping clustering by proposing an overlapping Kmeans algorithm where a fixed number of points is assigned to at least one cluster and a fixed number is not assigned at all. Alternatively, outliers can be identified by keeping track of the distances of points to their nearest center. Olukanmi and Twala (2017) determine such threshold iteratively.
A third alternative is to formulate the problem as a bi-criteria optimization, where a fixed penalty is paid (Charikar et al., 2001;Tseng, 2007) for each outlier point. Thus, the outliers can be considered to be assigned to a k + 1'th cluster, with a penalty added to the loss function (1). Tseng (2007) further give this bi-criterion approach a model-based interpretation, pointing out that the penalty term corresponds to assuming that the noise points are uniformly distributed, i.e., that they emerge from a homogeneous Poisson process.

Selection of the number of centers
Selecting the number of clusters is a widely studied topic in partitional clustering, with numerous approaches having been proposed over the years. Further, a number of reviews have been written on the topic. The review by Xu et al. (2016) focuses on clustering methods in general, while Pham et al. (2005) and Yuan and Yang (2019) focus on approaches to K-means.
In location-allocation, there are three typical approaches to determining the number of facilities. The first approach, referred to as the set covering problem, places a minimal number of facilities by guaranteeing each demand point a facility within a given distance threshold (Toregas and Revelle, 1973).
In the following two subsections, we will connect the two other approaches, namely the cost-effectiveness curve and the new facility opening cost, to the approaches utilized in clustering.

Opening cost
The number of facilities can be selected by adding a penalty to the objective function (1), referred to as the facility cost, for opening a new facility. The sums of weighted distances are then referred to as the service cost (Charikar et al., 2001). The obtained optimization problem can be referred to as a facility location problem or plant location problem (Arabani and Farahani, 2012). The cost of opening a new facility can depend on its location (see e.g. Azarmand and Neishabouri, 2009), but in many applications, it is independent of location, and a predetermined fixed price λ k is paid for each opened facility. Thus the opening cost penalty corresponds to adding a penalty term λ k k to the objective function (1).
Such a penalty for new centers has been proposed also in the partitional clustering framework (Manning et al., 2008), but the interpretation of the approach is different. The sum of distances (1) and the penalty term λ k k are then considered to correspond to distortion and model complexity, respectively. A large number of centers corresponds to a small distortion but high model complexity, and vice versa for a small number of centers. Thus, the optimization problem can be seen as scalarized bi-objective optimization, where the parameter λ k controls the trade-off between the two aims. A small λ k produces partitions with a large number of clusters, and a large λ k produces partitions with a small number of clusters. The problem is thus converted to choosing the value of λ k .
One possible approach is to utilize information-theoretic measures that trade off distortion against model complexity. For example, Akaike's information criterion (AIC) and Bayesian information criterion (BIC) are defined as AIC = argmin where l is the model maximum log-likelihood for k clusters, and q(k) is the number of parameters of a model with k clusters (Manning et al., 2008). For spatial K-means clustering with known variance this corresponds to where RSS * is the minimized residual sum of squares with k centers (8). Thus, the opening cost approach with the squared Euclidean distance metric corresponds to clustering with AIC and BIC, assuming that the data points of each cluster follow a Gaussian distribution with known σ 2 . The opening cost parameter for AIC is λ k = 4σ 2 and for BIC λ k = 2 log(n)σ 2 .
AIC implies a smaller penalty term than BIC. Correspondingly, AIC is known to favor more complex models than BIC, with a higher number of parameters (Theodoridis and Koutroumbas, 2008), corresponding here to more clusters. On the other hand, BIC tends to favor models that are too simple (Theodoridis and Koutroumbas, 2008), that is, have too few clusters.
In practice, a λ k based on neither AIC nor BIC can be blindly followed in the selection of λ k , with the ultimate choice of λ k remaining the responsibility of the model builder (Manning et al., 2008). Nevertheless, they can be used as rough guidelines for λ k , with AIC pointing towards the higher end in the spectrum of feasible λ k 's, and BIC towards the lower end. Further, mixture model based clustering can be connected to the opening cost approach if the number of clusters is assumed to be a random parameter that follows the Dirichlet distribution, i.e., if a Dirichlet process mixture model is applied Raykov et al., 2016). Furthermore, Kulis and Jordan (2012) show that if the data for each cluster is assumed to have emerged from isotropic Gaussian distributions with the same variance for each cluster, the small-variance asymptotic of a Gibbs sampler for Dirichlet process mixture model corresponds to the minimization of RSS + λ k k.

Cost-effectiveness curve
The MCLP approach draws a curve, referred to as the cost-effectiveness curve, where the horizontal axis represents the number of facilities and the vertical axis represents the percentage of population within a predefined covering area from a facility (Church and ReVelle, 1974). The number of facilities can then be chosen such that there is a reasonable trade-off between the covering percentage and the number of facilities.
This approach resembles the "elbow" approach, often used in the clustering literature (Kodinariya and Makwana, 2013). The elbow approach draws a curve with the number of clusters on the horizontal axis and the minimum of the objective function on the vertical axis. The optimal number of clusters is then selected visually as the "elbow" of the curve, where increasing the number of clusters does not appear to produce a considerable decrease in the value of the objective function.

The new tools
We introduce a novel approach that covers and unifies location-allocation and partitional clustering. The approach covers a wide range of extensions utilized in both methodologies. We refer to the proposed approach as PACK (PlAcement with Capacitated K-family 3 ). We have proposed earlier versions of PACK, with fewer extensions, in Lähderanta et al. (2020) and Lovén et al. (2020). For highlighting the novelty of this article, we summarize the extensions introduced in the different versions in Table 2.
(13) with the following constraints: 1. Center c j ∈ {p 1 , p 2 , ..., p s }, ∀j, Weights and location preference. PACK allows weighted data points and acknowledges that, regardless of their weight, it may be more important for some points than others to have a center close by. This assumption is incorporated into the problem by adding a parameter γ i > 0 to the weight of the preferred points i.e. w i = w i + γ i . Such parameter is referred to as the preference parameter. The larger the parameter γ i , the more strongly the point i attracts a center (Ackerman et al., 2012).

Capacity constraints
Continuous and discrete settings. Both settings are covered. To be precise, constraint 1 can be either omitted or applied.
Fixed and released centers. Applying constraint 2, 0 ≤ m ≤ k centers can be assigned to have predetermined locations f 1 , . . . , f m . Further, we propose a relaxation where a fixed center at a preassigned location can be released and relocated if a penalty λ f is paid. In location-allocation, λ f can be considered a price for relocating a previously constructed facility. Thus, if t fixed centers are released, the cost is λ f t. A center is released if the release decreases the value of the objective function more than λ f . Distance metrics. PACK is agnostic towards the distance metric used. The objective function (13) admits an arbitrary distance measure d.
Membership. PACK considers hard, fractional and overlapping memberships. In hard membership, y ij ∈ {0, 1} and in fractional membership y ij ∈ [0, 1] for all i, j. As was pointed out in Section 3.5, fractional membership differs from hard membership only if also capacity constraints are applied. Overlapping membership is considered from a reliability viewpoint, assigning points to multiple centres as Karatas et al. (2016). However, for avoiding ex-cess spare capacity, overlapping assignments Q i > 1 can be used only for a predetermined set of points that are considered critical in the partition.
Capacity. Both lower and upper capacity constraints can be applied for controlling the amount of weight assigned to the centers. Depending on the upper and lower limits, this constraint can be used for both balancing the data and constraining only either too large or too small clusters.
Further, as typical in location-allocation literature, PACK allows the capacity constraint and the loss function to have different constants (a i and w i , respectively). This allows, for example, applying a location preference for the objective function but not for the capacity limits.
Outliers. PACK identifies outliers by applying the bi-criterion approach (Tseng, 2007;Charikar et al., 2001). Instead of assigning a constant outlier penalty λ o for each point, we use a penalty that is proportional to the weight of the point, λ o w i . As a result, the capacity limits permitting, a point is assigned as an outlier if d(x i , c j ) > λ o . Hence, the selection depends on the distance to the clusters and not on the weight of the point.
Selecting the number of centers. PACK supports both the elbow approach as well as the opening cost approach for selecting the number of centers k.
In the elbow approach, the objective function (13) is minimized with a range of k's. The k for which the value of the objective function appears to have an "elbow" can subsequently be selected as the number of centers. However, if no capacity constraints are applied, the minimizing value of the objective function is a nonincreasing function of the number of centers. With a lower capacity limit, the objective function may however increase with a higher number of centers, as with a large number of centers, it may be difficult to obey the lower capacity limit.
The opening cost approach sets a price λ k for opening a new center. Accordingly, PACK adds a penalty λ k k to the loss function (13). PACK optimizes the penalized loss function first by minimizing the loss function without the opening cost penalty with a range of ks (see Section 4.3), and then choosing the k with the smallest value of the penalized objective function. The result depends on the value of λ k , but the relationship is not obvious or easy to predict. We propose to evade this difficulty by employing a wide range of λ k s and using the one that produces the lowest value of the objective function with most of the used values.

Location preference as a prior
Location preference refers to an additional weight, γ i , set on certain points. Such a weight can be interpreted from two points of view. First, it may be more important for some points to have a center close by than for the other points. Lähderanta et al. (2020) refer to this point-based perspective as location priority. Alternatively, centers may have a preference for some locations, due to for example some superior characteristic not represented by the weights of

Algorithm 1 PACK-algorithm
Input: x i , w i , k, n, λo, λ f , f j , j = 1, . . . , m Output: locations of the centers c * j and allocations of points to centers y * ij , j = 1, . . . , k + 1 1: for i = 1 to n do 2: Initialize c j , j = 1, 2, ..., k using K-means++ 3: while c j changes do 4: Allocation-step: minimize (13) with respect to y ij 5: Location-step: minimize (13) with respect to c j 6: S = the value of the objective function 7: end while 8: if S < S min or i = 1 then 9: S min ← S 10: c * j ← c j 11: y * ij ← y ij 12: end if 13: end for 14: return c * j , y * ij the data points. Such a preference can be seen as prior information on the location of a center.
Indeed, based on prior information, there are b locations z l , l = 1, . . . , b, which may have a center in their neighborhood. These locations may, for example, correspond to the locations of cluster centers in a previously analyzed data set. Note that b can be higher or lower than k. Assume further that the strength of this prior information can be measured with a Gaussian distribution N (z l , 1/γ 2 l I) for each l = 1, . . . , b. The higher the precision γ l , the stronger our prior belief that a center is located near z l .
Instead of building a full Bayesian posterior model, PACK employs approximate or semi-Bayesian analysis using data augmentation, where such prior information is incorporated into the model as pseudo-data (see e.g. Bedrick et al., 1996;Rhodes et al., 2016). Adding pseudo-points at z l and giving each such point a weight γ l , l = 1, . . . , b, we get where h lj indicates which cluster a pseudo-point l is assigned to. In practice, PACK does not separate the real and the pseudo-points, but rather augments the real data vector x i with the new points, and uses weights w i = w i + γ i for the augmented data, with w i = 0 for the pseudo-points not included in the real data.

The block coordinate descent algorithm with fixed k
PACK employs a block coordinate descent algorithm (Algorithm 1) for minimizing (13) with a fixed number of centers. On each iteration, a block coordinate descent algorithm minimizes the objective function with respect to a block of variables while holding other blocks of variables fixed at the values obtained in the previous iteration step (see e.g. Wright, 2015). The two main steps, i.e., the blocks, for the algorithm are the allocation-step (line 4), where the points are assigned to the centers, and the location-step (line 5), where the centers are relocated based on the points assigned to them. The steps are iterated until convergence is reached.
Further, as block-coordinate descent finds the local minima close to the initial values, PACK uses a number of initial values to find the global minimum.
The initial values are obtained using the K-means ++-algorithm (Arthur and Vassilvitskii, 2007) that spreads the initial locations of centers improving both the speed and the accuracy of the K-means method.
Allocation step. The allocation step minimizes the objective function (13) with respect to y ij . Constraints 3, 4 and 5 are applied, and the locations of the centers c j are assumed to be fixed. If no capacity constraints are applied, this step corresponds to assigning each point to the nearest facility measured in terms of the chosen distance metric. In addition a point is assigned as an outlier if its distance to the nearest center exceeds the limit λ o .
If capacity constraints are applied, this step is an NP-hard integer programming task when the hard membership constraint (2) is employed, or a linear programming task when the fractional membership constraint (12) is employed. Fractional membership can thus be considered a linear relaxation, facilitating computations. Thus, in some applications, fractional membership can be used as an approximation of the hard membership as it enables faster computations and the utilization of larger data sets.
Interpretation of the allocation step depends on the framework. In clustering, allocation corresponds to determining the cluster where the point appears to have emerged from. On the other hand, in location-allocation, allocation corresponds to determining which facility is the best for serving each demand point.
Location step. Location step minimizes the objective function (13) with respect to c j 's, while keeping the allocations y ij fixed. In other words, each center is separately relocated given the points assigned to it. This step is omitted for fixed centers.
In continuous space with the squared Euclidean distance metric, c j is the weighted mean of the points allocated to cluster j. With the Euclidean distance metric, c j is the geometric median and it is computed using Weiszfeld's algorithm (Weiszfeld, 1937;Kuhn, 1973). In the discrete setting, c j is the point which minimizes the sum of pairwise distances among points assigned to c j .
If releasing fixed centers is allowed, a center is released and relocated if the sum of distances to the assigned points is reduced more than λ f . Similarly, a previously released fixed center can be returned to its original location if the increase in the sum of distances is lower than the releasing penalty λ f .
In model-based clustering, location step corresponds to estimating the parameters of the underlying probability distributions and further, in partitional clustering, determining the best representator for the assigned points. In location-allocation, location step corresponds to determining the best location for a facility that serves a given set of demand points.
Implementation. The algorithm is implemented as an open source R package called rpack (Lähderanta et al., 2019) available at GitHub. The allocation step is run on Gurobi (Gurobi Optimization, 2018), a fast optimizer package with R bindings freely available for academic use. If Gurobi is not available, PACK uses the lpSolve-package for R (Berkelaar and others, 2015) for the allocation step.

Evaluation
We demonstrate the feasibility of PACK with three experiments.
The first experiment uses an artificial data set to demonstrate the effects of capacity limits and outliers while clustering with the Euclidean and the squared Euclidean distance metrics.
The second experiment demonstrates the use of the opening cost, predetermined server locations, releasing fixed centers, location preference, and selecting k. The experiment places edge computing servers in the city center of Shanghai, China Wang et al., 2019;Xu et al., 2019).
The third example compares the effects of hard and fractional memberships. The experiment places a predetermined number of recycling centers in the town of Partizanske in Slovakia (Cebecauer and Buzna, 2018).

Experiment 1: clustering with artificial data points
This example shows how the capacity constraints guarantee that the workload of different clusters is reasonable. Further, the example illustrates how outliers are identified. We also evaluate different clusterings with a known cluster structure, using the adjusted Rand index (ARI) (Rand, 1971) measure.
We simulate 500 data points which form 10 clusters as follows: 1. The number of data points in the simulated clusters is set as 20, 20, 40, 40, 50, 50, 60, 60, 80 and 80. The uneven cluster sizes illustrate how the results change with and without the capacity limits. 2. Data points are simulated from a set of multivariate gamma distributions constructed using a normal (Gaussian) copula utilizing the R package lcmix (Xue-Kun Song, 2000). The shape parameters of the gamma distributions vary between 0 and 15, and the scale parameters vary between 0 and 100. These parameter values will generate both asymmetric and symmetric clusters. 3. The weights of each data point are generated independently from the uniform distribution Unif[1, 100]. To make clusters less homogeneous in the terms of the weights, about half of the clusters have small weights at the cluster center and heavy weight points at the edge of clusters. These weights are directly proportional to the Euclidean distance between the cluster center (the mean point of the multivariate gamma distribution) and the data point. 4. Clusters that consist of independently generated data points, drawn from a set of multivariate gamma distributions, frequently overlap. Separating clusters from each other is thus difficult because cluster shapes are unclear.
To reduce the overlap of the resulting clusters, the mean of each multivariate gamma distribution is shifted randomly on a predefined grid. Moreover, we use an additional scaling parameter to draw data points closer to the corresponding cluster center to further reduce the overlap of data points. See Fig. 1 for a schematic illustration of this step. 5. Additional outlier points are sampled uniformly from the grid determined by the outermost data points. Like the weights of the "true" data points, the weights of the outlier points are sampled independently from the uniform distribution Unif [1,100]. Overall, 20 outlier points are added to the data.
The simulation algorithm we use to generate data points from multivariate Gamma distribution(s) is included in the rpack R-package (Lähderanta et al., 2019). We set the number of centers k as 10. The capacity limits of each cluster are set as [1587,3587], that is, [w − 1000,w + 1000], wherew = n i w i /k = 2587 is the empirical mean of the cluster weights. We set the outlier penalty parameter λ as 0.05 with the squared Euclidean distance metric and as 0.2 with the Euclidean distance metric to demonstrate how identifying outliers change the resulting clusters. For each data replicate we run the Algorithm 1 with 50 initial values, selected with the k-means++ algorithm. Clustering results are illustrated in Figs. 2 and 3. Fig. 2: Continuous clustering of 500 data points with 10 clusters and the squared Euclidean distance metric, with capacity limits (A) and without capacity limits (B). The black crosses denote estimated cluster centers.
Further, we illustrate how capacity constraints and the outlier penalty change the clustering of data points determined by latent factors, that is, different multivariate gamma distributions and their parametrizations. We run several data simulations and each time examine how the estimated clusters compare with the densities of the underlying multivariate gamma distributions providing the simulated data points. We call these densities the ground truth clusters. We use the adjusted Rand index (ARI), computed with the R package fpc (Vinh et al., 2010), as a diagnostic metric in this simulation example. The closer the ARI is to one, the closer the clusters are to the ground truth clusters. In this example, capacity limits are determined by the empirical mean of the cluster weights, and we use the same outlier penalty parameter values as described above. Averaged ARI measures are determined based on 50 simulation replications. The averaged ARI results are presented in Fig. 4. Finally, we note that the averaged ARI values cannot be exactly one due to the overlap of the ground truth clusters.

Artificial data points clustering results
Fig. 2 B clearly demonstrates that without capacity limits the resulting cluster sizes are drastically uneven: the smallest cluster size is 259, whereas the largest is 5211. Further, there is a very large cluster with a concentration of data points. Using capacity limits, the concentration of data points is divided into more even-sized clusters (Fig. 2 A).
Further, Fig. 3 demonstrates how excluding the oultier penalty parameter results in very broad clusters. For example, the identified cluster found in the right upper corner in Fig. 3 A is very wide and it would be reasonable to Fig. 3: Continuous clustering of 500 data points with 10 clusters and the standard Euclidean distance metric, with the outlier penalty (A) and without the outlier penalty (B). The black crosses denote estimated cluster centers.
consider the outermost data points as outliers. When the penalty is utilized, these outermost points are consistently determined as outliers.
The squared Euclidean distance metric, combined with the outlier penalty, seems to produce the largest averaged ARI values (Fig 4). There is no striking difference between the averaged ARI measures while capacity limits are used in the clustering algorithm and when the same data is clustered without the capacity limits. Overall, the inclusion of the outlier penalty and the capacity limits did not produce clusters which would be drastically in contradiction to the mechanism by which the data points are generated. This is a good property because data points are not forced into clusters which would not reflect the true data while capacity constraints are applied.
The difference between the squared Euclidean and the Euclidean distance metrics is very subtle in the continuous setting. Although the resulting clusters are almost identical, it seems the squared Euclidean metric favors more spherical cluster structures. With the Euclidean metric, large clusters of the original data are divided into smaller, less spherical like clusters.

Experiment 2: Edge server placement
The second experiment optimizes the telecommunications infrastructure in Shanghai, placing edge computing servers to reduce the communication latency experienced by the mobile phone users. Edge computing (Shi et al., 2016) refers to computing infrastructure that facilitates data processing directly at the points of interest. Deploying such infrastructure, characteristics such as radio coverage, the available telecommunications networks, and users' need 2 ) in the loss function are illustrated with filled boxplot and clusters determined using Euclidean distance (L 2 ) in the loss function are illustrated with empty boxplots. Notations "+ LU" and "+λ out " denote, respectively, if capacity limits (Lower and Upper) or the outlier penalty were applied in PACK.
for computational capacity, largely dictate the physical placement options of the servers .
The experiment is based on the Shanghai Telecom data set Wang et al., 2019;Xu et al., 2019), which contains the locations of 2732 base stations in the Shanghai region and the mobile user connections to those base stations in a six-month period.
The goal of edge server placement is to deploy a number of edge computing servers in the city region such that the latencies between the servers and the base stations assigned to each server are minimized. Each base station is assigned to exactly one edge server, and the workload (i.e., weight) of the base station, that is, its maximum number of concurrent users in the recorded data, is allocated to that server. Spatial distribution of the base stations and their workloads are shown in Fig. 5.
This minimization problem can be interpreted as a capacitated locationallocation problem, as each server, represented by a cluster center, has a limited computing capacity for the workload it can handle, represented by the sum of the weights of the base station locations in the corresponding cluster. Further, each edge server can only be placed at the same location as one of the base stations. Since the topology of the underlying network is unknown, we use geospatial distances to approximate the latencies between the base stations. Edge server placement provides an ideal setting for studying two particular PACK extensions, as these extensions have clear real-world interpretations in the scenarios included: 1. Choosing k, or the optimal number of edge servers to be placed on the region. 2. Predetermination of center locations, placing a predetermined number of edge servers in the region, with ten of those servers already deployed or having preferred locations.

Choosing the number of edge servers
We first determine the optimal number of edge servers to be deployed. We apply the opening cost approach by selecting five different values for the λ k penalties, one of which is zero, and examine the value of the objective function with each possible value of k. We let k range from 35 to 45 and with each k we use 100 different initial values, chosen with k-means++. We assume the edge servers are identical and have an upper limit for the workload they can manage. Further, we set a lower limit for each server to balance the workload between the servers and to ensure no server capacity is wasted. We set the limits of one server to [349,648] which corresponds to a wide capacity range for the server. These limits correspond to [x−150,x+150], withx = 498.875 the average workload for the servers in the k = 40 scenario.
As the value of the opening cost penalty λ k is unknown in this experiment, we use five different values of λ k for each value of k. The resulting values of objective function allow us to evaluate the robustness of the deployment with respect to the penalty parameter and to decide accordingly which amount of servers is best suited for the data.
The values of the objective function with the given variables, that is, the penalties and the number of centers, are listed in Fig. 6. As expected, when no opening-cost is set, the value of the objective function decreases as k increases. An exception, resulting from the lower capacity limit making it difficult to place servers optimally, can be seen at k = 45, with the objective function value increasing in all five penalty cases. An "elbow shape" can be seen at k = 38 with all penalty values, which indicates that number of servers provides a robust clustering. Another possible clustering choice is at k = 44, as it seems to have the lowest value of the objective function with small penalty values.

Location-allocation with predetermined server locations
In this scenario, we place a total of 38 edge servers in the region, with ten of the server locations predetermined in four alternative ways, corresponding to the following setups: 1. No pre-determined cluster center locations. 2. The locations of the ten predetermined server locations are fixed. This corresponds to the deployment of 28 servers, with ten servers already deployed. 3. There's a release cost to the ten servers such that if the benefit (i.e., reduction in the objective function (13)) of relocating any of those servers is greater than the release cost, then that server is relocated. 4. A location preference is set for the ten predetermined server location candidates such that those locations are considered "more important" and that they attract the servers more than the other locations. Location preference in edge server placement would correspond to a situation where, for example, mobile connectivity in a certain area would be above that of the surrounding areas and especially good at the center of that area. Such an area, and especially its center, would thus be a preferred location for an edge server. The value of the additional weight is γ i = 200 for each preferred server location i = 1, . . . , 10.   7 compares these four different setups. Overall, each setup produces deployment of servers with some important differences. First, the benefit of relocating the fixed servers is clear: some of the fixed servers in the southwest and in the south-east are located on the edge of the map which overall worsens the average distance to the center, that is, latency to the server, as can be seen in Table 4. This average latency is improved as we allow the servers to be relocated (setup 3). Setup 4 can be seen as a "relaxed" version of the setup 3, as it pulls the servers closer to the predetermined locations rather than forcing the servers to be located exactly at the predetermined locations. Furthermore, the mean and median distances in setup 4 are much smaller than in setups 2 and 3. Fig. 7: The optimal placement of 38 edge servers in the Shanghai region with ten fixed servers. The edge servers are marked as crosses and locations of fixed servers as squares. On the top left, the clustering is done with no predetermined server locations (setup 1). On the top right, predetermined server locations are fixed (setup 2). On the lower left, we allow the predetermined locations to relocate if the benefit is greater than the release cost (setup 3). On the lower right, we add a location priority to the predetermined server locations (setup 4).
All the four setups produce clusterings with no unnatural cluster structure. The differences in the resulting clusters are mainly due to the locations of the predetermined spots: forcing the servers to the fixed locations can cause worse average latency in the deployment, which can be solved by allowing servers at the fixed locations to relocate. Fig. 8: The optimal placement of 12 recycling centers in the region of Partizánske. Left: the data points and the center location candidates shown as red crosses. Right: The partitioning with fractional membership and capacity limits for each center set to [2963,4962].

Experiment 3: Recycling center placement
The third experiment places a fixed number of recycling centers in the Slovakian city of Partizánske such that the average road distance between residents and the centers is minimized. The data set (Cebecauer and Buzna, 2018) includes a distance matrix, listing the road distances between the demand points, and the locations of those demand points. Further, each demand point has a residential population associated with it, which corresponds here to the weight of the point.
The experiment applies two PACK extensions. First, as only the resulting recycling center locations are in focus, and not the allocations of the residents, we apply fractional membership to speed up the optimization. Second, since placing a recycling center anywhere is not possible, we restrict the possible placement to 50 predefined locations in the city. These locations as well as the demand points can be seen in the left panel of Fig. 8.
To study the effect of the fractional membership on the speed of the algorithm, we analyze the impact of two different capacity limits on the placement in both fractional and hard membership scenarios. We set the wide capacity limits to [2963,4962] and the tight capacity limits to [3962,3963].
The right panel of Fig. 8 illustrates the resulting clustering, obtained with fractional membership and wide capacity limits. The overall clustering seems to be reasonable and the recycling center locations are spread out. There is some overlapping of the clusters in the north-west region of the city. This is due to the underlying distance matrix, where some of the points have exactly the same distance to a certain point, even though the distances appear different  based on the map. In the city center, the centers are close to each other, which is partly due to the distribution of the population. The difference between the clusters resulting from fractional and hard membership is minor based on the mean, median an 95% quantile distances (Table 5) and based on the clustering results on the map (result not shown). However, the fractional membership scenarios are much faster to run than the hard membership scenarios. With the wide capacity limits, the algorithm was approximately 3 times faster with fractional membership and with tight limits approximately 10 times faster. Furthermore, we tested even tighter limits that allow only one size for clusters. This problem was easily solved with fractional membership, but it was not solvable with the hard membership constraint, as it was not possible to distribute the weights of the demand points evenly to each cluster without dividing individual weights.

Discussion
In this article, we mapped the connections and the different interpretations between location-allocation and partitional spatial clustering and highlighted their immediate links to model-based clustering. Specifically, we considered how these frameworks consider different extensions such as capacity constraints, outliers, different membership types, and the inclusion of prior information. In addition, we introduced the improved PACK algorithm as a flexible tool for both frameworks. We proposed new extensions to PACK such as outliers and releasing a fixed center. While PACK is a versatile tool that can be used in various domains, here we demonstrated its feasibility in the topical domains of location intelligence and telecommunications, and explored the effects of its different extensions with both real and artificial data sets.
With the artificial data set, we demonstrated that while the capacity constraints reduced the overall performance of the clustering, the resulting clusters were much more balanced than with no capacity constraints. Further, while obeying capacity constraints, PACK did not produce clusters that would drastically contradict the natural groupings of the data, with the squared Euclidean distance metric producing slightly more spherical clusters than the Euclidean distance metric.
Finally, PACK was able to identify true outliers even though they were not prominent in the data. Identifying the outliers did not change the cluster centers markedly, as the simulated clusters were relatively large. With smaller clusters, the effect of one point would be larger, implying that outliers would change the clustering results more.
With the edge server placement example, we demonstrated the selection of the number of centers, the effect of releasing a fixed center, and the prior for the locations of the centers. Testing a wide range of penalty values and applying the one that had the lowest, or relatively low, value of the objective function with most of the penalties, identified two prominent alternatives for the number of clusters, one of which was chosen as the final number of clusters. Such a range of penalty, tuning or smoothing parameter values has been proven to be a useful parameter selection method also in various other optimization and estimation models (Holmström and Pasanen, 2017;Lee et al., 2012;Hastie et al., 2009).
Both releasing fixed enters and the location preference prior provide means to incorporate prior information about the locations of the centers. However, the extensions had different effects on the partition, with location preference favoring the neighborhoods of the given points, while releasing fixed centers emphasized only the given points, either locating a center exactly on them or allowing a free placement for the center.
With the recycling center placement experiment we focused on the effect of fractional membership. Placing recycling centers, the resulting center locations with both hard and fractional memberships were relatively similar. Fractional membership provided a considerable speed up in the computations, especially with tight capacity limits. Further, with extremely tight capacity limits, the problem was not even solvable applying hard membership.
Given these computational benefits, fractional membership may be preferred as an approximative method also in settings where a hard membership would be more appropriate. Each point whose membership is divided between a number of centers could, in such a setup, be assigned to the center with the highest membership. However, the downside of such an approximation is that after the rounding, the capacity limits may not be satisfied.

Conclusion
Our survey revealed that there are very similar approaches proposed for different extensions and constraints in partitional clustering and location-allocation. Highlighting the intersection between partitional clustering and location-allocation, along with the drawn parallels to model-based clustering, gives valuable insights into the implicit assumptions of the used models and helps unify the seemingly different frameworks. Further, bringing forth these connections opens up new pathways to both partitional clustering and location-allocation research.
Extending our earlier work, we proposed a novel algorithm for clustering, covering an extensive number of extensions common both in partitional clustering and location-allocation, and implemented the algorithm as a readilyavailable software tool for conducting clustering analyses. We argue that the flexibility of the PACK algorithm, combined with the easy-to-use open-source R software package, provides a versatile toolbox for both spatial clustering and location-allocation.

Declarations
Funding. This research is supported by Academy of Finland 6Genesis Flagship (grant 318927), the Infotech Oulu research institute, the Future Makers program of the Jane and Aatos Erkko Foundation and the Technology Industries of Finland Centennial Foundation, by Academy of Finland PROFI5 funding for mathematics and AI: data insight for high-dimensional dynamics, and by the personal grant for Lauri Lovén on Edge-native AI research by the Tauno Tönning foundation.
Conflicts of interest. The authors declare that they have no conflict of interest.
Availability of data and material. Data sets similar to those of Experiment 1 can be reproduced based on sripts available at https://github.com/terolahderanta/rpack. The data set of Experiment 2 is available at http://sguangwang.com/TelecomDataset.html and the data set of Experiment 3 at http://frdsa.uniza.sk/∼buzna/page5/supplement5/data paper.html.