INFORMATICA Informatica 1822-88440868-49520868-4952 Vilnius University INFOR610 10.15388/25-INFOR610 Research Article Global Optimization Algorithm for the Multi-Weber Problem with Polyhedral Barriers KepalasMindaugasmindaugas.kepalas@mif.stud.vu.lt∗

M. Kepalas is a PhD student at the Institute of Data Science and Digital Technologies, Vilnius University, Lithuania. His research interests are global optimization, location problems with constraints, net-constrained clustering, mathematical programming, optimization algorithms and applications.

 ŽilinskasJuliusjulius.zilinskas@mif.vu.lt

J. Žilinskas is a research professor and a head of Global optimization group at the Institute of Data Science and Digital Technologies, Vilnius University, Lithuania. His research interests are global optimization, parallel computing, data analysis and visualization. He published around 100 international journal papers and 4 books. He is a member of editorial boards of 10 international journals. He is a vice president of the International Society of Global Optimization.

 Vilnius University, Faculty of Mathematics and Informatics, Institute of Data Science and Digital Technologies, Akademijos st. 4, Vilnius, LT-08412, Lithuania Corresponding author. 20252111202536487590282025112025 © 2025 Vilnius University2025 Open access article under the CC BY license.

In this paper, we consider the multi-Weber problem with polyhedral barriers. For this problem, a set of obstacles are introduced where travelling or placement is prohibited, which makes the distance metric non-convex and requires constructing a special graph for calculating the distances between pairs of points. For obtaining the global solution of the problem, we build a branch and bound algorithm with pruning criteria based on dividing clients into groups and analysing them separately. We have managed to obtain global solutions to several multi-Weber with polyhedral barriers problem size instances which to our knowledge have not been reported before.

 multi-Weber problem with barriers clustering problems with centre location constraints global optimization Introduction

This paper considers the problem of determining the global solution of a general continuous multi-locations problem, where a set of clients (i.e. points at given coordinates) have to be assigned to a set of facilities, for which the optimal locations have to be determined. In the paper, the distance from a client to a facility is assumed to be arbitrary in a sense that it is not necessarily measured by the standard Euclidean distance between two points; in particular, in our numerical part we consider the multi-Weber with polyhedral barriers problem. For this problem, a set of obstacles are introduced where travelling or facility placement is prohibited, which makes the distance metric non-convex. Also, we assume that the facilities could in addition be restricted to an arbitrary non-convex finite subset of the plane (e.g. this could be a union of a set of segments and triangles).

At optimality for the problems which are studied in this paper, a client always chooses the closest facility, thus the problem can be also seen as a clustering problem: the goal is to find the optimal partition of the clients into groups (clusters), and then the optimal facility can be determined for each cluster independently. Thus, for obtaining the global solution of the multi-locations problem, we employ a branch and bound with an all-possible-partitions enumeration tree and pruning criteria. The pruning criteria roughly uses the following simple idea: if we divide the clients into two groups and find the optimal solution within each group separately, summing the two terms we obtain a lower bound estimate for the partitions within a corresponding branch.

The Weber problem, i.e. the problem of determining the coordinates of the optimal point minimizing the sum of the Euclidean distances to other given points, has a long and rich history, many extensions and publications (Drezner et al., 2002). The extension we concentrate on in this paper was probably firstly introduced by Aneja and Parlar (1994). In this model, obstacles (barriers) in the plane are introduced; it is not allowed to travel through these barriers or to place a solution point inside. In real life, such barriers could be considered to represent lakes or buildings on a geographical map, or some prohibited places in printed circuit boards (Hamacher and Nickel, 1995).

The Weber location with barriers problem has attracted the attention of many scientists. However, most of the work concentrates on solving the single location problem, see e.g. Aneja and Parlar (1994), Klamroth (2002), Bischoff and Klamroth (2007), Oğuz et al. (2016). The first reason for this is the following motivation: having a “good” algorithm for the single location problem, we can implement a location-allocation algorithm in the style of Cooper (1964). For example, such an algorithm for multilocations with polyhedral barriers problem is considered by Bischoff et al. (2009). The second reason why the case of multilocations has not received much scientific effort might be related to the combinatorial complexity of this problem.

The only non-theoretical work dedicated to the multi-Weber with barriers problem we have found in the scientific literature is the thesis of Krau (1997). In his work, the author extends the column generation principle of Gilmore and Gomory (1961) to the problem. The column generation approach seems to provide a benchmark solution for several clustering-like problems such as the minimum sum of squares clustering (Aloise et al., 2012) or the (barrierless) multi-Weber problem (du Merle et al., 1999). Using this method, Krau (1997) solves the problem presented by Aneja and Parlar (1994) which consists of 18 points and 12 barriers, and also an extended version of this problem with 30 points, for cluster sizes K = 2 , … , 7.

As the results indicate, the cases with small K seem to be the hardest for a column generation based algorithm (Krau, 1997). However, with our approach we have solved the problem instances of 60 / 2, 45 / 3, 40 / 4 and 35 / 5 ( N / K means “N points into K clusters”) to global optimality (within ϵ-accuracy with ϵ = 1 %). To our knowledge, this is the first time the solution of instances with such a size to global optimality is reported. Moreover, we think our approach is very general. We believe the ideas could be also extended to other variants of multi-Weber problems.

Our paper is structured as follows. In Section 2, we formulate the general model of the problem we seek to solve in this paper. In Section 3, we describe the building blocks of our approach and give some mathematical proofs of intuitively rather obvious facts. Using these blocks, in Section 4 we present the sketch of the enumeration algorithm. Finally, in Section 5, we clarify algorithm details for the multi-Weber with polyhedral barriers problem and give an analysis of the performance of the algorithm.

 General Multi-Locations (Clustering) Problem

Let us formulate the following general multi-locations (clustering) problem: min ( C 1 , , μ 1 ) , … , ( C K , μ K ) ∑ k = 1 K ∑ i ∈ C k d ( p i , μ k ) s.t. μ k ∈ R ∪ , k = 1 , … , K . Problem input:

N “points”, “clients” p 1 , p 2 , … , p N in the problem space S: these are the elements we want to group into clusters. In our numerical part for the multi-Weber with barriers problem, these are points in R 2 , but the theoretical ideas could also be applied to points in an arbitrary R D with dimension D > 2, or even other elements from some more complex space S, e.g. some function space.

K – the number of clusters.

R ∪ – the constrained set for cluster centres. This set is a union of M “containers”, “cells” or “region-units” R 1 , R 2 , … , R M : R ∪ = ⋃ m = 1 M R m . See further comments on R ∪ in Section 2.1.

Distance metric (function) d : S × S ↦ R + to measure the distance between a “client” p i and a “facility” μ k . In this paper, we concentrate on a multi-location with barriers problem, where obstacles (barriers) in the plane are introduced. Any allowed path between points p i and μ k must not cross any of the barriers.

 Problem variables:

Cluster centres (“facilities”) μ k , constrained to be placed within the set R ∪ : μ k ∈ R ∪ , k = 1 , … , K .

Corresponding clusters1

A partition of the set N 1 N = { 1 , 2 , … , N }: ⋃ k = 1 K C k = N 1 N , C k ∩ C l = ∅, k ≠ l.

C k , k = 1 , … , K. These are the sets of “clients” assigned to “facility” at μ k . Somewhat, such variable is redundant, because in our problem formulation each “client” is assigned to the closest “facility”, thus the clusters are determined given μ k , k = 1 , … , K. However, this notion will be useful when considering the k-means-type (location-allocation) algorithm discussed further and also for the branch and bound enumeration algorithm presented in the paper.

 Our goal is to find the global solution to problem (1). Using notation P n 1 n 2 = { p n 1 , … , p n 2 } , for n 1 ⩽ n 2 , we label problem (1) with [ P 1 N / K ] | [ R ∪ ] for short.2

In words, “cluster the point set P 1 N into (/) K clusters given a condition (|) that cluster centres are restricted to be within R ∪ ”.

 Correspondingly, we also label the globally optimal (minimal) value of (1) with Loss opt ( [ P 1 N / K ] | [ R ∪ ] ).

Constrained Set <inline-formula id="j_infor610_ineq_038"><alternatives><mml:math> <mml:msup> <mml:mrow> <mml:mi mathvariant="script">R</mml:mi> </mml:mrow> <mml:mrow> <mml:mo>∪</mml:mo> </mml:mrow> </mml:msup></mml:math><tex-math><![CDATA[${\mathcal{R}^{\cup }}$]]></tex-math></alternatives></inline-formula>

In the most general case, we assume that the constrained set R ∪ is a (finite) union of arbitrary closed and bounded subsets (“region-units”) of the problem space S: R ∪ = ⋃ m = 1 M R m , R m ⊂ S , m = 1 , … , M . The bounded property is not restrictive for the problem, since the set { p 1 , … , p N }, as a finite set of points, is always bounded, and it can be shown that this implies that all optimal centres μ k for the problem (1) always lie within some bounded subset of S. In addition, the closed property for each region-unit R m ensures that for any p ∈ S, we can always find the nearest and the furthest points μ m inf ∈ R m and μ m sup ∈ R m , such that d ( p , μ m inf ) = inf μ ∈ R m d ( p , μ ) ⩾ 0 , d ( p , μ m sup ) = sup μ ∈ R m d ( p , μ ) < ∞ .

As a less general example (for the multi-locations with polyhedral barriers problem we study in the numerical part of this paper), we define R ∪ to be a union of (arbitrary) polygons. We note that any such set R ∪ of (possibly, overlapping) polygons can be decomposed into a union of (non-overlapping) triangles by some triangulation method,3

See for example the implementation in [R] sfdtc package (Sumner, 2025).

where by “non-overlapping” we mean that all the triangles will “touch” only along the sides. Such decomposition is also possible if some of the region-units R m , m = 1 , … , M are segments; then R ∪ can be decomposed into a union of “non-overlapping” segments and triangles.

Going to higher dimensions, e.g. for S = R 3 , R ∪ could be decomposed into a union of segments, triangles and tetrahedrons (Paulavičius and Žilinskas, 2013), etc., we can continue to higher dimensional Euclidean spaces. This “decomposition into simplices” is a separate research topic and we do not consider this question in our paper. Therefore, in our numerical part for the multi-Weber with barriers problem we assume that R ∪ = ⋃ m = 1 M R m is provided as a union of triangles R m , m = 1 , … , M.

 Alternative Problem Formulation and an Iterative Algorithm to Obtain a Local Solution

For an arbitrary cluster C k ⊆ { 1 , 2 , … , N } and an arbitrary centre μ k ∈ S, define the related loss for this cluster-centre pair ( C k , μ k ): Loss ( C k , μ k ) = ∑ i ∈ C k d ( p i , μ k ) . With this definition, the loss in (1) can be written this way: ∑ k = 1 K ∑ i ∈ C k d ( p i , μ k ) = ∑ k = 1 K Loss ( C k , μ k ) . As can be seen from the above equation, the optimal centres given a predefined partition of the clients into clusters can be found independently. Therefore, (1) can be rewritten this way: min C 1 , … , C K ∑ k = 1 K min μ k ∈ R ∪ Loss ( C k , μ k ) . Formulation (3) can be noticed to lead us to the classical location-allocation (Cooper, 1964) (k-means-type) algorithm:

Given the (fixed) clusters C 1 , … , C K , determine (update) optimal cluster centres: μ k = arg min μ ∈ R ∪ Loss ( C k , μ ) , k = 1 , … , K .

For the (fixed) cluster centrees μ 1 , … , μ K , determine (update) optimal clusters: C k = { i : d ( p i , μ k ) < d ( p i , μ l ) , ∀ l ≠ k } , k = 1 , … , K .

Repeating the procedure above, the loss (3) decreases with every step, and at the termination we are guaranteed to arrive at a local solution of problem (1), i.e. at a solution which cannot be improved by switching cluster memberships or changing cluster centres. For the multi-Weber with polyhedral barriers problem, see, for example, such a location-allocation algorithm by Bischoff et al. (2009).

 Building Blocks of the Algorithm Enumeration Tree

To find the global solution of problem (1), we could theoretically enumerate all the possible partitions of the index set N 1 N = { 1 , … , N } into K clusters C 1 , … , C K , find the optimal centres for each possible clustering, and pick the best solution.

Such an approach would be brutal and would allow to solve only very small problems, as there are S 2 ( N , K ) of such partitions into clusters.4

Here, S 2 ( ∙ , ∙ ) is Stirling’s number of the 2nd kind. This number satisfies the recursive relation S 2 ( N + 1 , K ) = K S 2 ( N , K ) + S 2 ( N , K − 1 ) , and from this formula one can calculate that, e.g. S 2 ( 18 , 4 ) = 2 , 798 , 806 , 985 – already a very big number. Furthermore, from (4) one can see that adding any additional element to the clustering problem (and keeping the number of clusters fixed at K) increases the number of partitions more than K times, i.e. this number grows exponentially.

Nevertheless, we have shown in our recent article (Kepalas and Žilinskas, 2024) on a related problem that having a good pruning criteria, problems of size 50 / 7 (i.e. N = 50 points into K = 7 clusters) become possible to solve using this approach. We will discuss the pruning criteria in Section 3.2. For now, consider the tree T ( N , K ) in Fig. 1, which visualizes such an enumeration procedure.

Clustering enumeration tree T ( N , K ) for the algorithm.

In the tree, the n-th level (distance from the root v 1 ) of the node corresponds to the n-th element in the set { p 1 , … , p N }, and the label λ of the corresponding node represents the cluster (index) to which this element is assigned along the path. As an example, the highlighted path v 1 v 2 v 3 v 4 v 5 in Fig. 1 results in a sequence λ ( v 1 ) = 1 , λ ( v 2 ) = 2 , λ ( v 3 ) = 2 , λ ( v 4 ) = 1 , λ ( v 5 ) = 3 , which corresponds to clustering C 1 = { 1 , 4 } , C 2 = { 2 , 3 } , C 3 = { 5 } . In general, the enumeration procedure is as follows:

Start with the first element in the first cluster C 1 = { 1 } and the remaining clusters empty: C k = ∅, k = 2 , … , K.

When assigning the next element to a cluster (i.e. when going to the next level), you are only allowed to put the corresponding element

into a non-empty cluster or

into an empty cluster with the lowest index.

By two rules, the 2nd element can be assigned to either C 1 or C 2 .

For the 3rd element, C 3 might not be available (in case C 1 = { 1 , 2 }, C 2 = ∅; see the nodes at the 3rd level of T ( N , K )).

For the 4th element, see the nodes at the 4th level of T ( N , K ) and etc.

Terminating at the leaf-node v N ∈ T ( N , K ), we obtain clustering C 1 = { 1 , … } , C 2 = { … } , … , C K = { … } by tracing the labels in the path from the root-node to the leaf-node.

 It can be shown by simple arguments that this way the nodes v n at the level n enumerate all the possible partitions of the set { 1 , 2 , … , n } into no more than K groups. We notate such partition associated with any vertex v n ∈ T ( N , K ) with P ( v n ), i.e. P ( v n ) = { C 1 ( v n ) , … , C K ( v n ) }.5

C k ( v n ) = C k ( v 1 v 2 … v n ) = { i ∈ { 1 , 2 , … , n } : λ ( v i ) = k }.

 When n < N, to emphasize that not all N points are yet assigned to some cluster, we sometimes call such a partition “partial” and add the symbol ∂ to our notation, i.e. ∂ P ( v n ) = P ( v n ) for n < N. Also, we notate the set of all the vertices at the level N with L ( N , K ) (this is the leaf-node set of the tree).

 Pruning Criteria

The pruning criterion is based on the following simple observation.

(<italic>Lower bound by division into two groups</italic>)<italic>.</italic>

Consider a division of the point set P 1 N = { p 1 , … , p N } into two groups P ˜ 1 , P ˜ 2 s.t. P ˜ 1 ∩ P ˜ 2 = ∅ and P ˜ 1 ∪ P ˜ 2 = P 1 N . Then, Loss opt ( [ P 1 N / K ] | [ R ∪ ] ) ⩾ Loss opt ( [ P ˜ 1 / K ] | [ R ∪ ] ) + Loss opt ( [ P ˜ 2 / K ] | [ R ∪ ] ) .

Define index sets N ˜ 1 = { i : p i ∈ P ˜ 1 }, N ˜ 2 = { i : p i ∈ P ˜ 2 }. Suppose the globally optimal value for problem [ P 1 N / K ] | [ R ∪ ] is attained for clustering C 1 , … , C K with centres at locations μ 1 , … , μ K . Also label C k † = C k ∩ N ˜ 1 , C k ‡ = C k ∩ N ˜ 2 , k = 1 , … , K . We have Loss opt ( [ P 1 N / K ] | [ R ∪ ] ) = ∑ k = 1 K ∑ i ∈ C k d ( p i , μ k ) = ∑ k = 1 K ∑ i ∈ C k † d ( p i , μ k ) + ∑ k = 1 K ∑ i ∈ C k ‡ d ( p i , μ k ) ⩾ ∑ k = 1 K [ min μ ∈ R ∪ ∑ i ∈ C k † d ( p i , μ ) ] + ∑ k = 1 K [ min μ ∈ R ∪ ∑ i ∈ C k ‡ d ( p i , μ ) ] ⩾ Loss opt ( [ P ˜ 1 / K ] | [ R ∪ ] ) + Loss opt ( [ P ˜ 2 / K ] | [ R ∪ ] ) . The last inequality is justified by the reasoning that if we “relax” the partition of the set N ˜ 1 from being “fixed” to C 1 † , … , C K † , we can only obtain an improvement. (Same applies to N ˜ 2 and C 1 ‡ , … , C K ‡ .)  □

Take some n < N. Then, Loss opt ( [ P 1 N / K ] | [ R ∪ ] ) ⩾ Loss opt ( [ P 1 n / K ] | [ R ∪ ] ) + Loss opt ( [ P n + 1 N / K ] | [ R ∪ ] ) .

Label Loss opt ( { p i , i ∈ C k } | R ∪ ) = min μ ∈ R ∪ ∑ i ∈ C k d ( p i , μ ) . Consider some v N ∈ L ( N , K ) and the corresponding partition P N = { C 1 ( v N ) , … , C K ( v N ) } . Define Loss opt ( [ P 1 N / K ] | [ R ∪ , P N ] ) = ∑ k = 1 K Loss opt ( { p i , i ∈ C k ( v N ) } | R ∪ ) . Now suppose we are at some vertex v n ∈ T ( N , K ) at the level n < N, i.e., we are in the path v 1 v 2 … v n which corresponds to a partial clustering ∂ C 1 ( v n ) , … , ∂ C K ( v n ). Consider the following problem, where v N is a leaf-node of the tree T ( N , K ), i.e. v N ∈ L ( N , K ): min v N ∈ L ( N , K ) ( ∑ k = 1 K Loss opt ( { p i , i ∈ C k ( v N ) } | R ∪ ) ) s.t. vertex v N is reached through v n . Looking at the enumeration tree in Fig. 1, in problem (7) we only consider clusterings which correspond to the branch starting at the vertex v n .

Consider the unique path v 1 v 2 … v n from v 1 to v n and label the corresponding partial partition with ∂ P n = ∂ P ( v n ). For short, we notate problem (7) with [ P 1 N / K ] | [ R ∪ , ∂ P n ], and the optimal solution with Loss opt ( [ P 1 N / K ] | [ R ∪ , ∂ P n ] ).

We now state the following theorem. (<italic>Optimal solution pruning criteria</italic>)<italic>.</italic>

Loss opt ( [ P 1 N / K ] | [ R ∪ , ∂ P n ] ) ⩾ Loss opt ( [ P 1 n / K ] | [ R ∪ , ∂ P n ] ) + Loss opt ( [ P n + 1 N / K ] | [ R ∪ ] ) .

The proof is rather analogous to Theorem 1. Suppose the optimal value for problem (7) is attained for clustering C 1 , … , C K with centres at locations μ 1 , … , μ K . Label C k † = ∂ C ( v n ) = C k ∩ N 1 n , C k ‡ = C k ∩ N n + 1 N , k = 1 , … , K . We have Loss opt ( [ P 1 N / K ] | [ R ∪ , ∂ P n ] ) = ∑ k = 1 K ∑ i ∈ C k d ( p i , μ k ) = ∑ k = 1 K ∑ i ∈ C k † d ( p i , μ k ) + ∑ k = 1 K ∑ i ∈ C k ‡ d ( p i , μ k ) . Now if we “relax” centre locations for the first term on the right-hand-side by allowing optimization, and, similarly, centre locations and partition for the second term, we can only obtain an improvement. These “relaxations” result precisely in the definitions of Loss opt ( [ P 1 n / K ] | [ R ∪ , ∂ P n ] ) and Loss opt ( [ P n + 1 N / K ] | [ R ∪ ] ).  □

An idea for pruning criteria. Suppose we have an upper bound L UB for the problem [ P 1 N / K ] | [ R ∪ ], i.e. suppose we know that Loss opt ( [ P 1 N / K ] | [ R ∪ ] ) ⩽ L UB . Such a bound is given by the optimal solution of any (fixed) particular clustering C 1 , C 2 , … , C K . A “good quality” lower bound can be obtained by a location-allocation algorithm as discussed in Section 2.2. Now, if we know that Loss opt ( [ P 1 n / K ] | [ R ∪ , ∂ P n ] ) + Loss opt ( [ P n + 1 N / K ] | [ R ∪ ] ) ⩾ L UB . Theorem 2 implies that Loss opt ( [ P 1 N / K ] | [ R ∪ , ∂ P n ] ) ⩾ L UB and the branch corresponding to ∂ P n can be cut, i.e. there is no need to explore it, as we will not find a better solution there.

 Bounds for Cluster Loss

Take some arbitrary cluster C k ⊆ { 1 , … , N } and consider the problem Loss opt ( { p i , i ∈ C k } | R ∪ ) = min μ ∈ R ∪ ∑ i ∈ C k d ( p i , μ ) , i.e. the problem of finding the optimal centre for a given set of points. For a distance metric d induced by a set of barriers, this is sufficiently complicated problem in itself (Klamroth, 2002). The book suggests solving (9) by formulating a mixed integer programming problem. Other suggested methods in the literature use heuristic algorithms (Aneja and Parlar, 1994; Bischoff and Klamroth, 2007).

Instead of trying to solve (9) directly, we will concentrate on finding a lower bound for the problem, i.e. we will define Loss LB ( { p i , i ∈ C k } | R ∪ ) (here, the superscript LB stands for the “lower bound”) which will satisfy Loss LB ( { p i , i ∈ C k } | R ∪ ) ⩽ Loss opt ( { p i , i ∈ C k } | R ∪ ) . The main thing we seek for this definition is that given the bounds for some cluster C k , it would be “easy” to find the bounds for cluster C k ∪ { j }, e.g. if we “increment” the cluster with “client” p j .

The Nearest and the Furthest Point in a Region-Unit

Remember that we have defined R ∪ as a union of a finite number of “region-units”: R ∪ = ⋃ m = 1 M R m . Consider an arbitrary “client” p n and an arbitrary “region-unit” R m ∈ { R 1 , … , R M }. Define the nearest and the furthest points μ n m inf and μ n m sup for p n in cell R m : μ n m inf = arg min μ ∈ R m d ( p n , μ ) , μ n m sup = arg max μ ∈ R m d ( p n , μ ) . The existence and finiteness of these points is ensured by the closed-and-bounded property of R m as discussed in Section 2.1. Also, we notate the minimum and maximum distances attained at these points: d n m inf = inf μ ∈ R m d ( p n , μ ) = d ( p n , μ n m inf ) , d n m sup = sup μ ∈ R m d ( p n , μ ) = d ( p n , μ n m sup ) .

The Lower and Upper Bounds for the Loss of a Cluster

Take some arbitrary cluster C k ⊆ { 1 , 2 , … , N } and some arbitrary region-unit R m ∈ { R 1 , … , R M }. Consider the following problem: Loss opt ( { p i , i ∈ C k } | R m ) = min μ ∈ R m ∑ i ∈ C k d ( p i , μ ) , the problem of finding the optimal “facility” location for cluster C k when the centre is restricted to region-unit R m .

The following bounds are valid: ∑ n ∈ C k d n m inf ⩽ Loss opt ( { p i , i ∈ C k } | R m ) ⩽ ∑ n ∈ C k d n m sup .

Define μ k opt [ m ] – the optimal centre position for cluster C k in region-unit R m , i.e. the solution of (13), which exists and is finite since R m is closed and bounded. We have Loss opt ( { p n , n ∈ C k } | R m ) = min μ ∈ R m ∑ n ∈ C k d ( p n , μ ) = ∑ n ∈ C k d ( p n , μ k opt [ m ] ) ⩾ ∑ n ∈ C k inf μ ∈ R m d ( p n , μ ) = ∑ n ∈ C k d n m inf . Analogously, Loss opt ( { p n , n ∈ C k } | R m ) = min μ ∈ R m ∑ n ∈ C k d ( p n , μ ) ⩽ ∑ n ∈ C k sup μ ∈ R m d ( p n , μ ) = ∑ n ∈ C k d n m sup .  □

Now consider the optimal loss for cluster C k in the “complete” constrained-set R ∪ , Loss opt ( { p i , i ∈ C k } | R ∪ ) – as defined in (9). Since R ∪ = ⋃ m = 1 M R m , it can be seen that Loss opt ( { p i , i ∈ C k } | R ∪ ) = min m = 1 , … , M Loss opt ( { p i , i ∈ C k } | R m ) .

For the optimal cluster loss Loss opt ( { p n , i ∈ C k } | R ∪ ), these bounds are valid: min m = 1 , … , M ∑ n ∈ C k d n m inf ⩽ Loss opt ( { p i , i ∈ C k } | R ∪ ) ⩽ min m = 1 , … , M ∑ n ∈ C k d n m sup .

Since the bounds ∑ n ∈ C k d n m inf ⩽ Loss opt ( { p i , i ∈ C k } | R m ) ⩽ ∑ n ∈ C k d n m sup are valid element-wise (see (14)), they must be also valid over the minimum values.  □

Here we note that our sub-task of finding bounds for C k ∪ { j } given bounds for C k becomes easy. For this, we have to “memorize” M sums ( ∑ n ∈ C k d n m inf ) and ( ∑ n ∈ C k d n m sup ), and using them we can obtain the updated bounds in O ( M ) time: Loss opt ( { p i , i ∈ C k ∪ { j } } | R ∪ ) ⩾ min m = 1 , … , M [ ( ∑ n ∈ C k d n m inf ) + d j m inf ] , Loss opt ( { p i , i ∈ C k ∪ { j } } | R ∪ ) ⩽ min m = 1 , … , M [ ( ∑ n ∈ C k d n m sup ) + d j m sup ] .

The Global Lower Bound Definition

Define Loss LB ( { p i , i ∈ C k } | R ∪ ) = min m = 1 , … , M ∑ n ∈ C k d n m inf . In Theorem 4 we have shown that this definition satisfies (10): Loss LB ( { p i , i ∈ C k } | R ∪ ) ⩽ Loss opt ( { p i , i ∈ C k } | R ∪ ) . Now take some arbitrary leaf-node v N ∈ L ( N , K ) and consider the related partition P N = { C 1 ( v N ) , … , C K ( v N ) } . For each cluster C k ( v N ), k = 1 , … , K, (19) holds. Summing over all the clusters, we obtain ∑ k = 1 K Loss LB ( { p i , i ∈ C k ( v N ) } | R ∪ ) ⩽ ∑ k = 1 K Loss opt ( { p i , i ∈ C k ( v N ) } | R ∪ ) . The right-hand-side of this inequality is the definition of Loss opt ( [ P 1 N / K ] | [ R ∪ , P N ] ). Notate the left-hand-side with Loss LB ( [ P 1 N / K ] | [ R ∪ , P N ] ) = ∑ k = 1 K Loss LB ( { p i , i ∈ C k ( v N ) } | R ∪ ) . Now define: Loss LB ( [ P 1 N / K ] | [ R ∪ ] ) = min v N ∈ L ( N , K ) Loss LB ( [ P 1 N / K ] | [ R ∪ , P ( v N ) ] ) . Since Loss opt ( [ P 1 N / K ] | [ R ∪ ] ) = min v N ∈ L ( N , K ) Loss opt ( [ P 1 N / K ] | [ R ∪ , P ( v N ) ] ) , we immediately see from definition (22) and inequality (20) that Loss LB ( [ P 1 N / K ] | [ R ∪ ] ) ⩽ Loss opt ( [ P 1 N / K ] | [ R ∪ ] ) .

Determination

Now consider the problem of finding Loss LB ( [ P 1 N / K ] | [ R ∪ ] ) as defined in (22). Suppose we want to find it using the enumeration tree T ( N , K ). We will need pruning criteria.

Take arbitrary v n ∈ T ( N , K ), the related partial partition ∂ P n = ∂ P ( v n ) and consider the following problem (in analogy with (7)): min v N ∈ L ( N , K ) ( ∑ k = 1 K Loss LB ( { p i , i ∈ C k ( v N ) } | R ∪ ) ) s.t. vertex v N is reached through v n . We label the solution to the above equation with Loss LB ( [ P 1 N / K ] | [ R ∪ , ∂ P n ] ). (<italic>Lower bound pruning criteria</italic>)<italic>.</italic>

Loss LB ( [ P 1 N / K ] | [ R ∪ , ∂ P n ] ) ⩾ Loss LB ( [ P 1 n / K ] | [ R ∪ , ∂ P n ] ) + Loss LB ( [ P n + 1 N / K ] | [ R ∪ ] ) .

Again, analogous to the proof of Theorem 1 and Theorem 2. Suppose the minimal value of problem (23) is attained for clustering C 1 , … , C K with optimal region-unit indexes m 1 , … , m K , e.g. Loss LB ( { p i , i ∈ C k } | R ∪ ) = min m = 1 , … , M ∑ n ∈ C k d n m inf = ∑ n ∈ C k d n m k inf , k = 1 , … , K . Define C k † = ∂ C k ( v n ) = C k ∩ N 1 n , C k ‡ = C k ∩ N n + 1 N , k = 1 , … , K . Now we enumerate the partitions of the index-set N n + 1 N . Consider the enumeration tree T ( N − n , K ) and its leaf-set L ( N − n , K ). Take any w N − n ∈ L ( N − n , K ) and consider the corresponding path w 1 w 2 … w N − n . Label C k ‡ ( w N − n ) = C k ‡ ( w 1 w 2 … w N − n ) = { i ∈ N n + 1 N : λ ( w i − n ) = k } . We have (in the equations below, d n m = d n m inf , as defined in (11a)) Loss LB ( [ P 1 N / K ] | [ R ∪ , ∂ P n ] ) = ∑ k = 1 K ∑ i ∈ C k d i m k = ∑ k = 1 K ∑ i ∈ C k † d i m k + ∑ k = 1 K ∑ i ∈ C k ‡ d i m k ⩾ ∑ k = 1 K ( min m = 1 , … , M ∑ i ∈ C k † d i m ) + min w ∈ L ( N − n , K ) ∑ k = 1 K ( min m = 1 … , M ∑ i ∈ C k ‡ ( w ) d i m ) = Loss LB ( [ P 1 n / K ] | [ R ∪ , ∂ P n ] ) + min w ∈ L ( N − n , K ) ∑ k = 1 K Loss LB ( { p i , i ∈ C k ‡ ( w ) } | R ∪ ) . The last term is the definition of Loss LB ( [ P n + 1 N / K ] | [ R ∪ ] ), see (21) and (22).  □

 MIP Formulation

Problem of finding Loss LB ( [ P n + 1 N / K ] | [ R ∪ ] ) can be seen as a discrete facility location problem. Indeed, we have to choose K “facilities” from { R 1 , … , R M } at positions m 1 , … , m K , so that the loss ∑ k = 1 K ∑ i ∈ C k d i m k inf is minimal. We formulate the corresponding mixed integer programming (MIP) problem. For this, introduce M binary variables { y 1 , … , y M }: y m = 1 , if “facility” at R m is opened, 0 otherwise. Also, introduce N × M binary variables { x n m , n = 1 , … , N , m = 1 , … , M }: x n m = 1 , if “client” p n is assigned to “facility” at R m , 0 otherwise. MIP formulation: min { y m } , { x n m } ∑ n = 1 N ∑ m = 1 M d n m inf x n m s.t. ∑ m = 1 M y m ⩽ K , x n m ⩽ y m , x n m ∈ { 0 , 1 } , y m ∈ { 0 , 1 } , n = 1 , … , N , m = 1 , … , M .

A Comment on the Error Between the Lower Bound and the Exact Solution

Define the “diameter” diam of an arbitrary region-unit R m as diam ( R m ) = sup x , y ∈ R m d ( x , y ) . In addition, suppose that distance metric d satisfies the triangle inequality: given any three points x , y , z ∈ S, it is always true that d ( x , z ) ⩽ d ( x , y ) + d ( y , z ) . In particular, this property holds for the distance metric d induced by a set of barriers, i.e. when d ( x , y ) is equal to the length of the shortest path between x and y (and the length of an arbitrary path from x to y is defined by the standard length of a curve in the Euclidean geometry).

Now take an arbitrary point p n ∈ { p 1 , p 2 , … , p N }. Consider the nearest and the furthest point for p n in region-unit R m , as was defined in (11). By the triangle inequality, we have d n m sup = d ( p n , μ n m sup ) ⩽ d ( p n , μ n m inf ) + d ( μ n m inf , μ n m sup ) ⩽ d n m inf + diam ( R m ) . Thus, if region-unit R m is sufficiently small ( diam ( R m ) ≈ 0), we can expect that d n m inf ≈ d n m sup . (If diam ( R m ) is not sufficiently small, we could divide R m into smaller parts and redefine R ∪ .)

Now consider some cluster C k and suppose that diam ( R m ) ⩽ ε for any R m ∈ { R 1 , R 2 , … , R M } for some ε > 0. By combining (16) and (25), Loss opt ( { p i , i ∈ C k } | R m ) ⩽ ∑ n ∈ C k d n m sup ⩽ (25) ∑ n ∈ C k ( d n m inf + diam ( R m ) ) ⩽ ∑ n ∈ C k d n m inf + | C k | ε , where | C k | is the number of items in cluster C k . Therefore, we have the following bounds: ∑ n ∈ C k d n m inf ⩽ Loss opt ( { p i , i ∈ C k } | R m ) ⩽ ∑ n ∈ C k d n m inf + | C k | ε . Similarly, for Loss opt ( { p i , i ∈ C k } | R ∪ ) we obtain min m = 1 , … , M ∑ n ∈ C k d n m inf ⩽ Loss opt ( { p i , i ∈ C k } | R ∪ ) ⩽ min m = 1 , … , M ∑ n ∈ C k d n m sup ⩽ min m = 1 , … , M ∑ n ∈ C k d n m inf + | C k | ε , and noting the definition of Loss LB ( { p i , i ∈ C k } | R ∪ ) in (18), inequalities (27) can be rewritten as Loss LB ( { p i , i ∈ C k } | R ∪ ) ⩽ Loss opt ( { p i , i ∈ C k } | R ∪ ) ⩽ Loss LB ( { p i , i ∈ C k } | R ∪ ) + | C k | ε . Summing over all the clusters of a particular partition P = { C 1 , … , C K }, we obtain (use the definitions in (6), (21) and inequality above) Loss LB ( [ P 1 N / K ] | [ R ∪ , P ] ) ⩽ Loss opt ( [ P 1 N / K ] | [ R ∪ , P ] ) ⩽ Loss LB ( [ P 1 N / K ] | [ R ∪ , P ] ) + N ε . In particular, the gap between the lower bound and the optimal solution for any partition P is at most N ε. Thus, we can control the level of this gap by subdividing R ∪ = ⋃ m = 1 M R m into sufficiently small region-units R m . These bounds also hold for the globally optimal solution: Loss LB ( [ P 1 N / K ] | [ R ∪ ] ) ⩽ Loss opt ( [ P 1 N / K ] | [ R ∪ ] ) ⩽ Loss LB ( [ P 1 N / K ] | [ R ∪ ] ) + N ε .

Summary

We have defined Loss LB ( { p i , i ∈ C k } | R ∪ ) in (18) and have proven in Theorem 4 that Loss LB ( { p i , i ∈ C k } | R ∪ ) ⩽ Loss opt ( { p i , i ∈ C k } | R ∪ ) . Also, we have defined Loss LB ( [ P 1 N / K ] | [ R ∪ ] ) in (22) which satisfies Loss LB ( [ P 1 N / K ] | [ R ∪ ] ) ⩽ Loss opt ( [ P 1 N / K ] | [ R ∪ ] ) . We have also proven in Theorem 5 the pruning criteria: Loss LB ( [ P 1 N / K ] | [ R ∪ , ∂ P n ] ) ⩾ Loss LB ( [ P 1 n / K ] | [ R ∪ , ∂ P n ] ) + Loss LB ( [ P n + 1 N / K ] | [ R ∪ ] ) , where ∂ P n = { ∂ C 1 , … , ∂ C K } is a (given) partition of index-set N 1 n .

 Overview of the Proposed Algorithm

Consider solving problem (1), i.e. clustering N points into K clusters. The algorithm we propose in this paper can be summarized with these steps:

Calculate d n m inf as defined in (11a). From programming point of view, this is the most difficult part of the algorithm which is problem specific. Fortunately, for multi-Weber problem with polyhedral barriers, calculation of d n m inf is possible (see Section 5.2).

Solve lower bound problems Loss LB ( [ P n + 1 N / K ] | [ R ∪ ] ) , for n = N − 1 , … , 0 , using the enumeration tree T ( N − n , K ) (illustrated in Fig. 1). Use pruning criteria (30): at step n, when looking for Loss LB ( [ P n + 1 N / K ] | [ R ∪ ] ), for some partial partition ∂ P of point set P n + 1 η , values Loss LB ( [ P n + 1 η / K ] | [ R ∪ , ∂ P ] ) , Loss LB ( [ P η + 1 N / K ] | [ R ∪ ] ) are known.

Suppose Loss LB ( [ P 1 N / K ] | [ R ∪ ] ) is attained for clustering C 1 LB , … , C K LB – as found in the previous step. For this clustering, for each cluster C k LB , k = 1 , … , K, solve the optimal centre problem Loss opt ( { p i , i ∈ C k LB } | R ∪ ) within ϵ-accuracy by subdividing the most promising region-units R m , m = 1 , … , M into smaller parts (e.g. as described in Section 5.3). This way, we obtain an upper bound L UB – as an optimal value of a particular solution corresponding to C 1 LB , … , C K LB .

Enumerate all v N ∈ L ( N , K ) s.t. Loss LB ( [ P 1 N / K ] | [ R ∪ , P ( v N ) ] ) ⩽ L UB . At any vertex v n ∈ T ( N , K ), prune the corresponding branch if Loss LB ( [ P 1 n / K ] | [ R ∪ , ∂ P ( v n ) ] ) + Loss LB ( [ P n + 1 N / K ] | [ R ∪ ] ) > L UB . Label the obtained set with V = { v N ∈ L ( N , K ) : v N satisfies (31) }.

Starting from the most promising v N ∈ V (e.g., with minimal values of Loss LB ( [ P 1 n / K ] | [ R ∪ , P ( v N ) ] )), solve each problem Loss opt ( [ P 1 N / K ] | [ R ∪ , P ( v N ) ] ) within ϵ-accuracy by partitioning region-units into smaller parts (e.g. as described in Section 5.3). Keep record of L UB – the upper bound of the best solution. If we get that Loss opt ( [ P 1 N / K ] | [ R ∪ , P ( v N ) ] ) > L UB , remove v N from V.

Numerical Example

We apply the ideas from previous sections to the multi-Weber with polyhedral barriers problem. An example of the problem and its ϵ-accuracy solution are illustrated in Fig. 2. The problem is defined within a unit-square with a 10 × 10 grid, with each square initially divided into 4 triangles. The blue polygons are the barriers; it is not allowed to travel through them – these barriers induce the distance metric d : S × S ↦ R + . There are 10 barriers consisting of 100 triangles in total. The region-union R ∪ consists of the remaining 300 gray-coloured triangles (region-units R m ): R ∪ = ⋃ m = 1 300 R m .

Location with barriers problem example. We seek to cluster 25 given points into 5 clusters. The solution is proved by our algorithm to be within ϵ = 0.01 accuracy to the global optimum (i.e. the globally optimal loss differs by at most 1 % when compared with the reported solution).

 Visibility Graph

Visibility polygons for selected points.

To compute the distances d ( p i , μ k ) between client-points and facility-points, we need to find the so-called visibility-graph G ( V , E ) (Viegas and Hansen, 1985). Each point p ∈ R ∪ defines a visibility-polygon as illustrated in Fig. 3. This polygon contains the set of points “seen” from the corresponding point directly through a ray; to reach the points outside this polygon, we must bypass some barriers. It is easy to convince oneself that (see Viegas and Hansen, 1985), in case point q ∈ R ∪ is not seen directly from point p, the shortest path d ( p , q ) consists of a set of segments with the end-points at the corners of the barrier-regions.6

Except the two segments “starting” at p and “ending” at q.

Therefore, the visibility-polygons for the barrier corner-points are particularly interesting (as in Fig. 3b). Barrier corner-points, together with client-points, define the set of vertices V for the visibility-graph G ( V , E ). The (directed) edges E are inserted as follows: if vertex v ∈ V is “seen” from a barrier-corner-vertex w (i.e. v belongs to the visibility polygon of w), we add the edge w ↦ v to graph (with edge weight equal to the Euclidean-length from point p ( v ) to point p ( w )).

To compute the shortest paths from a given facility-point μ to each client-vertex v ( p 1 ) , … , v ( p N ), we could insert an additional vertex v ( μ ) to the graph, find the set of vertices v ∈ V “seen” from μ, and insert the corresponding edges. Now, the shortest-paths can be found with Dijkstra’s algorithm (Dijkstra, 1959).

Visibility graph for our example problem is shown in Fig. 4.

Visibility graph for the example problem.

 Minimal Distances to a Triangle-Facility

To apply the ideas from the previous sections, we have to compute the lower-bound distances d n m inf , n ∈ N 1 N , m ∈ N 1 M , i.e. the shortest possible distance d n m inf from the client-point p n to the region-unit (triangle-facility in our case) R m .

For this task, we use the visibility-graph G ( V , E ) described in the previous subsection.

To find the lower-bound distances d n m inf , n ∈ N 1 N , we firstly insert an additional vertex v ( R m ) into the graph. Now, we have to take care of the edges. This can be done as follows: for each v ∈ V, intersect the visibility polygon VP ( v ) of v with triangle R m : I ( v , R m ) = VP ( v ) ∩ R m .7

Since the polygon VP ( v ) is non-convex, it might be that I ( v , R m ) is a not a single, but a set of polygons.

If the intersection I ( v , R m ) is empty, do nothing. Otherwise, find the nearest point proj ( v , R m ) from p ( v ) to I ( v , R m ).8

The nearest distance from a point to a polygon can be found by iterating the “projections” on the edges of the polygon. Iterating further over the polygons constituting I ( v , R m ), we obtain the result.

 Now, add to G ( V , E ) a directed-edge v ( R m ) ↦ v with weight ‖ p ( v ) − proj ( v , R m ) ‖ 2 . The distances d n m inf , n ∈ N 1 N for the triangle-facility R m can now again be found with Dijkstra’s algorithm (Dijkstra, 1959).

The output of this procedure is illustrated in Fig. 5.

To compute an upper bound for the sum of distances from triangle-facility R m to all the clients, we could pick any point μ ∈ R m and calculate the distances as outlined in Section 5.1, i.e. using the visibility graph. However, in our case where all the region-units are triangles, for finding an upper-bound faster, we can use the following simple idea.

Pick the longest edge of the triangle, compute its length l 1 .

Pick the middle point of this edge. Define l 2 as the distance from this middle-point to the vertex in front of the longest edge.

All the points within the triangle are at maximum distance l = max ( l 1 2 , l 2 ). The upper bound for the triangle-facility R m can be found by the formula ( ∑ n = 1 N d n m inf ) + N l.

Shortest distances to triangle-facilities. The triangle-facility in Fig. 5b has the minimal sum of distances ∑ n = 1 N d n m inf from all region-units R m , m = 1 , … , M.

 Finding the Optimal Centre Within a Predefined Accuracy <italic>ϵ</italic>

Illustration of the triangle-partitioning algorithm.

To find the globally-optimal centre for a set of points, i.e. optimal centree for some cluster C k , we use the triangle-partitioning algorithm in the style of Paulavičius and Žilinskas (2013). For a start, we place all the initial triangles (region-units R m ) in a MinPQ9

Minimum priority queue.

(Sedgewick and Wayne, 2011), ordered by the lower-bound value Loss LB ( { p i , i ∈ C k } | R m ) = ∑ i ∈ C k d i m inf , m = 1 , … , M . Initially, the MinPQ contains M elements. Now, pick from MinPQ the triangle-facility R min with the smallest lower-bound value (32). Delete this triangle from MinPQ. Next, partition R min into two smaller triangles R left and R right along the longest edge (see the partitions of the triangle-facility from Fig. 5b in Fig. 6a, Fig. 6b). Both these triangles have an improved (increased) lower-bound estimates Loss LB ( { p i , i ∈ C k } | R left ) , Loss LB ( { p i , i ∈ C k } | R right ) , because they were obtained by subdividing R min into two smaller pieces. Since the triangles have been shrunk, the gap between the lower bound and the upper bound for these triangles improves, i.e. we have shown that inequality (28) holds with ε left = diam ( R left ) and ε right = diam ( R right ). Insert R left and R right into MinPQ. Continue deleting/partitioning triangles in MinPQ until the distance between the best known lower-bound and the best known upper-bound is within a predefined accuracy ϵ, i.e. until condition Loss LB ( { p i , i ∈ C k } | R min ) Loss LB ( { p i , i ∈ C k } | R min ) + | C k | diam ( R min ) > 1 − ϵ holds for some “small” ϵ (take for example ϵ = 0.001). The procedure is illustrated in Fig. 6.

 Other Steps

Other steps in the algorithm of Section 4 given a routine to calculate d n m inf and a routine to find the optimal centre for a cluster (as described in Section 5.3) are rather straightforward. We just emphasize that during enumeration procedure at any vertex v n ∈ T ( N , K ), we use formula (17a) to find Loss LB ( { p i , i ∈ C λ ( v n ) ∪ { n } } | R ∪ ) in O ( M ) time.

 Algorithm Analysis Improvement by Initial Subdivision of Region-Units

Our initial analysis of the algorithm for the problem instance presented in Fig. 2 is summarized in Table 1.10

The numbers are particular for this problem instance, i.e. the one shown in Fig. 2 (if N < 25, we have used only points referring to 0 , … , N − 1). We did not seek to perform a statistical experiment for this part as we just wanted to illustrate that a denser subdivision into region-units significantly improves the algorithm.

In the experiment, we have analysed algorithm performance for the number of clients N = 20 , … , 25 and the number of clusters K = 3 , 4 , 5. Table 1a shows time statistics for lower bound calculations (Step 1 in the algorithm of Section 4), Table 1b shows time statistics for the enumeration of candidate solutions (Step 3 in the algorithm of Section 4), and Table 1c gives the statistics for the number of candidate solutions to the global optima.

Algorithm statistics for initial triangles (300 region-units). The question mark (?) indicates that the data is missing due to the time limit (see the last column in Table 1b and also compare with Table 2b).

While for the lower bound step times differ only slightly, the number of possible candidates to the global solution grows exponentially; the time to enumerate all these solutions is proportional. The reason to this is that our lower bound estimates are considerably (say, within 20 %-tolerance) smaller than the value of the globally optimal solution.

We try the following simple “remedy” to this issue: suppose that each triangle-unit R m , m = 1 , … , 300 is subdivided into 2 smaller triangles R m 1 , R m 2 .11

Giving R ∪ = ⋃ m = 1 300 ( R m 1 ∪ R m 2 ) – a union of 600 triangle-units.

Such a subdivision gives an improvement of the error between the lower bound and the exact solution (as discussed in Section 3.5). The updated statistics are presented in Table 2. Comparing with Table 1, we can see time increase for lower bound calculations (as expected), though this time is not crucial in the total time of the algorithm. But for the numbering the candidate solutions, the effect is clear: the number of candidates is reduced significantly.

Algorithm statistics when initial triangles are halved (600 region-units).

Continuing the experiment, we have halved the initial triangles 3 times, in total into 8 smaller triangles. The statistics are presented in Table 3. Comparing with Table 2, we can see that lower bound calculations take approximately 4 times more time, as expected – because we have 4 times more triangle-units. But as can be seen in Table 3c, the number of candidates to the global solution further decreases significantly.

Algorithm statistics when initial triangles are divided into 8 pieces (2400 region-units).

 Predicted time

Based on the observations of previous subsection, we decided to test our algorithm performance when initial triangle-units of R ∪ are subdivided 4 times into 16 smaller triangles. We have collected time statistics for various random problem instances. These problems were all generated within a unit-square with a 10 × 10 grid, with each square initially divided into 4 pieces (in total, 400 problem-triangles). 10 barriers were then randomly generated, with a total of 100 triangles classified as barrier-triangles. The remaining 300 were subdivided into 16 smaller triangles, resulting that R ∪ consists of M = 300 × 16 = 4800 triangle-units.

Models estimated using [R] quantreg package (Koenker, 2025).

Enumeration algorithm visited nodes log-proportion statistics. Number of clients N in the x-axis, number of clusters K in the title. Box-plots represent log-proportion statistics for a fixed N. Green line represents the linear model for the median log-proportion ( l 0.5 ( N , K ) = α 0.5 [ K ] + β 0.5 [ K ] N), red – the 90th percentile. 12

Algorithmic complexity of our approach is strongly related with the size of the tree T ( N , K ) of Fig. 1. Labeling with S 2 ( N , K ) the Stirling’s number of the second kind, and defining S 2 ( N , K ) = 0 when K > N, one can calculate that T ( N , K ) has N ( N , K ) = ∑ n = 1 N ∑ k = 1 K S 2 ( n , k ) nodes. Our algorithm performance depends on what proportion of the tree will be explored. During the experiment, we have recorded the number of tree-nodes visited in the enumeration of candidate solutions step (Step 3 in the algorithm of Section 4). The proportion statistics are represented in Fig. 7. Using the obtained linear models for log-proportion, we can predict the median time and the 90th percentile time outside of experiment data. On our machine,13

Processor 1.4 GHz Quad-Core Intel Core i5, Memory 8 GB 2133 MHz.

processing of a particular v n ∈ T ( N , K ) using formula (17a) takes approximately c = 9.098275 × 10 − 6 seconds.14

The algorithmic complexity of formula (17a) is roughly M additions and M conditionals; in the experiment, M is always fixed to 4800 triangle-units.

 Thus, we can predict algorithm time using the formulas below: T 0.5 ( N , K ) = c · exp { α 0.5 [ K ] + β 0.5 [ K ] N } · N ( N , K ) , T 0.9 ( N , K ) = c · exp { α 0.9 [ K ] + β 0.9 [ K ] N } · N ( N , K ) . In the above equations, N ( N , K ) is the number of nodes in the tree as defined in (33). The predicted times for various combinations of ( N , K ) when the number of triangle-units is fixed to M = 4800 are visualized in Fig. 8.

Estimated time of the algorithm for various problems [ P 1 N / K ] | [ R ∪ ] with K in the x-axis and N in the y-axis. Blue points represent the pairs for which time statistics were collected, red points have no time statistics corresponding to them and are predicted. In the figure, the letter next to a number means “s”-seconds, “m”-minutes, “h”-hours, “d”-days ( = 24 hours).

 Solved Example Instances

Some examples of problems solved within ϵ-accuracy are shown in Figs. 9 and 10. These are optimal solutions for largest problem instances we can solve in acceptable time. The results suggest themselves to be compared with the work of Rosing (1992), where the author had studied the unconstrained multi-Weber problem without barriers. This paper was published more than 30 years ago, but the results of their work are comparable to ours in the sense that, just like us (see Fig. 10), they reported global solutions for up to 35/5 and 30/6 problems. Our contribution is that we have managed to solve constrained multi-Weber problems with barriers of the same size, i.e. the ideas presented in this paper allow us to solve much more general class of uncapacitated multi-location problems.

Examples of within ϵ-accuracy solved problems ( ϵ = 1 %).

More examples of within ϵ-accuracy solved problems ( ϵ = 1 %).

 Conclusions

In this paper, we have outlined ideas how the globally-optimal solution of a very general multi-locations problem can be approached. The problem of finding the global solution perhaps is more of theoretical and intellectual than practical interest. Also, the question is more of proving that the best found solution is the global optimum. In practice, if one is faced with a barriers problem, a good-quality solution can be obtained by a location-allocation algorithm. The ideas of this paper can be easily extended for an implementation of such an algorithm.

We provided some numerical examples solving the multi-Weber with polyhedral barriers problem instances and predicted time for various sizes of problems. With our approach we have solved the problem instances of 60 / 2, 45 / 3, 40 / 4, 35 / 5 and 30 / 6 ( N / K means “N points into K clusters”) to global optimality (within ϵ-accuracy with ϵ = 1 %). To our knowledge, this is the first time the solution of instances with such a size to global optimality is reported.

 Acknowledgements

Our solver for the multi-locations with barriers problem was implemented with java programming language on top of the algorithms and data-structures from Sedgewick and Wayne (2011), which are available at https://algs4.cs.princeton.edu/home/. In particular, most of the figures were generated using the StdDraw.java class, and we have reused their implementations of Dijkstra’s algorithm (DijkstraSP.java) and various data-structures like MinPQ (MinPQ.java).

 References Aloise, D., Hansen, P., Liberti, L. (2012). An improved column generation algorithm for minimum sum-of-squares clustering. Mathematical Programming, 131, 195–220. Aneja, Y.P., Parlar, M. (1994). Algorithms for Weber facility location in the presence of forbidden regions and/or barriers to travel. Transportation Science, 28(1), 70–76. Bischoff, M., Klamroth, K. (2007). An efficient solution method for Weber problems with barriers based on genetic algorithms. European Journal of Operational Research, 177(1), 22–41. Bischoff, M., Fleischmann, T., Klamroth, K. (2009). The multi-facility location–allocation problem with polyhedral barriers. Computers & Operations Research, 36(5), 1376–1392. Cooper, L. (1964). Heuristic methods for location-allocation problems. SIAM Review, 6(1), 37–53. Dijkstra, E.W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, 1(1), 269–271. Drezner, Z., Klamroth, K., Schöbel, A., Wesolowsky, G. (2002). The Weber Problem. Facility Location: Applications and Theory, 1–36. du Merle, O., Villeneuve, D., Desrosiers, J., Hansen, P. (1999). Stabilized column generation. Discrete Mathematics, 194(1), 229–237. Gilmore, P.C., Gomory, R.E. (1961). A linear programming approach to the cutting-stock problem. Operations Research, 9(6), 849–859. Hamacher, H.W., Nickel, S. (1995). Restricted planar location problems and applications. Naval Research Logistics (NRL), 42(6), 967–992. Kepalas, M., Žilinskas, J. (2024). Solving net-constrained clustering problem. Journal of Nonlinear and Variational Analysis, 8, 987–1012. Klamroth, K. (2002). Single Facility Location Problems with Barriers. Springer, New York. Koenker, R. (2025). quantreg: Quantile Regression. Krau, S. (1997). Extensions du problème de Weber. PhD thesis, École Polytechnique de Montréal. Oğuz, M., Bektaş, T., Bennell, J.A., Fliege, J. (2016). A modelling framework for solving restricted planar location problems using phi-objects. The Journal of the Operational Research Society, 67(8), 1080–1096. Paulavičius, R., Žilinskas, J. (2013). Simplicial Global Optimization. Springer, New York. Rosing, K.E. (1992). An optimal method for solving the (generalized) multi-Weber problem. European Journal of Operational Research, 58(3), 414–426. Sedgewick, R., Wayne, K. (2011). Algorithms, 4th ed. Addison-Wesley Professional, Princeton. Sumner, M. (2025). sfdct: Constrained Triangulation for Simple Features. Viegas, J., Hansen, P. (1985). Finding shortest paths in the plane in the presence of barriers to travel (for any lp – norm). European Journal of Operational Research, 20(3), 373–381.