Revisiting several problems and algorithms in continuous location with lp norms

arXiv:1312.7473v1 [math.OC] 28 Dec 2013 Revisiting several problems and algorithms in continuous location with ℓpnorms V´ıctor Blanco, Justo Puerto, Safae El Haj Ben Ali Abstract This paper addresses the general continuous single facility location problems in finite dimension spaces under possibly different ℓpnorms in the demand points. We analyze the difficulty of this family of problems and revisit convergence properties of some well-known algorithms. The ultimate goal is to provide a common approach to solve the family of continuous ℓpordered median location problems in dimension d(including of course the ℓpminisum or Fermat-Weber location problem for any p≥1). We prove that this approach has a polynomial worse case complexity for monotone lambda weights and can be also applied to constrained and even non-convex problems. Keywords: Continuous location ·Ordered median problems ·Semidefinite programming ·Moment problem. 1 Introduction Location analysis is a very active topic within the Operations Research community. It has giving rise to a number of nowadays standard optimization problem some of them in the core of modern mathematical programming. One of its branches is continuous location a family of models directly related to important areas of mathematics such as linear and non-linear programming, convex analysis and global optimization (see e.g. [20] and the references therein). It is widely agreed that modern continuous location started with the paper by Weber (1909) [51] who first considers the minimization of weighted sums of distances as an economical goal to locate industries. This problem is currently known as Fermat-Weber, also because of the three points Fermat problem (s. XVII) firstly solved by Torricelli in 1659. The algorithmic part of this history starts at 1937 with the paper by Weiszfeld [52] who proposed an iterative gradient type algorithm to find or to approximate the solutions of the above mentioned Fermat-Weber problem. For several decades this algorithm remains forgotten but in 1973 Kuhn [29] rediscovered it and proved its convergence, under some conditions, in the Euclidean case. One year later Katz [30] gives another convergence result. Several years later, a number of authors considered the weighted minisum 1 problem under different norms mainly ℓpor polyhedral (see e.g. [20] for a detailed literature review) and Chandrasekaran and Tamir [12] raise several interesting questions concerning resolubility of Weiszfeld algorithm. Eventually, starting in the nineties, several authors were very interested in proving the convergence of some modifications of the Weiszfeld algorithm, usually called modified Weiszfeld or generalized iterative procedure for minisum location problems. The convergence for Euclidean distances (p=2) was studied later by [29, 30], among others. Since then, we can find in the literature many references concerning this algorithm, as for instance the generalization to ℓpdistances with p∈[1, 2][38] or the analysis of its local and global convergence [4, 5, 6, 10]. Also, these results were extended to more general problems: on Banach spaces [23, 43, 44], on the sphere [53], with regional demand [13, 49], with sets as demand facilities and using closest Euclidean distances [7] or with radial distances [14, 15, 16, 31, 39]. In addition, one can find in the literature papers where the convergence is accelerated using alternative step sizes [17, 18, 25, 30, 42] or some related properties concerning the termination of the algorithm in any of the demand points after a finite number of iterations [6, 9, 8, 11, 12, 29]. The influence of Weiszfeld algorithm in Location Analysis has been rather important very likely due to its very easy implementation. For several years, it was a very effective method to solve minisum continuous location problems, even though its theoretical convergence was not proven. Thus, locators have devoted a lot of effort to prove its convergence. The global convergence result of this algorithm for ℓpp∈[1, 2]was proved in [5] and recently [48] has given a proof to close the cases (p > 2) that were not yet justified. This has been an important effort from a mathematical point of view. Nevertheless, pursuing this goal locators did not focus on the origin of the problem, namely to search for alternative, efficient algorithms to solve the ℓpminisum and some more general families of location problems. The situation is even harder if we consider a more general family of location problems that have attracted a lot of attention in the field in the last years, namely continuous ordered median location problems [41]. Ordered median problems represent as special cases nearly all classical objective functions in location theory, including the Median, CentDian, center and k-centra. More precisely, the 1-facility ordered median problem can be formulated as follows: A vector of weights (λ1,...,λn)is given. The problem is to find a location for a facility that minimizes the weighted sum of distances where the distance to the closest point to the facility is multiplied by the weight λn, the distance to the second closest, by λn−1, and so on. The distance to the farthest point is multiplied by λ1. Many location problems can be formulated as the ordered 1-median problem by selecting appropriate weights. For example, the vector for which all λi=1is the unweighted 1-median problem, the problem where λ1=1and all others are 2 equal to zero is the 1-center problem, the problem where λ1=...=λk=1 and all others are equal to zero is the k-centrum. Minimizing the range of distances is achieved by λ1=1,λn= −1and all others are zero. Despite its full generality, the main drawback of this framework is the difficulty of solving the problems with a unified tool. There have been some successful approaches that are now available whenever the framework space is either discrete (see [3, 36]) or a network (see [27], [28], [40] or [45]). Nevertheless, the continuous case has been, so far, only partially covered even under the additional hypothesis of convexity. There have been some attempts to overcome this drawback and there are nowadays some available methodologies to tackle these problems, at least in the plane and with Euclidean norm. In Drezner [19] and Drezner and Nickel [21, 22] the authors present two different approaches. The first one uses a geometric branch and bound method based on triangulations (BTST) and the second one on a D-C decomposition for the objective function that allow solving problems on the plane. Espejo et al [24] also address the unconstrained convex ordered median location problem on the plane and Rodriguez-Chia et al. [47] attacks the k-centrum problem using geometric arguments and developing a better algorithm applicable only for that unconstrained problem on the plane and Euclidean distances. More recently, Blanco et al. [2] have presented a new methodology based on a hierarchy of SDP relaxations that can be used to solve (approximate) the optimal solutions of the general ordered median location problems which main drawback is the size of the SDP objects that have to be used to get good accuracy in high dimension. The above discussion points out that there exists a lack of a unified resolution approach to those problems as well as effective algorithms for the general cases. Our goal in this paper is to design a common approach to solve the family of continuous ℓpordered median location problems in dimension d(including of course the ℓpminisum or Fermat-Weber location problem for any p≥1). We prove that this approach has a polynomial worse case complexity for monotone lambda weights and can be also applied to constrained and to approximate even non-convex problems. Thus, providing a unifying new algorithmic paradigm for this class of location problems. First, for convex location problems it avoids the drawback of limit convergence proven for the Weiszfeld type algorithms. Then, it can be applied to any convex ordered median problem, even with mixed norms, in any dimension and with rather general convex constraints. Moreover, we show an explicit reformulation of these problems as SDP problems which enables the usage of standard free source solvers (SEDUMI, SDPT3,...) to solve them up to any degree of accuracy. Finally, we also show how to adapt this approach to approximate up to any degree of accuracy non-convex constrained location problems using a hierarchy of convergent relaxed problems. The paper is organized in 5 sections. In Section 2 we provide a compact representation, valid for any unconstrained convex ordered location prob- 3 lem, by means of a new formulation that reduces these problems to semidefinite problems. This approach allows us to ensure that all these problems are polynomially solvable in finite dimension. Section 3 is devoted to extend the results of Section 2 to the case of constrained problems under the condition of SDP-representability. Then, we handle the general case of nonconvex constrained ordered median location problem for which we construct a hierarchy of SDP relaxations that converges to the optimal solution of the original problem. Our Section 4 is devoted to the computational experiments. We report results in four different problems type, namely minisum (Weber), minimax (center), k-centrum (minimizing the sum of the k-largest distances) and general ordered median problems. The paper ends, in Section 5, with some conclusions and an outlook for further research . 2 A compact representation of the convex ordered median problem In this section we present the convex ordered median problem in dimension dwhere the distances are measured with a general ℓτ-norm being τ∈Q. We are given a set of demand points S={a1, ..., an}and two sets of scalars Ω:= {ω1, ..., ωn},ωi≥0, ∀i∈{1,...,n}and Λ:= {λ1, ..., λn}where λ1≥... ≥λn≥0. The elements ωiare weights corresponding to the importance given to the existing facilities ai, i ∈{1, ..., n}and depending on the choice of the elements of Λwe get different classes of problems. We denote by Pnthe set of permutations of the first nnatural numbers. Given a permutation σ∈ Pnsatisfying ωσ(1)kx−aσ(1)kτ≥...≥ωσ(n)kx−aσ(n)kτ, the unconstrained ordered median problem (see [41]) consists of min x∈Rd n X i=1 λiωσ(i)kx−aσ(i)kτ.(1) We start by showing a compact reformulation of the above problem that will be later useful in our approach. Special Theorem 1. Let τ=r sbe such that r, s ∈N\ {0},r > s and gcd(r, s) = 1. For any set of lambda weights satisfying λ1≥... ≥λn, 4 Problem (1) is equivalent to min n X k=1 vk+ n X i=1 wi(2) s.t vi+wk≥λkzi,∀i, k =1, ..., n, (3) yij −xj+aij ≥0, i =1, . . . , n, j =1, ..., d, (4) yij +xj−aij ≥0, i =1, . . . , n, j =1, ..., d, (5) yr ij ≤us ijzr−s i, i =1, . . . , n, j =1, ..., d, , (6) ω r s iPd j=1uij ≤zi, i =1, . . . , n, (7) uij ≥0, i =1, . . . , n, j =1, . . . , d. (8) Proof. Because of the condition λ1≥... ≥λn, Problem (1) can be equivalently written as min x∈Rdmax σ∈Pn n X i=1 λiωσ(i)kx−aσ(i)kτ,(9) Let us introduce auxiliary variables zi,i=1,...,n to which we impose that zi≥ωikx−aikτ, to model the problem in a convenient form. Now, for any permutation σ∈ Pn, let zσ= (zσ(1),...,zσ(n)). Moreover, let us denote by (·)the permutation that sorts any vector in nonincreasing sequence, i.e. z(1)≥z(2)≥... ≥z(n). Using that λ1≥... ≥λnand since zi≥0, for all i=1,...,n then n X i=1 λiz(i)=max σ∈Pn n X i=1 λizσ(i). The permutations in Pncan be represented by the following binary variables pik =1, if zigoes in position k, 0, otherwise, imposing that they verify the following constraints:            n X i=1 pik =1, ∀k=1, ..., n, n X k=1 pik =1, ∀i=1, ..., n. (10) Next, combining the two sets of variables we obtain that the objective 5 function of (9) can be equivalently written as::                          n X i=1 λiz(i)=max n X i=1 n X k=1 λkzipik s.t n X i=1 pik =1, ∀k=1, ..., n, n X k=1 pik =1, ∀i=1, ..., n, pik ∈{0, 1}. (11) Now, we point out that for fixed z1, ..., zn, the above problem is an assignment problem and its constraint matrix is totally unimodular, so that solving a continuous relaxation of the problem always yields an integral solution vector [1], and thus a valid permutation. Moreover, the dual of the linear programming relaxation of (11) is strong and also gives the value of the original binary formulation of (11). Hence, for any vector z∈Rn, by using the dual of the assignment problem (11) we obtain the following expression      n X i=1 λiz(i)=min n X k=1 vk+ n X i=1 wi s.t vi+wk≥λkzi,∀i, k =1, ..., n. (12) Finally, we replace (12) in (9) and we get          min n X k=1 vk+ n X i=1 wi s.t vi+wk≥λkzi,∀i, k =1, ..., n, zi≥ωikx−aikτ, i =1, ..., n. (13) It remains to prove that each inequality zi≥ωikx−aikτ, i =1, ..., n can be replaced by the system: yij −xj+aij ≥0, j =1, ..., d. yij +xj−aij ≥0, j =1, ..., d. yr ij ≤us ijzr−s i, j =1, ..., d. ω r s i d X j=1 uij ≤zi, uij ≥0, ∀j=1, . . . , d. Indeed, set ρ=r r−s, then 1 ρ+s r=1. Let (¯ x, ¯ zi)fulfill the inequality 6 zi≥ωikx−aikτ. Then we have ωik¯ x−aikτ≤¯ zi⇐⇒ ωi  d X j=1 |¯ xj−aij|r s  s r ≤¯ z s r i¯ z 1 ρ i ⇐⇒ ωi  d X j=1 |¯ xj−aij|r s¯ z r s(− r−s r) i   s r ≤¯ z s r i ⇐⇒ ω r s i d X j=1 |¯ xj−aij|r s¯ z−r−s s i≤¯ zi(14) Then (14) holds if and only if ∃ui∈Rd,uij ≥0, ∀j=1, ..., d such that |¯ xj−aij|r s¯ z−r−s s i≤uij,satisfying ω r s i d X j=1 uij ≤¯ zi, or equivalently, |¯ xj−aij|r≤us ij¯ zr−s i, ω r s i d X j=1 uij ≤¯ zi.(15) Set ¯ yij =|¯ xj−aij|and ¯ uij =|¯ xj−aij|τ¯ z−1/ρ i. Then, clearly (¯ x, ¯ zi,¯ y, ¯ u)satisfies (4)-(8). Conversely, let (¯ x, ¯ zi,¯ y, ¯ u)be a feasible solution of (4)-(8). Then, ¯ yij ≥ |¯ xij −aij|for all i, j and by (6) ¯ uij ≥¯ y(r/s) ij z−r−s s i≥|¯ xj−aij|τ¯ z−r−s s i. Thus, ωr/s i d X j=1 |¯ xj−aij|r/s¯ z−r−s s i≤ωr/s i d X j=1 ¯ uij ≤¯ zi, which in turns implies that ωr/s iPd j=1|¯ xj−aij|r/s ≤¯ zi¯ z r−s s iand hence, ωik¯ x− aik ≤ ¯ zi.  Problem (2) -(8) is an exact representation of Problem (1) in any dimension and for any ℓτ-norm. For the case τ=1, the above problem reduces to a linear programming problem. 7 Corollary 2. If τ=1the reformulation given by Problem (1) is                  min n X k=1 vk+ n X i=1 wi s.t vi+wk≥λkzi,∀i, k =1, ..., n, zi≥ωiPd j=1uij, i =1, ..., n, xj−aij ≤uij i=1, ..., n, j =1, . . . , d, −xj+aij ≤uij i=1, ..., n, j =1, . . . , d. (16) The reader may observe that the representation given in Theorem 1 is new and different from the one used in [2]. On the one hand, this new formulation is more efficient than the one presented in [2] and specially tailored for the case of non-increasing monotone lambda weights, see e.g. [2, Lemma 8]. For the sake of readability we include it in the following. Let Sk(x) := k X j=1 z(j),(17) where z(j)is such that z(1)≥... ≥z(n). That formulation applied to the setting of this paper reads as: min Pn k=1(λk−λk+1)Sk(x)(18) tk+rkj ≥zj(x), j, k =1, . . . , n, rkj ≥0, j, k =1, . . . , n, zj≥ωjkx−ajkτ, j =1, ..., n. It is easy to see that formulation (18) has O(n2+d)variables and O(n2) constraints whereas the new one written in similar terms as presented in (13) has O(n+d)variables and O(n2)constraints. Our goal is to show that for any τ∈Q, Problem (2) -(8) also admits a compact formulation within another easy class of polynomially solvable mathematical programming problems. In order to get that we need to prove a technical lemma. Let #Adenote the cardinality of the set A. Lemma 3. Let τ=r sbe such that r, s ∈N\ {0},r > s and gcd(r, s) = 1. Let x, u and tbe non negative and satisfying xr≤ustr−s.(19) Assume that 2k−1< r ≤2kwhere k∈N\ {0}such that x2k≤ustr−sx2k−r,(20) 8 and s=αk−12k−1+αk−22k−2+...+α121+α020,(21) r−s=βk−12k−1+βk−22k−2+...+β121+β020,(22) 2k−r=γk−12k−1+γk−22k−2+...+γ121+γ020,(23) where αi, βi, γi∈{0, 1}. Then, if (x, t, u)is a feasible solution of (19) there exists wsuch that either 1. (x, t, u, w)is a solution of System (24), if αi+βi+γi=1, for all 0<i<k−1.   w2 1≤uα0tβ0xγ0, w2 i+1≤wiuαitβixγi, i =1,...,k−2 x2≤wk−1d, (24) where d=wk−2if αk−1βk−1+γk−1=0 uαk−1tβk−1uγk−1if αk−1βk−1+γk−1=1 2. (x, t, u, w)is a solution of System (25), if there exist ijand il(j),j= 1,...,c such that: 1. 0 < i1< i2< . . . < ic≤k−2, 2. ij< il(j)< ij+1, 3. αij+βij+γij=3,αil(j)+βil(j)+γil(j)=0and αh+βh+γh=2 for h=ij+1,...,il(j)−1.                                            w2 1≤uα0tβ0xγ0, w2 i+1≤wiuαitβixγi, i ∈{1,...,i1−1} −−−−−−−−for each j=1,...,c−−−−−−−−− w2 θ(j)≤ut, w2 θ(j)+1≤wθ(j)−1x w2 θ(j)+2∗s≤wθ(j)+2(s−1)aij+s w2 θ(j)+2∗s+1≤wθ(j)+2s−1bij+s,s=1,...,il(j)−ij−1 and aij+s+bij+s=uαij+stβij+sxγij+s, w2 θ(j)+2(il(j)−ij)≤wθ(j)+2(il(j)−ij−1)wθ(j)+2(il(j)−ij−1)+1,if m −1 > θ(j) + 2(il(j)−ij−1), w2 θ(j)+2(il(j)−ij)+s≤wθ(j)+2(il(j)−ij)+s−1uαil(j)+stβil(j)+sxγil(j)+s,for all s=1, . . . , ij+1−il(j)−1 −−−−−−−−−−−−−−−−−−−−−−−−−−−−−− x2≤wmd. (25) where θ= (θ(j))c j=1such that θ(j) = 2#{i:αi+βi+γi≥2, 1 < i ≤ ij}+#{i:αi+βi+γi≤1, 1 < i ≤ij}for j=1, ..., c,m=1+2#{i: αi+βi+γi≥2, 1 < i < k}+#{i:αi+βi+γi=1, 1 < i < k}≤2k and d=wm−1if αk−1+βk−1+γk−1=0 uαk−1tβk−1uγk−1if αk−1+βk−1+γk−1=1. 9 form a2≤bc, involving 3 variables, a, b, c with b, c non negative. Finally, it is well-known by Schur complement that a2≤bc ⇔  b+c 0 2a 0 b +c b −c 2a b −c b +c  0, b +c≥0. Hence, Problem (1) is a SDP because it has a linear objective function, n2+ n(2d+1)linear inequalities and at most 2nd log(r)linear matrix inequalities.  There is an interesting observation that follows from the above result. It was already known that continuous convex ordered location problems with ℓ1norm were reducible to linear programming (see e.g. [41]). This paper proves that most continuous convex ordered location problems with ℓpnorms are reducible to SDP programming showing the similarities existing between all this class of problems and moreover that convex continuous single facility location problems are among the “easy” optimization problems. Moreover, Theorem 6 allows us to apply the general theory of SDP to derive a general result of convergence for solving the family of continuous convex ordered single facility location problems: Problem (1) is polynomially solvable in dimension d∈Nand for any set of nonincreasing lambda weights. Moreover, we can be more precise and can state the following result: Special Theorem 7. Let ε > 0 be a prespecified accuracy and (X0, S0) be a feasible primal-dual pair of initial solutions of Problem (30-35). An optimal primal-dual pair (X, S)satisfying X·S≤εcan be obtained in at most O(αlog X0·S0 ε)iterations and the complexity of each iteration is bounded above by O(αβ3, α2β2, α2)being α=3n +2nd(1+log r)and β=p, the dimension of the dual matrix variable Sp. The reader may observe that this result is mainly of theoretical interest because the bound is based on general results on primal-dual algorithms, such as the modification of Kouleai and Terlaki [33] to the Mehrota type algorithm [37] applied to SDP problems. Nevertheless, it is important to realize that it states an important difference with respect to any other known result in the area of continuous location where convergence results, when available, are only proven limit of sequences and never in finite number of steps nor accuracy ensured. In this case, one can ensure a prespecified accuracy of the solution in a known number of iterations. 3 Constrained Ordered median problems This section extends the above results to constrained location problems. Therefore, we address now the restricted case of Problem (1). Let {g1,...,gl}⊂ 16 R[x]be real polynomials and K:= {x∈Rd:gj(x)≥0, j =1,...,l} a basic closed, compact semialgebraic set with nonempty interior satisfying the Archimedean property. Recall that the Archimedean property is equivalent to impose that for some M > 0 the quadratic polynomial u(x) = M−Pd i=1x2 ihas a representation on Kas u=σ0+Pℓ j=1σjgj, for some {σ0,...,σl}⊂R[x]being each σjsum of squares [46]. We remark that the assumption on the Archimedean property is not restrictive at all, since any semialgebraic set K⊆Rdfor which is known that Pd i=1x2 i≤M holds for some M > 0 and for all x∈K, admits a new representation K′=K∪{x∈Rd:gl+1(x) := M−Pd i=1x2 i≥0}that trivially verifies the Archimedean property. In our framework the compactness assumption which is usually assumed in location analysis implies that this condition always holds. In this framework we assume that the domain Kis compact and has nonempty interior. We observe that we can extend the results in Section 2 to a broader class of convex constrained problems. In order to do that we need to introduce some notation. Let κ= (κα)be a real sequence indexed in the monomial basis (xβzγvδwζuηyθ)of R[x, z, v, w, u, y](with α= (β, γ, δ, ζ, η, θ)∈Nd×Nn×Nn×Nn×Nn×d× Nn×d). Let D=3n + (2n +1)ddenote the dimension of the space of variables. Define Υ= (x, z, v, w, u, y)to be the vector of indeterminates so that Υα=xβzγvδwζuηyθ. For any integer Nconsider the monomial vector [ΥN] = [(x, z, v, w, u, y)N] = [1, x1,...xd, z1,...zn,...,ynd, x2 1, x1x2,...,yN nd]t. Then, [ΥN][ΥN]tis a square matrix and we write [ΥN][ΥN]t=X 0≤|α|≤2N A0 αΥα for some symmetric 0/1 matrices A0 α. Here, for a vector α,|α|stands for the sum of its components. For any sequence, κ= (κα)α∈ND⊂R, indexed in the canonical monomial basis B, let Lκ:R[Υ]→Rbe the linear functional defined, for any f= Pα∈NdfαΥα∈R[Υ], as Lκ(f) := Pαfακα. The moment matrix MN(κ)of order Nassociated with κ, has its rows and columns indexed by (Υα)and MN(κ)(α, α′) := Lκ(Υα+α′) = κα+α′, for |α|,|α′|≤N. Therefore, MN(κ) = X 0≤|α|≤2N A0 ακα Note that the moment matrix of order Nhas dimension D+N D×D+N D and that there are D+2N Dκαvariables. 17 For g∈R[Υ] (= Pν∈NMgνΥν), the localizing matrix MN(gκ)of order Nassociated with κand g, has its rows and columns indexed by (Υα)and MN(gκ)(α, α′) := Lκ(Υα+α′g(Υ)) = Pνgνκν+α+α′, for |α|,|α′|≤N. Therefore, MN(gκ) = X 0≤|α|≤2N Ag ακα, for some symmetric 0/1 matrices Ag αthat depend on the polynomial g. Also for convenience, we shall denote by Ag eithe matrix associated with κeithe moment variable linked to the monomial xei=x1 i. (The interested reader is referred to [34] and [35] for further details on the moment approach applied to global optimization.) Special Theorem 8. Consider the restricted problem: min x∈K⊂Rd n X i=1 λiωσ(i)kx−aσ(i)kτ.(36) Assume that the hypothesis of Theorem 6 holds. In addition, any of the following conditions hold: 1. gi(x)are concave for i=1,...,ℓ and −Pℓ i=1µi∇2gi(x)≻0for each dual pair (x, µ)of the problem of minimizing any linear functional ctx on K(Positive Definite Lagrange Hessian (PDLH)). 2. gi(x)are sos-concave on Kfor i=1,...,ℓ or gi(x)are concave on Kand strictly concave on the boundary of Kwhere they vanish, i.e. ∂K∩∂{x∈Rd:gi(x) = 0}, for all i=1,...,ℓ. 3. gi(x)are strictly quasi-concave on Kfor i=1,...,ℓ. Then, there exists a constructive finite dimension embedding, which only depends on τand gi,i=1,...,ℓ, such that (36) is a semidefinite problem. Proof. The unconstrained version of Problem (36) can be equivalently written as a SDP using the result in Theorem 6. Therefore, it remains to prove that under the conditions 1, 2 or 3 the constraint set x∈Kis also exactly represented as a finite number of semidefinite constraints or equivalently that it is semidefinite representable (SDr). Let us begin with condition 1. Consider the system of linear matrix inequalities: A(k) 0+ ℓ X i=1 Agk eixi+X 1≤α≤2N Agk ακα0, k =0, . . . , ℓ. (37) Under the hypothesis of condition 1, the set Ksatisfies the Putinar-Prestel’s Bounded Degree Nonnegative Representation property (PP-BDNR), see [26, 18 Theorem 6]. This condition ensures that there exists a finite Nsuch that the set ^ SN={(x, κ) : satisfying inequalites (37)} projects via the xcoordinate onto the set K. Hence, an exact lifted representation of Problem (36) is the one provided by Theorem 6 augmented with the additional linear matrix inequalities in (37). Let us assume now that condition 2 holds. Consider the set ^ SN={(x, κ) : MN(κ)0, Lκ(gi)≥0, i =1 . . . , ℓ, Lκ(xj) = xj, j =1, . . . , d, κ0=1}. Under our hypothesis, Theorem 11.11 in [35] ensures that there exists a finite Nsuch that ^ SNprojects via the xvariables onto the set K. Hence, we obtain another lifted exact SDP formulation for Problem (36) using the formulation induced by Theorem 6 augmented with the inequalities MN(κ) 0, Lκ(gi)≥0, i =1 . . . , ℓ, Lκ(xj) = xj, j =1, . . . , d, κ0=1. Finally, let us consider the case in condition 3. If giare strictly quasiconcave on K, Proposition 10 in [26] implies that one can find some new polynomials −pithat have positive definite Hessian in K. Let us denote P:= {x∈Rd:pi(x)≥0, i =1,...,ℓ}. Thus, in some open set Ucontaining Kit holds P∩U=K. Next, define the set ^ SN={(x, κ) : satisfying inequalites (38)-(41)} where the set of linear matrix inequalities (38)-(41) are given by: A(k) 0+ ℓ X i=1 A(k) eixi+X 1≤α≤2N Apk ακα0, k =0,...,ℓ (38) Lκ(pk)≥0, k =0,...,ℓ (39) Lκ(xj) = xj, j =1,...,d (40) κ0=1. (41) Under the hypothesis of condition 3, Theorem 24 in [26] ensures that there exists a finite Nsuch that ^ SNprojects via the xvariables onto K. Hence, we obtain the third lifted exact SDP formulation for Problem (36) using the formulation induced by Theorem 6 augmented with the inequalities (38)-(41). We observe that according to Theorem 29 in [26], since we assume the Archimedean property holds in all these cases, Ncan be bounded above by some finite constant that only depends on the polynomials gi,i=1,...,ℓ.  19 We shall finish this section with another convergence result applicable to the case of non-convex constrained location problems. Again, let {g1,...,gl}⊂R[x]and K:= {x∈Rd:gk(x)≥0, k =1,...,ℓ}a basic, compact, closed semialgebraic set satisfying the Archimedean property, with nonempty interior and such that Kdoes not satisfies the hypothesis of Theorem 8, in particular some of the gjmay not be concave. Now, we can prove a convergence result that allows us to approximate, up to any degree of accuracy, the solution of the class of problems defined in (36) when the hypothesis of Theorem 8 fails. Let ξk:= ⌈(deg gk)/2⌉ where {g1,...,gℓ}are the polynomial constraints that define K. For N≥ N0:= max{max k=1,...,ℓ ξk, 1}, we introduce the following hierarchy of semidefinite programs: (QN) : min n X k=1 vk+ n X i=1 wi(42) s.t. vi+wk≥λkzi,∀i, k =1, ..., n, (43) yij −xj+aij ≥0, ∀i=1, ..., n, j =1, ..., d. (44) yij +xj−aij ≥0, ∀i=1, ..., n, j =1, ..., d. yr ij ≤us ijzr−s i,∀i=1, ..., n, j =1, ..., d, (45) ω r s iPd j=1uij ≤zi,∀i=1, ..., n, (46) MN(κ)0, (47) MN−ξk(gk, κ)0, k =1, . . . , ℓ, (48) Lκ(xj) = xj, j =1, . . . , d, Lκ(zi) = zi, i =1, . . . , n, Lκ(vi) = vi, i =1, . . . , n, Lκ(wi) = wi, i =1, . . . , n, Lκ(uij) = uij, i =1, . . . , n, j =1, . . . , d, Lκ(yij) = yij, i =1, . . . , n, j =1, . . . , d, κ0=1 uij ≥0, ∀i=1, ..., n, j =1, . . . , d. (49) with optimal value denoted min QN. Special Theorem 9. Consider ρλdefined as the optimal value of the problem: ρλ=min x∈K⊂Rd n X i=1 λiωσ(i)kx−aσ(i)kτ.(50) Then, with the notation above: (a) min QN↑ρλas N→∞. 20 (b) Let κNbe an optimal solution of Problem (QN). If rank MN(κN) = rank MN−N0(κr) = ϑ then min QN=ρλand one may extract ϑpoints (x∗(i), z∗(i), v∗(i), w∗(i), u∗(i), y∗(i))ϑ i=1⊂K, all global minimizers of Problem (50). Proof. First of all, we observe that an optimal solution of Problem (50) does exist by the compactness assumption on K. Moreover, the convergence of the semidefinite sequence of problems (QN)follows from a result by Lasserre [35, Theorem 5.6] that it is applied here to the SDP problem (30-35) on the closed semialgebraic set K. The second assertion on the rank condition, for extracting optimal solutions, follows from applying [35, Theorem 5.7] to Problem (QN).  4 Computational Experiments A series of computational experiments have been performed in order to evaluate the behavior of the proposed methodology. Programs have been coded in MATLAB R2010b and executed in a PC with an Intel Core i7 processor at 2x 2.93 GHz and 16 GB of RAM. The semidefinite programs have been solved by calling SDPT3 4.0[32]. We run the algorithm for several well-known continuous single facility convex ordered location problems: Weber, center, k-center and general ordered median problem with random non-increasing monotone lambda. For each of them, we obtain the CPU times for computing solutions as well as the accuracy given by the solver SDPT3 4.0. In addition, to illustrate the application of the result in Theorem 9, we also report results on a problem which consists of minimizing the range of distances in R3with two additional non-convex constraints. In this case, we include running times and gap with respect to upper bounds obtained with the battery of functions in optimset of MATLAB which only provide approximations on the exact solutions (optimality cannot be certified). In this last case, in order to compute the accuracy of an obtained solution, we use the following measure for the error (see [50]): ǫobj =|the optimal value of the SDP −fopt| max{1, fopt},(51) 21 Dim 2 3 10 τn Time(Ave) Gap(Ave) Time(Ave) Gap(Ave) Time(Ave) Gap(Ave) τ=1.5 10 0.20 0.00000001 0.28 0.00000001 0.86 0.00000000 100 1.71 0.00000001 3.16 0.00000001 10.89 0.00000001 500 10.78 0.00000001 15.84 0.00000001 51.23 0.00000000 1000 21.22 0.00000001 30.67 0.00000001 103.17 0.00000001 5000 103.50 0.00000001 178.50 0.00000001 566.64 0.00000001 10000 210.22 0.00000001 455.05 0.00000001 1330.36 0.00000001 τ=2 10 0.06 0.00000000 0.12 0.00000000 0.46 0.00000000 100 0.40 0.00000000 0.99 0.00000000 4.63 0.00000000 500 1.50 0.00000000 5.83 0.00000000 21.63 0.00000000 1000 3.27 0.00000000 11.08 0.00000000 44.04 0.00000000 5000 17.06 0.00000000 58.16 0.00000000 218.82 0.00000000 10000 33.19 0.00000000 118.91 0.00000000 455.30 0.00000000 τ=3 10 0.19 0.00000001 0.31 0.00000001 0.99 0.00000001 100 1.88 0.00000001 3.60 0.00000001 12.71 0.00000001 500 10.82 0.00000001 17.87 0.00000001 57.91 0.00000001 1000 21.73 0.00000001 33.99 0.00000001 118.11 0.00000001 5000 110.87 0.00000001 181.17 0.00000001 646.46 0.00000001 10000 245.66 0.00000001 477.38 0.00000001 1616.26 0.00000001 τ=3.5 10 0.33 0.00000001 0.47 0.00000001 1.75 0.00000001 100 3.87 0.00000001 5.44 0.00000001 19.71 0.00000001 500 18.62 0.00000001 26.99 0.00000001 92.64 0.00000001 1000 37.06 0.00000001 51.50 0.00000001 192.77 0.00000001 5000 280.27 0.00000001 304.94 0.00000001 1178.17 0.00000001 10000 964.18 0.00000001 872.29 0.00000001 2431.58 0.00000001 Table 1: Computational results for Weber problem with different norms and different dimensions. where fopt is the approximated optimal value obtained with the functions in optimset. The reader may note that we solve relaxed problems that give lower bounds. Therefore, the gap of our lower bounds is computed with respect to upper bounds which implies that actually may be even better than the one reported. (See Table 5.) We have organized our computational experiments in three different problems types. Our test problems are set of points randomly generated on the [0, 10000]hypercubes of the d-dimensional space, d=2, 3, 10. For Weber, center and k-centrum problems, we could solve instances with at least 10000 points and for different ℓτ-norms, τ=1.5, 2, 3, 3.5. The general case with random lambda weights is harder and we only solved in all cases instances up to 1000 points. 22 Our goal is to present the results organized per problem type, framework space (Rd,d=2, 3, 10) and norm (ℓτ,τ=1.5, 2, 3, 3.5). Tables 1, 2 and 3 report our results on the problems of minimizing the weighted sum of distances (Weber), the maximum distance (center) and the sum of the n/2 largest distances (n/2-centrum). In all cases, the accuracy and resolutions times needed for the solver are rather good, even for 10000 points in dimension 10 and rather complicated norms (e.g. ℓ3.5). The reader can see that the hardest type is the k-centrum. In this problem type CPU times increase one order of magnitude because the structure of the problem does not allow to reduce the size of the formulation. Dim 2 3 10 τn Time(Ave) Gap(Ave) Time(Ave) Gap(Ave) Time(Ave) Gap(Ave) τ=1.5 10 0.22 0.00000103 0.39 0.00000000 1.02 0.00010830 100 2.13 0.00000000 8.75 0.00000000 43.28 0.00002092 500 12.59 0.00000000 63.16 0.00000000 237.99 0.00000835 1000 27.11 0.00000000 115.35 0.00000000 327.04 0.00000615 5000 150.34 0.00000001 357.49 0.00000000 1231.78 0.00000539 10000 371.39 0.00000000 1297.23 0.00000000 2762.51 0.00000714 τ=2 10 0.11 0.00000005 0.17 0.00000001 0.44 0.00002446 100 1.34 0.00000012 2.09 0.00000001 7.68 0.00024196 500 8.80 0.00000008 13.69 0.00000001 50.25 0.00001027 1000 20.85 0.00000051 30.27 0.00000001 115.64 0.00001137 5000 119.90 0.00000133 212.21 0.00000001 912.44 0.00384177 10000 287.13 0.00000265 467.08 0.00000001 1510.41 0.00823965 τ=3 10 0.23 0.00000021 0.37 0.00000001 1.08 0.00010908 100 2.32 0.00000000 7.48 0.00000000 37.66 0.00001214 500 14.47 0.00000001 52.27 0.00000018 209.46 0.00000624 1000 28.93 0.00000001 119.96 0.00000056 293.48 0.00002443 5000 160.96 0.00000001 456.11 0.00010340 1223.41 0.00024982 10000 434.79 0.00000000 6829.70 0.00002703 2663.15 0.00026867 τ=3.5 10 0.33 0.00000000 0.56 0.00000570 1.80 0.00004732 100 4.47 0.00000013 13.72 0.00000036 68.91 0.00002253 500 21.93 0.00000002 90.65 0.00000020 373.21 0.00000236 1000 44.82 0.00000002 179.48 0.00001955 603.41 0.00001127 5000 244.04 0.00000002 551.16 0.00000219 2279.93 0.00020991 10000 510.25 0.00000002 2618.87 0.00001699 4814.26 0.00006055 Table 2: Computational results for the Center problem with different norms and different dimensions. Table 4 reports our results on the general ordered median problem with non-increasing monotone lambda weights. For this family of problems we 23 Dim 2 3 10 τn Time(Ave) Gap(Ave) Time(Ave) Gap(Ave) Time(Ave) Gap(Ave) τ=1.5 10 0.29 0.00000001 0.42 0.00000000 1.01 0.00019138 100 4.67 0.00000000 8.01 0.00000000 30.45 0.00000003 500 39.92 0.00000001 49.26 0.00000000 205.37 0.00000002 1000 68.11 0.00000001 109.52 0.00000000 437.95 0.00000002 5000 476.95 0.00000000 668.08 0.00000000 4738.70 0.00000002 10000 1242.57 0.00000001 2016.01 0.00000000 15348.57 0.00000002 τ=2 10 0.13 0.00000000 0.18 0.00000000 0.45 0.00008221 100 1.36 0.00000000 2.03 0.00000000 7.56 0.00000001 500 8.40 0.00000000 15.63 0.00000000 53.87 0.00000000 1000 22.38 0.00000000 30.90 0.00000000 108.22 0.00000000 5000 128.18 0.00000000 195.56 0.00000000 815.34 0.00000000 10000 337.17 0.00000000 460.98 0.00000000 2373.23 0.00000000 τ=3 10 0.30 0.00000001 0.42 0.00000000 1.09 0.00003397 100 5.65 0.00000001 10.20 0.00000000 40.08 0.00000016 500 50.36 0.00000001 72.35 0.00000000 225.13 0.00000004 1000 100.17 0.00000001 145.24 0.00000000 463.74 0.00000004 5000 582.84 0.00000002 894.95 0.00000000 4067.13 0.00000002 10000 1715.00 0.00000001 2565.21 0.00000002 13649.88 0.00000005 τ=3.5 10 0.44 0.00000006 0.60 0.00000000 1.81 0.00018065 100 10.90 0.00000001 16.97 0.00000000 60.45 0.00000004 500 80.28 0.00000002 124.50 0.00000002 379.20 0.00000004 1000 171.62 0.00000002 252.79 0.00000002 852.59 0.00000004 5000 1033.28 0.00000002 1700.71 0.00000002 8510.86 0.00000006 10000 2345.25 0.00000002 4682.55 0.00000002 27723.99 0.00000004 Table 3: Computational results for the 0.5-centrum problem with different norms and different dimensions. could solve with our general formulation, in all cases, problem sizes of 1000 points. Accuracy is rather good and the bottleneck here is the size of the SDP object to be handle since the fact that all lambda are non-null makes it impossible to simplify the representation. Finally, we also report in Table 5, for the sake of illustration, an example of application of the result in Theorem 9. This problem consists of the minimization of the difference between the maximum and minimum distances of a number of demand points (nranging between 10 and 1000) with respect to a solution point that must belong to a non-convex feasible region defined by the following two non-convex constraints x2 1−2x2 2−2x2 3≥0and −2x2 1+5x2 2+4x2 3≥0within the unit cube. Clearly, this case is more difficult to solve since this problem is non-convex and thus, we need to resort to the 24 Dim 2 3 10 τn Time(Ave) Gap(Ave) Time(Ave) Gap(Ave) Time(Ave) Gap(Ave) τ=1.5 10 0.24 0.00000000 0.40 0.00000000 1.19 0.00000256 100 4.03 0.00000001 6.73 0.00000001 22.22 0.00000000 500 159.42 0.00000001 190.77 0.00000001 380.99 0.00000000 1000 1270.76 0.00000013 1730.61 0.00000002 2379.03 0.00000000 τ=2 10 0.14 0.00000003 0.18 0.00000001 0.62 0.00002783 100 5.11 0.00000255 6.94 0.00000002 17.07 0.00000003 500 427.54 0.00000314 443.47 0.00000266 1092.70 0.00000054 1000 2079.97 0.00000315 7702.03 0.00000241 9235.59 0.00000073 τ=3 10 0.51 0.00000005 0.72 0.00000010 1.89 0.00000405 100 64.14 0.00000483 58.10 0.00000119 152.45 0.00000148 500 1532.17 0.00018236 2269.39 0.00003466 7950.09 0.00000931 1000 4546.73 0.00025375 5678.17 0.00008893 18011.79 0.00003434 τ=3.5 10 0.63 0.00000130 1.48 0.00001077 7.05 0.00000197 100 33.32 0.00000097 302.76 0.00009247 596.20 0.00024303 500 1555.08 0.00000524 2774.06 0.00035570 8705.00 0.00014431 1000 7625.95 0.00001702 7681.10 0.00059695 18845.92 0.00020324 Table 4: Results for convex ordered median problem with general λ, different norms and dimensions. hierarchy of relaxations introduced in Theorem 9. Nevertheless, we have obtained good results in this case even with the first relaxation order. Dim 3 τn Time(Ave) Gap(Ave) τ=2 10 0.46 0.00001623 100 9.45 0.00457982 500 80.56 0.00030263 1000 204.96 0.00094492 Table 5: Computational results for the Range problem with two nonconvex constraints. 5 Conclusions We develop a unified tool for minimizing convex ordered median location problems in finite dimension and with general ℓτ-norms. We report computational results that show the powerfulness of this methodology to solve medium size continuous location problems. 25