- Email: [email protected]
A new method for solving capacitated location problems based on a set partitioning approach
CAOR 760 SHANKAR BALA CHRIS
Computers & Operations Research 29 (2002) 365}386
A new method for solving capacitated location problems based on a set partitioning approach Roberto Baldacci , Eleni Hadjiconstantinou , Vittorio Maniezzo, Aristide Mingozzi * Imperial College, Management School 53 Princes's Gate, Exhibition Road, London SW7 2PG, UK Department of Mathematics, University of Bologna via Sacchi, 3, 47023 Cesena, Italy Received 1 February 1999; received in revised form 1 February 2000
Abstract We consider the capacitated p-median problem (CPMP) in which a set of n customers must be partitioned into p disjoint clusters so that the total dissimilarity within each cluster is minimized and constraints on maximum cluster capacities are met. The dissimilarity of a cluster is computed as the sum of the dissimilarities existing between each entity of the cluster and the median associated to such cluster. In this paper we present an exact algorithm for solving the CPMP based on a set partitioning formulation of the problem. A valid lower bound to the optimal solution cost is obtained by combining two di!erent heuristic methods for solving the dual of the LP-relaxation of the exact formulation. Computational tests on problems proposed in the literature and on new sets of test problems show the e!ectiveness of the proposed algorithm. Scope and purpose A basic location problem consists of locating a number of facilities or depots to supply a set of customers. The objective is to minimize the cost of locating the facilities and assigning the customers to them. This problem has been extensively studied in the literature and is commonly referred to as the plant location problem, or facility location problem. When each potential facility has a constraint on the maximum demand that it can supply and the number of facilities to locate is speci"ed, the problem is known as the Capacitated p-median problem (CPMP). The purpose of this paper is to present a new exact algorithm for the CPMP and evaluate its computational performance on a set of test problems taken from the literature and on a new set of test problems. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Capacitated p-median; Facility location; Set partitioning problem; Dual heuristic solution
* Corresponding author Tel.: #39-0547-642812; fax: #39-0547-610054. E-mail address: [email protected] (A. Mingozzi). 0305-0548/02/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 5 - 0 5 4 8 ( 0 0 ) 0 0 0 7 2 - 1
CAOR 760 366
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
1. Introduction The capacitated p-median problem (CPMP) is a particular location problem in which a set of n customers must be partitioned into p disjoint clusters so that the total dissimilarity within each cluster is minimized and constraints on maximum cluster capacities are met. The dissimilarity of a cluster is computed as the sum of the dissimilarities existing between each customer of the cluster and the median associated to the cluster. This problem, which appears also under the names of capacitated warehouse location problem, sum-of-stars clustering problem and others, is NP-hard [1]. The CPMP arises in several real-world situations. The most common application occurs in the design of a distribution network, where a set of customers must be supplied goods from plants or warehouses. The main strategic choice becomes that of deciding the location, among the feasible ones, of the suppliers of goods. However, the CPMP also arises in other contexts, for example in the topological design of computer communication networks. In this case it deals with the design process of dividing network nodes into groups, and selecting a concentrator location for each group so that all the nodes in a group can be assigned to the same concentrator without violating its capacity constraints [2]. Moreover, the CPMP arises as a subproblem in several important real-world applications, such as problems in telecommunication network design, vehicle routing and in general location theory (see the book edited by Mirchandani and Francis [3]). The CPMP has been extensively studied in clustering and location theory. The following exact algorithms have been proposed in the literature for the CPMP. Pirkul [2] describes a branch and bound method which uses the Lagrangean relaxation of the partitioning constraints. Also, exact techniques based on a set partitioning with side constraints formulation of the CPMP have been investigated by Neebe and Rao [4] and Hansen et al. [5]. Heuristic algorithms have been proposed by Mulvey and Beck [6] and Pirkul [2]. Metaheuristic approaches are described in Golden and Skiscim [23], in Osman and Christo"des [7] and in Maniezzo et al. [8]. In real-world applications, additional requirements may be imposed on the clusters such as bounds on the clusters cardinality or incompatibility between customers [9]. Other important constraints are the distance constraints, which require that each customer is assigned to a median within a speci"ed distance limit. This prevents solutions in which some customers are too far from the medians. An analysis of network location problems with distance constraints has been presented by Moon and Chaudhry [10], by Choi and Chaudhry [11], by Moon and Papayanopoulos [12] and by Chaudhry et al. [13]. Another problem closely related to the CPMP accepts multiple partial assignments of entities to clusters. A mixed integer formulation of this problem and exact algorithms have been proposed by Christo"des and Beasley [14], Leung and Magnanti [15] and Aardal [16]. Heuristic methods for this problem have been investigated by Van Roy [17] and Beasley [18]. In this paper the CPMP is formulated as a set partitioning problem with a side constraint (SP). Each column of the SP corresponds to a feasible cluster and the additional constraint forces any feasible solution to contain exactly p clusters. In order to obtain a valid lower bound to the CPMP, we use a heuristic procedure to solve the dual problem (DSP) of the linear programming relaxation of formulation SP. This heuristic procedure, called HDSP, combines two di!erent heuristic algorithms, H and H, each one "nding a feasible solution to DSP without requiring the entire set of the dual constraints. Procedure H is
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
367
based on a Lagrangean relaxation and subgradient optimization, while procedure H is based on linear programming optimization and requires the generation of a limited number of variables of formulation SP starting from the dual solution produced by H. The dual solution obtained by H and a valid upper bound to the CPMP are then used to reduce the number of clusters (i.e., the variables of the SP formulation) which may form an optimal solution. The size of the reduced integer problem might still be too large for solving it by a branch and bound method. In this case we propose a procedure, called EHP, that consists in reducing the number of variables of the integer program so that the resulting problem can be solved by an integer programming solver. This method may terminate without having found an optimal solution. However, the EHP procedure provides a means to estimate the maximum deviation from optimality of the CPMP solution obtained. Furthermore, the EHP procedure can easily be extended to deal with additional constraints, such as those introduced above. The paper is structured as follows. Section 2 gives two mathematical formulations of the CPMP and presents the pricing procedure for reducing the variables of the SP model. Section 3 describes the bounding procedure HDSP, while Section 4 presents the method used for generating the clusters. Section 5 describes the EHP method for solving the CPMP and Pirkul's branch and bound method, which we implemented for computational comparisons. Moreover, we describe simple extensions of procedure HDSP and EHP to deal with CPMP problems with additional constraints. Computational results are presented for a number of problems drawn from the literature and for new sets of problems with additional constraints. The results show the e!ectiveness of the proposed method in solving problems up to 200 customers. Conclusions are presented in Section 6. 2. Mathematical formulations of the CPMP Let N"1,2, n be a set of n customers and [d ] be a n;n matrix indicating the dissimilarities GH between pairs of customers of the set N. We assume that d *0 and d "0 for all i, j3N. A positive GH GG integer weight q is associated with each customer i, i3N. Any subset B-N is called cluster. Given G a cluster B, the customer jH3B such that ∀j3BjH d H ) d GH GH GZ GZ is called the median of B and will be denoted with (B). A positive integer weight Q is associated with each customer j, j3N, which denotes the capacity H of j when it is used as the median of a cluster. A cluster B is feasible if q )Q , G L GZ where Q denotes the capacity of median (B). For a given integer p, 2)p)n, a feasible CPMP L solution is represented by a partition S"B , B ,2, B of N into p feasible clusters and its cost is N given by the sum of the cluster dissimilarities, that is N d l , z(S)" GL l GZ l
CAOR 760 368
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
where l denotes the median of cluster Bl . An optimal CPMP solution corresponds to a partitioning of the customer set N into p feasible clusters of minimum cost. Let be a (0}1) variable that is one if and only if customer i is assigned to a cluster whose GH median is j. We assume that "1 means that customer j is chosen to be a median of a cluster. HH A well-known mathematical formulation of the CPMP is as follows. (F) z(F)"Min s.t.
d GH GH GZ, HZ,
(1)
"1, ∀i3N, GH HZ,
(2)
q )Q , ∀j3N, G GH H HH GZ,
(3)
"p, HH HZ, 30,1, ∀i, j3N. GH
(4) (5)
Constraints (2) force each customer to be assigned to a cluster, constraints (3) impose that the total capacity of a median must not be exceeded, constraint (4) speci"es that the total number of clusters must be equal to p and constraints (5) are the integrality constraints. Di!erent relaxations of formulation F have been proposed in the literature to derive lower bounds to CPMP. Mulvey and Beck [6] proposed a Lagrangean relaxation of constraints (2), while Beasley [18] used a Lagrangean relaxation of constraints (2)}(4). The CPMP can also be formulated as a set partitioning problem with an additional constraint as follows: E Let P be the index set of all feasible clusters whose median is customer j, j3N, and let H B"P P 2P . L E Let B be the index set of all clusters containing customer i3N. G E Let cl , Bl and l indicate the cost, the subset of customers and the median of cluster l3B, respectively. E Let xl be a (0}1) variable that is equal to 1 if and only if cluster l3B belongs to the optimal solution. The resulting mathematical formulation, called SP, of problem CPMP is as follows: cl xl Z
(SP) z(SP)"Min
(6)
l B
s.t.
xl "1, ∀i3N, Z
(7)
xl "p, Z xl 30, 1, ∀l3B.
(8)
l B G
l B
(9)
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
369
Eqs. (7) impose that each customer i3N is assigned to one cluster. Eq. (8) forces the solution to contain p clusters. The formulation SP is easily extendible to deal with additional cluster constraints simply by removing from B any infeasible cluster. Hansen et al. [5] have shown that solving the LP-relaxation of formulation F cannot yield a larger lower bound than solving the LP-relaxation of formulation SP. Problem SP cannot be solved directly since the number of columns may be enormous even for CPMP instances of moderate size. In the following, the dual problem (called DSP) of the linear relaxation of SP, is used in order to generate a valid lower bound to the CPMP. The method used to solve DSP is heuristic (which of course does not a!ect the optimality of the "nal CPMP solution) and it does not require the explicit generation of the cluster-index set B. Moreover, the dual solution obtained is used to reduce the set B by removing those clusters that cannot belong to any optimal CPMP solution. Let u , i3N, be the dual variables associated with constraints (7) and w be the dual variable of G constraint (8). The dual of the LP-relaxation of SP, called DSP, is as follows: (DSP) z(DSP)"Max
s.t.
u #pw G GZ,
(10)
u #w)cl , ∀l3B, G GZ l u unrestricted, ∀i3N, G w unrestricted.
(11)
(12)
Like problem SP, problem DSP is impractical to solve, since the number of constraints is equal to the number of the variables in SP. 2.1. Variable reduction of problem SP Let (u, w) be a feasible solution of DSP of cost z(DSP) and let x be a feasible solution of SP of cost z(SP). It is well known, from linear programming duality theory, that z(DSP))z(SP) and, consequently, any feasible solution of DSP provides a valid lower bound to SP. Let cl be the reduced cost of cluster l3B corresponding to the dual solution (u, w), that is cl "cl ! u !w. G GZ l
(13)
This dual solution (u, w) can be used to reduce the number of variables of SP as it is established in the following. Theorem 1. Let S"l: l3B, s.t. xl "1. The following relationship holds: z(SP)"z(DSP)# cl . l Z1Y
(14)
CAOR 760 370
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
Proof. From Eq. (13), we have cl " cl ! u ! w. G l l l Z1Y Z1Y Z1Y GZ l Z1Y Since x is a feasible solution of CPMP, we have l
u # w" u #pw G l G Z1Y GZ l Z1Y GZ, and hence, l
(15) cl " cl ! u !pw. G l GZ, Z1Y Z1Y Noticing that z(SP)" l cl and that z(DSP)" u #pw, from Eq. (15) we obtain Eq. Z1Y GZ, G (14). 䊐 l
Corollary 1. Let z(;B) be the cost of a feasible CPMP solution and (u, w) be a feasible solution of DSP of cost z(DSP). Any optimal solution of SP of cost less than z(;B) cannot contain any cluster l3B whose reduced cost is greater or equal to z(;B)!z(DSP). Proof. It follows directly from Theorem 1. 䊐 Corollary 1 states that an optimal CPMP solution can be obtained by replacing in problem SP the set B with the subset B de"ned as follows: B"l: l3B, s.t. cl (z(UB)!z(DSP),
(16)
therefore, an optimal solution of problem SP can be obtained by replacing the set B with B. However, the size of B might still be too large for solving problem SP, even if the gap z(UB)!z(DSP) is small. In this case we propose to solve problem SP by using a subset F-B containing a limited number of clusters so that the resulting problem is solvable by an integer programming solver (e.g. CPLEX [19]). The optimal solution obtained for the resulting problem SP is not guaranteed to be an optimal CPMP solution. However, the method used to choose the subset F allows us to estimate the distance of the solution obtained from optimality.
3. A heuristic procedure for solving problem DSP In this section we describe a heuristic procedure (called HDSP) for "nding a feasible solution to DSP that is based on the following general idea. A feasible solution w"w#w#2#wI of the linear program: (LP) Max wb s.t.
wA)c, w unrestricted,
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
371
can be obtained by successively solving a sequence of k linear programs LP, LP,2, LPI by means of k di!erent heuristic procedures H, H,2, HI. Procedure HP "nds a feasible solution wP of the linear program LPP de"ned as follows: (LPP) Max wPb s.t.
wPA)cP, wP unrestricted,
where cP"c!(w#w#2#wP\)A and w"0. Assuming w"0, we have c"c and, therefore, the linear program LP corresponds to the original problem LP. This method has been applied by Mingozzi et al. [20] for solving the vehicle routing problem, by Bianco et al. [21] for the multiple depot vehicle scheduling, by Mingozzi et al. [24] for the crew scheduling problem and by Mingozzi et al. [22] for the vehicle routing problem with backhauls. The procedure HDSP involves two heuristics. The "rst, called H, solves problem DSP(,DSP) and does not require the generation of the set B while the second procedure, called H, solves DSP and requires the generation of a limited subset of the set B. 3.1. Procedure H The procedure H is based on the lower bound obtained from formulation F by relaxing the set partitioning constraints (2) in a Lagrangean fashion. Hansen et al. [5] have shown that the optimal solutions of this relaxation and of the LP-relaxation of the SP formulation have the same value. Let "( , ,2, ) be the Lagrangean multipliers associated with constraints (2). The relaxed L problem, called LR(), is as follows: GZ, s.t. GZ,
(LR()) z(LR())"Min
(d ! ) # GH G GH G HZ, GZ, q )Q , ∀j3N, G GH H HH
(3)
"p, HH HZ,
(4)
30, 1, ∀i, j3N. GH
(5)
For a given set of multiplier vector , the problem of solving LR() separates into n independent knapsack problems, called KP (), j"1,2, n, of the form H (KP ()) z(KP ())"Min (d ! ) H H GH G GH GZ, (3) s.t. q )Q , G GH H GZ, 30, 1, ∀i3N. GH
(5)
CAOR 760 372
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
An optimal KP () solution corresponds to a cluster BlHH such that lH"argminl PH cl ! l . H H GZ G Z An optimal solution of LR() is obtained as follows. Let z(KP ()), z(KP ()),2, z(KP N ()) be the p least cost solutions of the n knapsack problems H H H KP (), j"1,2, n. It is easy to see that the optimal LR() solution is given by H 1, k"1,2, p, HI I " HH 0, k"p#1,2, n
and HI "1 if i3BlHHI , HI "0 otherwise, GH GH
∀i3N, k"1,2, p
and HI "0, ∀i3N, k"p#1,2, n. GH The value of the objective function for this solution is z(LR())" N z(KP I ())# . I H GZ, G Theorem 2 shows that the solution of LR(), for any vector , can be transformed into a feasible DSP solution. Theorem 2. Let H be the optimal solution of LR() of cost z(LR()) for a given vector . A feasible solution (u, w) of DSP of cost z(DSP)"z(LR()) is obtained by setting, for every customer i3N:
if H"0, GG u" G G #z(KP ())! if H"1, G G GG w" ,
(17a) (17b)
where "maxz(KP ())H : j3N. H HH Proof. It is su$cient to show that, for any j3N: u#w)cl , ∀l3P . G H GZ l We must consider two cases:
(18)
(a) H "0. HH Using Eqs. (17a) and (17b), inequalities (18) can be written as )cl ! , ∀l3P . G H GZ l As H "0, from the de"nition of , we have HH z(KP ())* . H Moreover, since z(KP ()) is the optimal solution cost of problem KP () we have H H z(KP ()))cl ! , ∀l3P . G H H GZ l
(19)
(20)
(21)
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
373
Hence, from (20) and (21) we obtain inequalities (19). (b) H "1. HH Using Eqs. (17a) and (17b), inequalities (18) can be written as #z(KP ())! # )cl , G H GZ l which is equivalent to
∀l3P H
z(KP ()))cl ! , ∀l3P G H H GZ l and inequality (22) is valid since z(KP ()) is an optimal solution of problem KP (). 䊐 H H
(22)
Algorithm H for solving DSP is an iterative procedure that "nds a feasible solution of the problem Max [z(LR())]. An iteration of H consists of computing a new vector and of "nding H a new solution of the resulting problem LR(). The method used for updating , at each iteration, is as follows. Let H be the optimal solution of LR(). Let h be the number of medians to which customer i3N G is assigned in the solution H, i.e. h " H . In any feasible CPMP solution we have h "1, i3N, G HZ, GH G hence, a subgradient optimization method can be used to change as follows: z(UB)!z(LR()) (h !1), ∀i3N " !
G G G (h !1) IZ, I where z(UB) is a valid upper bound to the optimal CPMP solution cost and is the step size (and is a parameter). A feasible DSP solution (u, w) of cost z(DSP) is given by Eqs. (17) using the vector H that has produced, in a priori "xed number of iterations, the best approximate solution of problem Max [z(LR())]. H 3.2. Procedure H Procedure H is a heuristic procedure based on linear programming that "nds a feasible solution of the following problem: (DSP) z(DSP)"Max s.t.
u #pw G GZ, u #w)cl , ∀l3B, G GZ l u unrestricted, ∀i3N, G w unrestricted,
where cl is the reduced cost of cluster l3B computed according to the dual solution (u, w) produced by procedure H, that is cl "cl ! u!w. G GZ l
CAOR 760 374
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
Problem DSP could not be solved directly as it might involve a huge number of constraints. In this section we describe a procedure called H, for reducing the number of constraints of DSP so that the resulting problem, called RD, can be solved directly and any solution of RD is a feasible DSP solution. Problem RD is obtained from DSP as follows: (i) The number of constraints (11) are reduced by replacing B with a subset F of limited size; (ii) Constraints are added to force any RD solution to satisfy constraints (11) for any l3BF. Reduced problem RD. Let F be a subset of P , j3N, that satis"es the following two conditions: H H cl (z(UB)!z(DSP), ∀l3F H
(23a)
Max [cl ]) Min W c X. l l F l P F Z H Z H H
(23b)
A procedure for computing the sets F , j3N, is presented in Section 4. The reduced dual problem H RD is obtained from DSP by replacing the cluster set B with the subset F" F . HZ, H Problem RD is as follows: (RD) z(RD)"Max s.t.
u #pw G GZ,
u #w)cl , ∀l3F, G GZ l u ); , G G
∀i3N,
w)0,
(24) (25) (26) (27)
where the upper bound ; , i3N, are chosen such that G ; )cl , G GZ l
∀l3BF.
(28)
Theorem 3. Any feasible solution of RD is also a feasible solution of DSP. Proof. Constraints (26) and (28) imply that every feasible solution to RD satis"es the following inequalities: u #w) ; #w)cl #w, ∀l3BF. G G GZ l GZ l The validity of Theorem 3 follows from inequalities (27) and (29). 䊐
(29)
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
375
Problem RD can be considered to be the dual of the following problem RP: cl xl # ; y G G GZ, Z s.t. xl #y "1, ∀i3N, G l F Z G xl #y "p, l F Z xl *0, ∀l3F,
(RP) z(RP)"Min
l F
y *0, ∀i3N0. G Procedure H "nds an optimal solution (xH, yH) of RP of cost zH(RP) and the corresponding optimal dual solution (uH, wH) of RD. Hence, we have z(DSP)"zH(RP) and u"uH, w"wH. An optimal CPMP solution. Procedure HDSP "nds a solution (u, w) of DSP of cost z(DSP)"z(DSP)#z(DSP) by setting u"u#u, w"w#w. The cases where the optimal solution (xH,yH) of RP corresponds to an optimal solution of problem SP are as follows: (a) xH integer, yH"0. This solution is also an optimal CPMP solution of cost zH(SP)"z(DSP). (b) xH not integer, ; "R, ∀i3N, yH"0. G In this case the clusters of any optimal CPMP solution are in the set F, hence, the solution (u, w) is an optimal solution of problem DSP and an optimal CPMP solution can be obtained by solving problem SP by replacing the set B with F. (c) yHO0. The RP solution achieved is not feasible for CPMP even if the variables xH l , l3F, have integer values. In fact, yH'0, for some i3N, implies that vertex i is not fully covered by the clusters G that are in the optimal RP solution, i.e. the set F might not contain an optimal or even a feasible CPMP solution. 3.2.1. The computation of ; , i3N G A valid method for computing ; , i3N, in order to satisfy constraints (28) is as follows. Let G c "Maxl FH cl , ∀j3N and let jH be such that Z H c H /Q H "Min [c /Q ]. H H H H HZ, ; , i3N, is assigned the following value: G (30) ; "q c H /Q H . H G G H It is easy to show that the values assigned to ; , i3N according to expression (30) satisfy G inequalities (28). In fact, for any l3BF, ; " q c H /Q H , G G H H GZ l GZ l
(31)
CAOR 760 376
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
since c H /Q H )c l /Q l )cl /Q l . H H L L L From (31), we obtain ; ) q cl /Q l . G G L GZ l GZ l Moreover, as l q /Q l )1, from inequalities (32) we have GZ G L 䊐 ; )cl , ∀l3BF. G GZ l
(32)
4. Generation of the set F In this section, a procedure, called GEN(F ), that generates, for a given median j, a cluster H set F that satis"es conditions (23a) and (23b) is described. GEN(F ) is derived from an H H exact dynamic programming procedure for generating the set P by limiting the size of the state H space graph in such a way that the states generated at the last stage of the recursion correspond to a set F . H Let c(B)" (d !u)!w be the reduced cost of cluster B for median (B) with respect GZ GL G to the dual solution (u, w) produced by procedure H. Let NM "(i , i ,2, i ) be the ordered set L\ of the (n!1) customers obtained from N by removing the median customer j. From a computational point of view, the vector NM is ordered so that (d I !uI ))(d I> !uI> ), k"1,2, n!2. G H G G GH Consider a feasible cluster B- j, i , i ,2, i and let f (B) be the reduced cost of the feasible I I cluster of minimum reduced cost that can be obtained by expanding B with the customer subset i ,i , ,i . The value of f (B) can be computed as follows: I> I> 2 L\ I f (B)"c(B)#g (B), (33) I I> where L\ g (B)"Min (d P !uP )y I> GH G P PI> L\ s.t. q P y ) Q ! q , H G G P PI> GZ y 30, 1, r"k#1,2, n!1. P We assume g (B)"0 so that f (B)"c(B). It is easy to verify that function f (B) has the following L L\ I properties:
f (B))f (B), k"1,2, n!2, (34) I I> (P2) f (B))f (Bi ), k"1,2, n!2. (35) I I> I> Let R be the set of all feasible clusters that can be obtained by using the "rst k customers I i , i ,2, i of NM and such that I f (B)(z(UB)!z(DSP), ∀B3R . I I (P1)
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
377
We assume that j3B and q )Q , ∀B3R . Every cluster B such that (B)Oj is removed GZ G H I from R . It is evident that the resulting set R is a subset of P , that is: L\ L\ H R "B3P : c(B)(z(UB)!z(DSP), and therefore, R can be used instead of P for L\ H L\ H generating F . However, the size of each R , k"1,2, n!1, might be too large even if the gap H I z(UB)!z(DSP) is small. To overcome this problem, we propose a procedure called GEN(F ) H that generates, for a given median j, a sequence of subsets RM , k"1,2, n!1, satisfying the I following two conditions: RM ) , I (36) Max [f (B)]) Min [f (B)]. I I ZRI RM I ZRM I The cluster set F corresponds to RM after having removed any cluster B such that (B)Oj. H L\ Procedure GEN(F ) makes use of a temporary set ¹ representing, at each stage k, a subset of H I R such that Max I [f (B)])Min RI I [f (B)]. The set RM is then extracted from ¹ . At each I Z2 I Z 2 I I I stage k of GEN(F ), RM -¹ -R . H I I I Algorithm GEN (F ) H Step 0. Set ¹ " j and k"0. Step 1. (Dexne the subset RM -¹ ) I I If ¹ ) , then set RM "¹ and de"ne h "z(UB)!z(DSP). I I I I If ¹ ' , then let RM be the largest subset of ¹ such that RM ) and I I I I ". Max RM I [f (B)])Min I RM I [f (B)]. Set h "Max RM I [f (B)] and ¹ I Z2 I I Z I I> Z Step 2. (Generate ¹ ) I> For any B3RM , consider the two clusters B and B"Bi as possible elements of I I> ¹ . L> If f (B))h , then add B to ¹ . I> I I> If q )Q and f (B))h then add B to ¹ . GZ Y G H I> I I> Step 3. Set k"k#1. If k(n!1, then go to Step 1. Step 4. (Dexne F ) H Remove from RM every cluster B such that (B)Oj. L\ Set F be the largest subset of RM such that F ) and H L\ H "Max FH [c(B)]. Max FH [c(B)])Min RM L\ FH [c(B)]. Set h Z L\ Z Z It is easy to note that h *h *2*h since f (B))h , ∀B3¹ while RM -¹ . L\ I> I I> I> I> In order to prove that every RM , k"1,2, n!1, satis"es conditions (36), it is su$cient to show I that f (B)*h , I I
∀B3R ¹ . I I
(37)
In fact, due to Step 1 if condition (37) holds, f (B)*h , ∀B3R RM and f (B)*h , I I I I I I ∀B3¹ RM . I I Assume that R satis"es condition (36). This is true for k"2 if *2 since R "2. By I\ contradiction, assume that there exists a cluster BH3R ¹ such that f (BH)(h . There are two I I I I
CAOR 760 378
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
cases: 1. i 3BH. For property P2, f (BHi ))f (BH). Since f (BH)(h and h *h , f (BH I I\ I I I I I\ I I\ i )(h and (BHi )3RM . Consider Step 2 at stage k!1. The expansion of cluster I I\ I I\ (BHi )3RM produces BH that is added to ¹ since f (BH)(h and this show the I I\ I I I\ contradiction. 2. i , BH. Due to property P1 and h *h , f (BH)(h and, hence, BH3RM . Therefore, I I\ I I\ I\ I\ at Step 2, the cluster BH is added to ¹ and this shows the contradiction. I 5. Methods for solving the CPMP In this section, two methods for solving the CPMP are described. The "rst, called EHP, consists in reducing the number of variables of the integer program SP so that the resulting problem can be solved by an integer programming solver (e.g. CPLEX). The second is the branch and bound method proposed by Pirkul [2] which has been implemented in order to compare the computational performance of algorithm EHP. 5.1. The EHP procedure Whenever the procedure H ends without having found the optimal solution, it is necessary, as described in Section 2, to solve the following problem SP: (SP) cl xl Z Y s.t. xl "1, ∀i3N, l B B Z G Y xl "p, l B Z Y xl 30, 1, ∀l3B,
z(SP)"Min
l B
where the set B is computed according to expression (16). However, the size of B may be still too large. In this case, the set B is replaced with a subset F-B in such a way that the resulting problem becomes solvable by an integer programming code. The solution achieved might not be an optimal CPMP solution, but it is possible to evaluate its distance from optimality. Let (u, w) be the heuristic solution of DSP of cost z(DSP) produced by the bounding procedure HDSP. And let cl "cl ! l u !w be the reduced cost of cluster l3B with GZ G respect to (u, w). Let F be the subset of P , j3N, that satis"es conditions (23), where the H H reduced costs cl are replaced with cl and z(DSP) is substituted with z(DSP). The subsets F , H j3N, are computed by means of procedure GEN(F ) using the dual solution (u, w) instead of H (u, w).
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
379
Let xH be an optimal integer solution, of cost zH(SP), of the set partitioning problem SP where the set B is replaced by the set F (we assume zH(SP)"R if the set F does not contain any feasible CPMP solution). Notice that, if zH(SP)(R, the solution xH is a feasible CPMP solution. The following recognizes if xH is an optimal CPMP solution. Let ¸ be a lower bound to the reduced cost of the least reduced cost cluster in the set P F , H H H j3N and assume that ¸ "R if F ( . Let GAPMIN"Min [¸ ]. The following two cases H H HZ, H exist: (a) zH(SP))z(DSP)#GAPMIN: xH is optimal for the CPMP, since any feasible CPMP solution involving some cluster of the set B/F would have a cost greater or equal than z(DSP)#GAPMIN. (b) zH(SP)'z(DSP)#GAPMIN: xH might not be an optimal CPMP solution, and z(DSP)#GAPMIN is a valid lower bound to the optimal CPMP solution cost.
5.2. Branch and bound method BB In order to evaluate the computational performance of algorithm EHP we implemented the Branch and Bound procedure (BB) proposed by Pirkul [2]. BB makes use of the lower bound LR() obtained from formulation F by relaxing the set partitioning constraints (2) in a Lagrangean fashion (see Section 3.1). The tree-search is a two level binary tree-search in which the "rst level of the tree is formed by "xing variables , j"1,2, n, which de"ne the HH medians of the solution. Whenever a leaf node of this top level tree is reached, it can be treated as the root of a new sub-tree which is explored by "xing variables , i, j"1,2, n, iOj, which GH correspond to assigning the customers to the medians de"ned at the top level of the tree. At each tree node, a lower bound is computed by means of the subgradient procedure applied to the Lagrangean relaxation LR(). The values of the multipliers are initialized with the penalties associated with the lower bound found at the predecessor node. Fifty subgradient iterations are carried out at each tree node, with the exception of the root node where 300 subgradient iterations are performed.
5.3. Dealing with real-world constraints In constrained clustering, as described in Section 1, additional requirements are imposed on the clusters such as bounds on clusters cardinality, incompatibility between customers, distance constraints, etc. Additional constraints can be easily incorporated both in procedure EHP, rejecting any infeasible cluster, and in method BB, changing the branching strategy of the second level of the tree by rejecting the branches that lead to infeasible solutions. However, these simple modi"cations alone could lead to severely degraded performance if no corresponding change is introduced in the bound. The adjustment of the bound in BB can be very involved for some kind of additional constraints, such as customer incompatibility and maximum cluster cardinality. On the other hand, the bounding procedure HDSP can easily be amended to consider any additional
CAOR 760 380
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
cluster constraint by changing Step 2 of algorithm GEN(F ) to avoid the construction of infeasible H clusters.
6. Computational results The algorithms EHP and BB described in Section 5 have been coded in Fortran 77 and run on a Silicon graphics Indy machine (MIPS R4400/200 MHz processor) on six classes of test problems, called A, B, C, D, E, and F, respectively. Classes A and B are drawn from the literature and classes C, D, E and F are four new sets of test problems generated in order to evaluate the performance of EHP on CPMP instances with large dimensions or additional constraints. CPLEX 4.0 is used as the LP-solver in procedure H and as the integer programming solver in EHP. The problems of classes A and B correspond to the CPMP instances used by Osman and Christo"des [7]; set A contains 10 problems of size n"50 and p"5 while set B contains 10 problems of size n"100 and p"10. In these two sets of problems the dissimilarity matrices correspond to Euclidean distance matrices. Problems of classes C and D have the same structure of problems of classes A and B but are of bigger dimension: n"150 and p"15 for class C and n"200 and p"20 for class D. The problems of classes E and F are derived from problems of class A by imposing additional constraints on the clusters. These constraints are bounds on the cluster cardinality and incompatibilities between entities. The incompatibilities are de"ned by an incompatibility matrix [n ] where IH n "1 if entity k cannot be in the same cluster of entity j, 0 otherwise. We impose a maximum IH cluster cardinality of 6 for problems of class E and 11 for problems of class F. The incompatibility matrices are generated by randomly de"ning "ve incompatibilities in such a way that the optimal CPMP solution of the corresponding class A problem becomes infeasible. Table 1 reports the incompatibilities de"ned for problems of classes E and F. For solving problems of classes E and F, algorithm BB has been modi"ed by changing the branching strategy of the second level of the tree, rejecting the branches that lead to unfeasible solutions. Table 1 Problem classes E and F: incompatibilities between customers Problems
Incompatibilities
E1, F1 E2, F2 E3, F3 E4, F4 E5, F5 E6, F6 E7, F7 E8, F8 E9, F9 E10, F10
3, 1, 2, 2, 3, 2, 2, 3, 2, 1,
7 8, 43 11, 41 19, 47 33, 48 7 3, 19 10, 37 11, 32 20, 50 25 5, 49 8, 39 14, 45 27, 50 42 9, 46 10, 17 11, 16 24, 48 12 6, 48 7, 9 18, 31 21, 44 8 4, 14 10, 46 12, 47 18, 45 6 4, 18 10, 36 15, 49 27, 50 8 4, 6 9, 49 22, 44 23, 46 15 8, 42 9, 41 13, 50 17, 40 6 5, 41 10, 46 15, 50 28, 45
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
381
Fig. 1. Example of CPMP solutions with additional constraints.
Fig. 1 shows examples of CPMP solutions with additional constraints. Fig. 1(a) presents the optimal solution of problem CCPX5 of class A. Figs. 1(b) and (c) display CPMP solutions obtained by considering "ve incompatibilities between customers and by imposing a maximum cardinality of 6 and of 11 for problem E5 of class E (Fig. 1(b)) and for problem F5 of class F (Fig. 1(c)), respectively. In the "gures, the customers with the same letter (A, B, C, D or E) are incompatible, i.e. they cannot be in the same cluster. A time limit of 3600 CPU seconds was imposed on both the EHP and BB algorithms for problems of classes A, B, E and F, while a time limit of 7200 CPU seconds was imposed on the same procedures for problems of classes C and D. The parameter , used in GEN(F ), was set to 2000 in H both procedures H and EHP for all test problems. The results obtained are presented in Tables 2}7. The columns in these tables are de"ned as follows: Probl.: a problem instance identi"er. z(UB): cost of the CPMP solution found by the heuristic algorithm of Maniezzo et al. [8], increased of 1. This was done to evaluate if the EHP and BB algorithms are able, within the imposed time limits, to "nd a solution of cost at least equal to the cost of the solution found by the heuristic. zH(SP): cost of the optimal CPMP solution found by algorithm EHP (or cost of the best solution found). z(DSP): lower bound produced by procedure H after 300 subgradient iterations. t : computing time of bounding procedure H. & %E : percentage error of the lower bound z(DSP) (i.e. %E "100z(DSP)/zH(SP)). "1. "1. z(DSP): "nal lower bound produced by procedure HDSP. t : total computing time of procedure HDS. &"1. %E : percentage error of the lower bound z(DSP) (i.e. %E "100z(DSP)/zH(SP)). &"1. &"1. LS: "z(DSP)#GAPMIN, where GAPMIN is the value de"ned in Section 5.1. (The solution of cost zH(SP) produced by EHP is optimal if zH(SP))LS). F : number of clusters generated in EHP.
CAOR 760 382
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
Table 2 Computational results on problems of class A HDSP
EHP
Probl.
z(UB) z(DSP) t
z(DSP) t
CCPX1 CCPX2 CCPX3 CCPX4 CCPX5 CCPX6 CCPX7 CCPX8 CCPX9 CCPX10
714 741 752 652 666 779 788 821 716 830
705.0 740.0 749.0 651.0 664.0 778.0 778.3 772.1 712.6 818.4
Averages
704.2 739.5 748.2 650.2 663.1 777.6 778.1 770.8 712.1 815.1
%E & "1. 0.9 98.8 0.3 99.9 0.7 99.6 0.4 99.9 0.5 99.9 0.4 100.0 1.7 98.9 1.9 94.0 1.2 99.6 2.3 98.3 1.0
98.9
&"1. 1.4 0.4 0.9 0.6 0.7 0.5 2.6 26.7 1.6 4.5 4.0
%E &"1. 98.9 100.0 99.7 100.0 100.0 100.0 98.9 94.2 99.7 98.7 99.0
BB
zH(SP) LS
F
713 740 751 651 664 778 787 821 715 829
570 10 64 23 19 41 1915 30491 355 2511
715 742 753 653 667 780 788 789 717 834
t
#&. 1.7 0.4 1.0 0.6 0.7 0.6 4.7 3600.0 1.8 12.2 362.4
zH(BB)
t
713 740 751 651 664 778 787 820 715 829
5.4 0.8 2.4 0.7 4.1 2.0 29.3 1990.3 14.3 288.5 233.8
Optimal solution obtained by procedure HDSP. No solution found of cost smaller than z(UB).
t : total computing time of procedure EHP including t . #&. &"1 zH(BB): cost of the optimal CPMP solution found by algorithm BB (or cost of the best solution found). t : computing time of algorithm BB. Table 2 shows that most of the problems of class A are relatively easy, with the only exception of problem CCPX8. In fact, bound z(DSP) is already close to the optimum and bound z(DSP) is capable of "nding four optimal solutions. The good quality of the bound enables both EHP and BB to be very e!ective on all problems, except CCPX8. Note that on this set, BB has been able, within the time limit of 3600 s, of proving the optimality of all instances, while EHP was unable to "nd an optimal solution of CCPX8. Table 3 shows the results of problems of class B. In this problems, z(DSP) does not improve much over z(DSP), as z(DSP) is already very close to the optimal solution cost on all instances. Algorithm EHP turns out to be superior to BB in fact it is able to prove the optimality of eight solutions, while BB is able to solve to optimality only four problems. Furthermore EHP "nds an improved solution of problem CCPX18 without proving its optimality, since its cost is higher than the corresponding LS value. Algorithm BB cannot improve the upper bound of 6 out of the 10 instances. Moreover, when both algorithms "nd an optimal solution, EHP does it in a much smaller computing time. Tables 4 and 5 have been included to show the e!ectiveness of our approach on large size instances. Problems of sizes n"150 and 200 are, in fact, beyond the dimensions currently reported in the literature. The quality of the bound is good, being the instances of these classes structured as those of classes A and B. Algorithm BB fails to improve the initial upper bound solution on 8 over 10 problems, even though for these bigger problems we allowed 7200 CPU seconds. Optimal
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
383
Table 3 Computational results on problems of class B HDSP Probl.
z(UB) z(DSP)
CCPX11 CCPX12 CCPX13 CCPX14 CCPX15 CCPX16 CCPX17 CCPX18 CCPX19 CCPX20
1007 967 1027 983 1092 955 1035 1044 1032 1006
EHP t
& 3.6 4.1 2.8 4.2 4.8 3.2 4.7 4.3 4.3 6.7
1000.1 958.2 1021.1 971.1 1079.1 950.2 1024.1 1031.2 1025.1 972.2
Averages
%E "1. 99.4 99.2 99.5 98.9 98.9 99.6 99.0 98.9 99.4 96.6
z(DSP) t 1001.1 958.5 1021.6 971.8 1079.8 951.3 1025.3 1031.7 1026.3 972.8
4.3 99.0
&"1. 6.4 11.6 3.7 12.0 16.6 4.9 11.8 15.9 7.5 77.3
%E &"1. 99.5 99.2 99.6 99.0 99.0 99.7 99.2 98.9 99.5 96.7
16.8
99.0
BB
zH(SP) LS
F
1006 966 1026 982 1091 954 1034 1043 1031 1006
2551 17458 995 17362 19394 1160 10252 23211 3105 62872
1008 968 1028 984 1092 956 1036 1040 1033 984
t
#&. 11.7 832.8 5.1 304.3 558.4 6.4 280.5 1137.8 22.2 3600.0
zH(BB) t 1006 967 1026 983 1092 954 1035 1044 1031 1013
675.9
1099.3 3600.0 88.3 3600.0 3600.0 1911.6 3601.0 3600.0 3115.8 3600.0 2781.8
This solution was not proved to be optimal. No solution found of cost smaller than z(UB). Table 4 Computational results on problems of class C HDSP
EHP
Probl.
z(UB) z(DSP) t
CCPX21 CCPX22 CCPX23 CCPX24 CCPX25 CCPX26 CCPX27 CCPX28 CCPX29 CCPX30
1290 1292 1220 1236 1189 1228 1270 1181 1260 1243
1270.1 1281.2 1192.1 1222.2 1186.1 1226.1 1260.1 1174.2 1236.2 1229.2
& 12.1 9.8 12.3 9.1 5.1 4.4 8.9 6.4 5.9 9.7
%E "1. 98.5 99.2 97.7 98.9 99.8 99.9 99.3 99.5 98.1 98.9
z(DSP) t &"1. 1270.8 119.5 1282.0 24.7 1192.5 180.2 1222.6 36.6 1188.0 1.2 1227.0 0.9 1260.9 23.1 1174.5 4.9 1236.9 119.3 1230.0 46.3
8.4 99.0
55.6
%E &"1. 98.5 99.3 97.7 98.9 100.0 100.0 99.4 99.5 98.2 99.0 99.0
BB
zH(SP) LS
F
1290 1291 1220 1236 1188 1227 1269 1180 1260 1243
111235 36091 124682 55178 203 210 35235 12425 111573 82852
1279 1291 1199 1229 1191 1229 1268 1181 1246 1236
t
#&. 7200 4597 7200 7200 6 5 3641 35 7200 7200 1657
zH(BB) t 1290a 1292 1220 1236 1188 1227 1270 1181 1260 1243
7200.0 7200.0 7200.0 7200.0 928.1 211.0 7200.0 7200.0 7200.0 7200.0 569.2
No solution found of cost smaller than z(UB). This solution was not proved to be optimal.
solutions can be found by BB only when the lower bound is less than 0.2% from the optimal solution value. In no case was BB able to solve a problem with n"200. On the other hand, EHP is able to solve to optimality four instances with n"150 and two instances with n"200. Furthermore, it "nds the optimal solution, without proving its optimality, of one additional instance in both sets.
CAOR 760 384
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
Table 5 Computational results on problems of class D HDSP
EHP
Probl.
z(UB)
z(DSP) t
CCPX31 CCPX32 CCPX33 CCPX34 CCPX35 CCPX36 CCPX37 CCPX38 CCPX39 CCPX40
1447 1352 1391 1395 1401 1384 1399 1462 1427 1393
1438.2 1335.1 1389.1 1362.1 1380.1 1371.1 1370.1 1438.2 1415.1 1390.2
& 15.5 16.7 6.8 18.2 17.7 16.6 18.4 20.1 16.8 8.1
%E "1. 99.5 98.8 99.9 97.6 98.5 99.1 97.9 98.4 99.2 99.9
15.5
98.9
z(DSP) 1438.7 1335.4 1390.0 1363.2 1380.8 1372.0 1370.9 1438.7 1415.5 1391.5
t &"1. 28.5 56.0 1.2 199.0 80.2 51.2 148.2 145.8 44.4 1.8 75.6
%E &"1. 99.5 98.8 100.0 97.7 98.6 99.1 98.0 98.4 99.2 100.0
zH(SP)
LS
F
1446 1352 1390 1395 1401 1384 1399 1462 1427 1392
1444 1340 1392 1369 1386 1376 1376 1443 1420 1395
42529 72627 299 94585 75159 64271 90667 95636 64984 498
98.9
t
#&. 4272.1 7200.0 8.0 7200.0 7200.0 7200.0 7200.0 7200.0 7200.0 10.1 1430.1
This solution was not proved to be optimal. No solution found of cost smaller than z(UB). Table 6 Computational results on problems of class E HDSP Probl.
z(UB) z(DSP)
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10
477 519 504 464 521 564 536 479 497 544
Averages
423.2 472.2 440.1 422.1 474.1 513.2 487.3 436.1 452.3 495.1
EHP t
& 1.5 1.8 1.7 1.9 2.0 1.8 2.6 1.9 2.8 2.5
%E "1. 88.9 97.6 87.5 94.0 96.8 92.5 92.1 93.4 91.9 94.5
2.0
92.9
z(DSP) t 435.0 480.3 468.6 441.8 490.0 544.9 501.0 461.9 489.7 512.5
&"1. 72.9 2.6 35.2 8.7 3.3 20.9 26.0 5.6 21.1 5.2
20.1
%E
&"1. 91.4 99.2 93.2 98.4 100.0 98.2 94.7 98.9 99.5 97.8 97.1
BB
zH(SP) LS
F
476 484 503 449 490 555 529 467 492 524
46392 240 26695 809 91 1312 16714 298 269 709
476 486 499 468 492 572 531 469 513 526
t
#&. 2644.6 2.7 74.8 9.9 3.3 23.5 369.2 5.7 21.4 5.7 316.1
zH(BB) t 477 484 504 464 491 564 536 479 497 544
3600.0 909.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3330.9
This solution was not proved to be optimal. No solution found of cost smaller than z(UB).
Table 6 presents the results for the instances of class E that include additional constraints. Algorithm EHP is able to prove the optimality of nine solutions over 10, while BB can solve only one problem. Notice that the quality of lower bound z(DSP) is much worse than that on the class A problems (from which the class E problems are drawn) and that z(DSP) gives an important contribution in "lling the gap zH(SP)!z(DSP).
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
385
Table 7 Computational results on problems of class F HDSP Probl.
z(UB) z(DSP)
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10
819 835 857 734 776 878 878 870 804 922
Averages
704.3 739.3 747.4 650.1 663.2 777.3 777.2 770.4 712.1 816.2
EHP t
& 1.5 1.6 1.7 1.6 2.2 1.6 2.2 1.9 2.0 4.2
%E "1. 86.1 93.8 87.3 90.5 85.6 91.6 91.2 90.3 92.7 91.2
z(DSP) t 738.7 758.7 768.7 682.3 723.9 810.9 818.4 801.6 736.4 845.9
4.2 90.0
&"1. 126.1 68.4 138.6 81.5 108.3 87.2 113.6 108.0 85.1 111.1
%E &"1. 90.3 96.3 89.8 95.0 93.4 95.5 96.1 94.0 95.9 94.5
102.8
94.1
BB
zH(SP) LS
F
818 788 856 718 775 849 852 853 768 895
39487 39837 46250 37240 39930 41040 46663 44253 40311 46395
763 771 787 697 742 831 832 818 746 869
t
#&. 478.2 85.9 249.0 94.7 126.1 172.9 131.4 224.3 170.8 568.5
230.2
zH(BB) t 819 835 857 734 776 878 878 870 804 922
3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0
This solution was not proved to be optimal. No solution found of cost smaller than z(UB).
Table 7 shows the results for problems of class F. The problems of class F are harder than the corresponding ones of class E for both algorithms EHP and BB. This re#ects on a worse quality of lower bound H even though z(DSP) continues to contribute signi"cantly over z(DSP). EHP "nds for all problems a better solution than the initial upper bound, but cannot prove the optimality of any of the solutions obtained. Algorithm BB, on the other hand, is unable to improve over the initial upper bound in any of the 10 instances.
7. Conclusions In this paper, a new method for the capacitated p-median problem (CPMP) based on a set partitioning formulation of the problem has been presented. A valid lower bound to the optimal solution cost is obtained by combining two di!erent heuristic methods for solving the dual of the LP-relaxation of the exact formulation. The dual solution obtained is used for generating a reduced set partitioning problem that can be solved by an integer programming solver. The solution achieved might not be an optimal CPMP solution, however the new method allows to estimate its maximum distance from optimality. Computational tests on problems proposed in the literature and on new sets of test problems with large dimensions or additional constraints show the e!ectiveness of the proposed algorithm.
References [1] Garey MR, Johnson DS. Computers and intractability: a guide to the theory of NP completeness. San Francisco: Freeman, 1979.
CAOR 760 386
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
[2] Pirkul H. E$cient algorithms for the capacitated concentrator location problem. Computers and Operations Research 1987;14(3):197}208. [3] Mirchandani PB, Francis RL. Discrete location theory. Chichester: Wiley, 1990. [4] Neebe AW, Rao MR. An algorithm for the "xed charge assigning users to sources problem. Journal of the Operational Research Society 1983;34(11):1107}13. [5] Hansen P, Jaumard B, Sanlaville E. Weight constrained minimum sum-of-star clustering. Gerad Technical Report G-93-38, 1994. [6] Mulvey JM, Beck MP. Solving capacitated clustering problems. European Journal of Operational Research 1984;18:339}48. [7] Osman IH, Christo"des N. Capacitated clustering problems by hybrid simulated annealing and tabu search. International Transactions in Operational Research 1994;1(3):317}36. [8] Maniezzo V, Mingozzi A, Baldacci R. A bionomic approach to the capacitated p-median problem. Journal of Heuristic 1998;4:263}80. [9] Hansen P, Jaumard B. Cluster analysis and mathematical programming. Mathematical Programming 1997;79:191}215. [10] Moon I, Chaudhry S. An analysis of network location } problems with distance constraints. Management Science 1984;30:290}307. [11] Choi I, Chaudhry SS. The p-median problem with maximum distance constraints: a direct approach. Location Science 1993;1(3):235}43. [12] Moon I, Papayanopoulos L. Facility location on a tree with maximum distance constraints. Computers and Operations Research 1995;22:905}14. [13] Chaudhry SS, Choi I, Smith DK. Facility location with and without maximum distance constraints through the p-median problem. International Journal of Operations and Production Management 1995;15(10):76}82. [14] Christo"des N, Beasley JE. Extensions to a Lagrangean relaxation approach for the capacitated warehouse location problem. European Journal of Operational Research 1983;12:19}28. [15] Leung JMY, Magnanti TL. Valid inequalities and facets of the capacitated plant location problem. Mathematical Programming 1989;44:271}91. [16] Aardal K. Capacitated facility location: separation algorithms and computational experience. Technical Report, Department of Econometrics, Tilburg University, 1994. [17] Van Roy TJ. A cross decomposition algorithm for capacitated facility location. Operations Research 1985;34:145}63. [18] Beasley JE. An algorithm for solving large capacitated warehouse location problems. European Journal of Operational Research 1988;33:314}25. [19] CPLEX Optimization Inc. Using the CPLEX callable library and CPLEX mixed integer library, 1996. [20] Mingozzi A, Christo"des N, Hadjiconstantinou E. An exact algorithm for the vehicle routing problem based on the set partitioning formulation. Internal Report Department of Mathematics, University of Bologna, Italy, 1994. [21] Bianco L, Mingozzi A, Ricciardelli S. A set partitioning approach to the multiple depot vehicle scheduling problem. Optimization Methods and Software 1994;3:163}94. [22] Mingozzi A, Baldacci R, Giorgi S. An exact method for the vehicle routing problem with backhauls. Transportation Science 1998;33:329}415. [23] Golden B, Skiscim C. Using simulated annealing to solve routing and location problems. Naval Research Logistic Quarterly 1986;33:261}79. [24] Mingozzi A, Boschetti MA, Bianco L, Ricciarelli S. A set partitioning approach to the crew scheduling problem. Operations Research 1999;47(6):873}88.
Computers & Operations Research 29 (2002) 365}386
A new method for solving capacitated location problems based on a set partitioning approach Roberto Baldacci , Eleni Hadjiconstantinou , Vittorio Maniezzo, Aristide Mingozzi * Imperial College, Management School 53 Princes's Gate, Exhibition Road, London SW7 2PG, UK Department of Mathematics, University of Bologna via Sacchi, 3, 47023 Cesena, Italy Received 1 February 1999; received in revised form 1 February 2000
Abstract We consider the capacitated p-median problem (CPMP) in which a set of n customers must be partitioned into p disjoint clusters so that the total dissimilarity within each cluster is minimized and constraints on maximum cluster capacities are met. The dissimilarity of a cluster is computed as the sum of the dissimilarities existing between each entity of the cluster and the median associated to such cluster. In this paper we present an exact algorithm for solving the CPMP based on a set partitioning formulation of the problem. A valid lower bound to the optimal solution cost is obtained by combining two di!erent heuristic methods for solving the dual of the LP-relaxation of the exact formulation. Computational tests on problems proposed in the literature and on new sets of test problems show the e!ectiveness of the proposed algorithm. Scope and purpose A basic location problem consists of locating a number of facilities or depots to supply a set of customers. The objective is to minimize the cost of locating the facilities and assigning the customers to them. This problem has been extensively studied in the literature and is commonly referred to as the plant location problem, or facility location problem. When each potential facility has a constraint on the maximum demand that it can supply and the number of facilities to locate is speci"ed, the problem is known as the Capacitated p-median problem (CPMP). The purpose of this paper is to present a new exact algorithm for the CPMP and evaluate its computational performance on a set of test problems taken from the literature and on a new set of test problems. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Capacitated p-median; Facility location; Set partitioning problem; Dual heuristic solution
* Corresponding author Tel.: #39-0547-642812; fax: #39-0547-610054. E-mail address: [email protected] (A. Mingozzi). 0305-0548/02/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 5 - 0 5 4 8 ( 0 0 ) 0 0 0 7 2 - 1
CAOR 760 366
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
1. Introduction The capacitated p-median problem (CPMP) is a particular location problem in which a set of n customers must be partitioned into p disjoint clusters so that the total dissimilarity within each cluster is minimized and constraints on maximum cluster capacities are met. The dissimilarity of a cluster is computed as the sum of the dissimilarities existing between each customer of the cluster and the median associated to the cluster. This problem, which appears also under the names of capacitated warehouse location problem, sum-of-stars clustering problem and others, is NP-hard [1]. The CPMP arises in several real-world situations. The most common application occurs in the design of a distribution network, where a set of customers must be supplied goods from plants or warehouses. The main strategic choice becomes that of deciding the location, among the feasible ones, of the suppliers of goods. However, the CPMP also arises in other contexts, for example in the topological design of computer communication networks. In this case it deals with the design process of dividing network nodes into groups, and selecting a concentrator location for each group so that all the nodes in a group can be assigned to the same concentrator without violating its capacity constraints [2]. Moreover, the CPMP arises as a subproblem in several important real-world applications, such as problems in telecommunication network design, vehicle routing and in general location theory (see the book edited by Mirchandani and Francis [3]). The CPMP has been extensively studied in clustering and location theory. The following exact algorithms have been proposed in the literature for the CPMP. Pirkul [2] describes a branch and bound method which uses the Lagrangean relaxation of the partitioning constraints. Also, exact techniques based on a set partitioning with side constraints formulation of the CPMP have been investigated by Neebe and Rao [4] and Hansen et al. [5]. Heuristic algorithms have been proposed by Mulvey and Beck [6] and Pirkul [2]. Metaheuristic approaches are described in Golden and Skiscim [23], in Osman and Christo"des [7] and in Maniezzo et al. [8]. In real-world applications, additional requirements may be imposed on the clusters such as bounds on the clusters cardinality or incompatibility between customers [9]. Other important constraints are the distance constraints, which require that each customer is assigned to a median within a speci"ed distance limit. This prevents solutions in which some customers are too far from the medians. An analysis of network location problems with distance constraints has been presented by Moon and Chaudhry [10], by Choi and Chaudhry [11], by Moon and Papayanopoulos [12] and by Chaudhry et al. [13]. Another problem closely related to the CPMP accepts multiple partial assignments of entities to clusters. A mixed integer formulation of this problem and exact algorithms have been proposed by Christo"des and Beasley [14], Leung and Magnanti [15] and Aardal [16]. Heuristic methods for this problem have been investigated by Van Roy [17] and Beasley [18]. In this paper the CPMP is formulated as a set partitioning problem with a side constraint (SP). Each column of the SP corresponds to a feasible cluster and the additional constraint forces any feasible solution to contain exactly p clusters. In order to obtain a valid lower bound to the CPMP, we use a heuristic procedure to solve the dual problem (DSP) of the linear programming relaxation of formulation SP. This heuristic procedure, called HDSP, combines two di!erent heuristic algorithms, H and H, each one "nding a feasible solution to DSP without requiring the entire set of the dual constraints. Procedure H is
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
367
based on a Lagrangean relaxation and subgradient optimization, while procedure H is based on linear programming optimization and requires the generation of a limited number of variables of formulation SP starting from the dual solution produced by H. The dual solution obtained by H and a valid upper bound to the CPMP are then used to reduce the number of clusters (i.e., the variables of the SP formulation) which may form an optimal solution. The size of the reduced integer problem might still be too large for solving it by a branch and bound method. In this case we propose a procedure, called EHP, that consists in reducing the number of variables of the integer program so that the resulting problem can be solved by an integer programming solver. This method may terminate without having found an optimal solution. However, the EHP procedure provides a means to estimate the maximum deviation from optimality of the CPMP solution obtained. Furthermore, the EHP procedure can easily be extended to deal with additional constraints, such as those introduced above. The paper is structured as follows. Section 2 gives two mathematical formulations of the CPMP and presents the pricing procedure for reducing the variables of the SP model. Section 3 describes the bounding procedure HDSP, while Section 4 presents the method used for generating the clusters. Section 5 describes the EHP method for solving the CPMP and Pirkul's branch and bound method, which we implemented for computational comparisons. Moreover, we describe simple extensions of procedure HDSP and EHP to deal with CPMP problems with additional constraints. Computational results are presented for a number of problems drawn from the literature and for new sets of problems with additional constraints. The results show the e!ectiveness of the proposed method in solving problems up to 200 customers. Conclusions are presented in Section 6. 2. Mathematical formulations of the CPMP Let N"1,2, n be a set of n customers and [d ] be a n;n matrix indicating the dissimilarities GH between pairs of customers of the set N. We assume that d *0 and d "0 for all i, j3N. A positive GH GG integer weight q is associated with each customer i, i3N. Any subset B-N is called cluster. Given G a cluster B, the customer jH3B such that ∀j3BjH d H ) d GH GH GZ GZ is called the median of B and will be denoted with (B). A positive integer weight Q is associated with each customer j, j3N, which denotes the capacity H of j when it is used as the median of a cluster. A cluster B is feasible if q )Q , G L GZ where Q denotes the capacity of median (B). For a given integer p, 2)p)n, a feasible CPMP L solution is represented by a partition S"B , B ,2, B of N into p feasible clusters and its cost is N given by the sum of the cluster dissimilarities, that is N d l , z(S)" GL l GZ l
CAOR 760 368
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
where l denotes the median of cluster Bl . An optimal CPMP solution corresponds to a partitioning of the customer set N into p feasible clusters of minimum cost. Let be a (0}1) variable that is one if and only if customer i is assigned to a cluster whose GH median is j. We assume that "1 means that customer j is chosen to be a median of a cluster. HH A well-known mathematical formulation of the CPMP is as follows. (F) z(F)"Min s.t.
d GH GH GZ, HZ,
(1)
"1, ∀i3N, GH HZ,
(2)
q )Q , ∀j3N, G GH H HH GZ,
(3)
"p, HH HZ, 30,1, ∀i, j3N. GH
(4) (5)
Constraints (2) force each customer to be assigned to a cluster, constraints (3) impose that the total capacity of a median must not be exceeded, constraint (4) speci"es that the total number of clusters must be equal to p and constraints (5) are the integrality constraints. Di!erent relaxations of formulation F have been proposed in the literature to derive lower bounds to CPMP. Mulvey and Beck [6] proposed a Lagrangean relaxation of constraints (2), while Beasley [18] used a Lagrangean relaxation of constraints (2)}(4). The CPMP can also be formulated as a set partitioning problem with an additional constraint as follows: E Let P be the index set of all feasible clusters whose median is customer j, j3N, and let H B"P P 2P . L E Let B be the index set of all clusters containing customer i3N. G E Let cl , Bl and l indicate the cost, the subset of customers and the median of cluster l3B, respectively. E Let xl be a (0}1) variable that is equal to 1 if and only if cluster l3B belongs to the optimal solution. The resulting mathematical formulation, called SP, of problem CPMP is as follows: cl xl Z
(SP) z(SP)"Min
(6)
l B
s.t.
xl "1, ∀i3N, Z
(7)
xl "p, Z xl 30, 1, ∀l3B.
(8)
l B G
l B
(9)
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
369
Eqs. (7) impose that each customer i3N is assigned to one cluster. Eq. (8) forces the solution to contain p clusters. The formulation SP is easily extendible to deal with additional cluster constraints simply by removing from B any infeasible cluster. Hansen et al. [5] have shown that solving the LP-relaxation of formulation F cannot yield a larger lower bound than solving the LP-relaxation of formulation SP. Problem SP cannot be solved directly since the number of columns may be enormous even for CPMP instances of moderate size. In the following, the dual problem (called DSP) of the linear relaxation of SP, is used in order to generate a valid lower bound to the CPMP. The method used to solve DSP is heuristic (which of course does not a!ect the optimality of the "nal CPMP solution) and it does not require the explicit generation of the cluster-index set B. Moreover, the dual solution obtained is used to reduce the set B by removing those clusters that cannot belong to any optimal CPMP solution. Let u , i3N, be the dual variables associated with constraints (7) and w be the dual variable of G constraint (8). The dual of the LP-relaxation of SP, called DSP, is as follows: (DSP) z(DSP)"Max
s.t.
u #pw G GZ,
(10)
u #w)cl , ∀l3B, G GZ l u unrestricted, ∀i3N, G w unrestricted.
(11)
(12)
Like problem SP, problem DSP is impractical to solve, since the number of constraints is equal to the number of the variables in SP. 2.1. Variable reduction of problem SP Let (u, w) be a feasible solution of DSP of cost z(DSP) and let x be a feasible solution of SP of cost z(SP). It is well known, from linear programming duality theory, that z(DSP))z(SP) and, consequently, any feasible solution of DSP provides a valid lower bound to SP. Let cl be the reduced cost of cluster l3B corresponding to the dual solution (u, w), that is cl "cl ! u !w. G GZ l
(13)
This dual solution (u, w) can be used to reduce the number of variables of SP as it is established in the following. Theorem 1. Let S"l: l3B, s.t. xl "1. The following relationship holds: z(SP)"z(DSP)# cl . l Z1Y
(14)
CAOR 760 370
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
Proof. From Eq. (13), we have cl " cl ! u ! w. G l l l Z1Y Z1Y Z1Y GZ l Z1Y Since x is a feasible solution of CPMP, we have l
u # w" u #pw G l G Z1Y GZ l Z1Y GZ, and hence, l
(15) cl " cl ! u !pw. G l GZ, Z1Y Z1Y Noticing that z(SP)" l cl and that z(DSP)" u #pw, from Eq. (15) we obtain Eq. Z1Y GZ, G (14). 䊐 l
Corollary 1. Let z(;B) be the cost of a feasible CPMP solution and (u, w) be a feasible solution of DSP of cost z(DSP). Any optimal solution of SP of cost less than z(;B) cannot contain any cluster l3B whose reduced cost is greater or equal to z(;B)!z(DSP). Proof. It follows directly from Theorem 1. 䊐 Corollary 1 states that an optimal CPMP solution can be obtained by replacing in problem SP the set B with the subset B de"ned as follows: B"l: l3B, s.t. cl (z(UB)!z(DSP),
(16)
therefore, an optimal solution of problem SP can be obtained by replacing the set B with B. However, the size of B might still be too large for solving problem SP, even if the gap z(UB)!z(DSP) is small. In this case we propose to solve problem SP by using a subset F-B containing a limited number of clusters so that the resulting problem is solvable by an integer programming solver (e.g. CPLEX [19]). The optimal solution obtained for the resulting problem SP is not guaranteed to be an optimal CPMP solution. However, the method used to choose the subset F allows us to estimate the distance of the solution obtained from optimality.
3. A heuristic procedure for solving problem DSP In this section we describe a heuristic procedure (called HDSP) for "nding a feasible solution to DSP that is based on the following general idea. A feasible solution w"w#w#2#wI of the linear program: (LP) Max wb s.t.
wA)c, w unrestricted,
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
371
can be obtained by successively solving a sequence of k linear programs LP, LP,2, LPI by means of k di!erent heuristic procedures H, H,2, HI. Procedure HP "nds a feasible solution wP of the linear program LPP de"ned as follows: (LPP) Max wPb s.t.
wPA)cP, wP unrestricted,
where cP"c!(w#w#2#wP\)A and w"0. Assuming w"0, we have c"c and, therefore, the linear program LP corresponds to the original problem LP. This method has been applied by Mingozzi et al. [20] for solving the vehicle routing problem, by Bianco et al. [21] for the multiple depot vehicle scheduling, by Mingozzi et al. [24] for the crew scheduling problem and by Mingozzi et al. [22] for the vehicle routing problem with backhauls. The procedure HDSP involves two heuristics. The "rst, called H, solves problem DSP(,DSP) and does not require the generation of the set B while the second procedure, called H, solves DSP and requires the generation of a limited subset of the set B. 3.1. Procedure H The procedure H is based on the lower bound obtained from formulation F by relaxing the set partitioning constraints (2) in a Lagrangean fashion. Hansen et al. [5] have shown that the optimal solutions of this relaxation and of the LP-relaxation of the SP formulation have the same value. Let "( , ,2, ) be the Lagrangean multipliers associated with constraints (2). The relaxed L problem, called LR(), is as follows: GZ, s.t. GZ,
(LR()) z(LR())"Min
(d ! ) # GH G GH G HZ, GZ, q )Q , ∀j3N, G GH H HH
(3)
"p, HH HZ,
(4)
30, 1, ∀i, j3N. GH
(5)
For a given set of multiplier vector , the problem of solving LR() separates into n independent knapsack problems, called KP (), j"1,2, n, of the form H (KP ()) z(KP ())"Min (d ! ) H H GH G GH GZ, (3) s.t. q )Q , G GH H GZ, 30, 1, ∀i3N. GH
(5)
CAOR 760 372
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
An optimal KP () solution corresponds to a cluster BlHH such that lH"argminl PH cl ! l . H H GZ G Z An optimal solution of LR() is obtained as follows. Let z(KP ()), z(KP ()),2, z(KP N ()) be the p least cost solutions of the n knapsack problems H H H KP (), j"1,2, n. It is easy to see that the optimal LR() solution is given by H 1, k"1,2, p, HI I " HH 0, k"p#1,2, n
and HI "1 if i3BlHHI , HI "0 otherwise, GH GH
∀i3N, k"1,2, p
and HI "0, ∀i3N, k"p#1,2, n. GH The value of the objective function for this solution is z(LR())" N z(KP I ())# . I H GZ, G Theorem 2 shows that the solution of LR(), for any vector , can be transformed into a feasible DSP solution. Theorem 2. Let H be the optimal solution of LR() of cost z(LR()) for a given vector . A feasible solution (u, w) of DSP of cost z(DSP)"z(LR()) is obtained by setting, for every customer i3N:
if H"0, GG u" G G #z(KP ())! if H"1, G G GG w" ,
(17a) (17b)
where "maxz(KP ())H : j3N. H HH Proof. It is su$cient to show that, for any j3N: u#w)cl , ∀l3P . G H GZ l We must consider two cases:
(18)
(a) H "0. HH Using Eqs. (17a) and (17b), inequalities (18) can be written as )cl ! , ∀l3P . G H GZ l As H "0, from the de"nition of , we have HH z(KP ())* . H Moreover, since z(KP ()) is the optimal solution cost of problem KP () we have H H z(KP ()))cl ! , ∀l3P . G H H GZ l
(19)
(20)
(21)
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
373
Hence, from (20) and (21) we obtain inequalities (19). (b) H "1. HH Using Eqs. (17a) and (17b), inequalities (18) can be written as #z(KP ())! # )cl , G H GZ l which is equivalent to
∀l3P H
z(KP ()))cl ! , ∀l3P G H H GZ l and inequality (22) is valid since z(KP ()) is an optimal solution of problem KP (). 䊐 H H
(22)
Algorithm H for solving DSP is an iterative procedure that "nds a feasible solution of the problem Max [z(LR())]. An iteration of H consists of computing a new vector and of "nding H a new solution of the resulting problem LR(). The method used for updating , at each iteration, is as follows. Let H be the optimal solution of LR(). Let h be the number of medians to which customer i3N G is assigned in the solution H, i.e. h " H . In any feasible CPMP solution we have h "1, i3N, G HZ, GH G hence, a subgradient optimization method can be used to change as follows: z(UB)!z(LR()) (h !1), ∀i3N " !
G G G (h !1) IZ, I where z(UB) is a valid upper bound to the optimal CPMP solution cost and is the step size (and is a parameter). A feasible DSP solution (u, w) of cost z(DSP) is given by Eqs. (17) using the vector H that has produced, in a priori "xed number of iterations, the best approximate solution of problem Max [z(LR())]. H 3.2. Procedure H Procedure H is a heuristic procedure based on linear programming that "nds a feasible solution of the following problem: (DSP) z(DSP)"Max s.t.
u #pw G GZ, u #w)cl , ∀l3B, G GZ l u unrestricted, ∀i3N, G w unrestricted,
where cl is the reduced cost of cluster l3B computed according to the dual solution (u, w) produced by procedure H, that is cl "cl ! u!w. G GZ l
CAOR 760 374
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
Problem DSP could not be solved directly as it might involve a huge number of constraints. In this section we describe a procedure called H, for reducing the number of constraints of DSP so that the resulting problem, called RD, can be solved directly and any solution of RD is a feasible DSP solution. Problem RD is obtained from DSP as follows: (i) The number of constraints (11) are reduced by replacing B with a subset F of limited size; (ii) Constraints are added to force any RD solution to satisfy constraints (11) for any l3BF. Reduced problem RD. Let F be a subset of P , j3N, that satis"es the following two conditions: H H cl (z(UB)!z(DSP), ∀l3F H
(23a)
Max [cl ]) Min W c X. l l F l P F Z H Z H H
(23b)
A procedure for computing the sets F , j3N, is presented in Section 4. The reduced dual problem H RD is obtained from DSP by replacing the cluster set B with the subset F" F . HZ, H Problem RD is as follows: (RD) z(RD)"Max s.t.
u #pw G GZ,
u #w)cl , ∀l3F, G GZ l u ); , G G
∀i3N,
w)0,
(24) (25) (26) (27)
where the upper bound ; , i3N, are chosen such that G ; )cl , G GZ l
∀l3BF.
(28)
Theorem 3. Any feasible solution of RD is also a feasible solution of DSP. Proof. Constraints (26) and (28) imply that every feasible solution to RD satis"es the following inequalities: u #w) ; #w)cl #w, ∀l3BF. G G GZ l GZ l The validity of Theorem 3 follows from inequalities (27) and (29). 䊐
(29)
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
375
Problem RD can be considered to be the dual of the following problem RP: cl xl # ; y G G GZ, Z s.t. xl #y "1, ∀i3N, G l F Z G xl #y "p, l F Z xl *0, ∀l3F,
(RP) z(RP)"Min
l F
y *0, ∀i3N0. G Procedure H "nds an optimal solution (xH, yH) of RP of cost zH(RP) and the corresponding optimal dual solution (uH, wH) of RD. Hence, we have z(DSP)"zH(RP) and u"uH, w"wH. An optimal CPMP solution. Procedure HDSP "nds a solution (u, w) of DSP of cost z(DSP)"z(DSP)#z(DSP) by setting u"u#u, w"w#w. The cases where the optimal solution (xH,yH) of RP corresponds to an optimal solution of problem SP are as follows: (a) xH integer, yH"0. This solution is also an optimal CPMP solution of cost zH(SP)"z(DSP). (b) xH not integer, ; "R, ∀i3N, yH"0. G In this case the clusters of any optimal CPMP solution are in the set F, hence, the solution (u, w) is an optimal solution of problem DSP and an optimal CPMP solution can be obtained by solving problem SP by replacing the set B with F. (c) yHO0. The RP solution achieved is not feasible for CPMP even if the variables xH l , l3F, have integer values. In fact, yH'0, for some i3N, implies that vertex i is not fully covered by the clusters G that are in the optimal RP solution, i.e. the set F might not contain an optimal or even a feasible CPMP solution. 3.2.1. The computation of ; , i3N G A valid method for computing ; , i3N, in order to satisfy constraints (28) is as follows. Let G c "Maxl FH cl , ∀j3N and let jH be such that Z H c H /Q H "Min [c /Q ]. H H H H HZ, ; , i3N, is assigned the following value: G (30) ; "q c H /Q H . H G G H It is easy to show that the values assigned to ; , i3N according to expression (30) satisfy G inequalities (28). In fact, for any l3BF, ; " q c H /Q H , G G H H GZ l GZ l
(31)
CAOR 760 376
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
since c H /Q H )c l /Q l )cl /Q l . H H L L L From (31), we obtain ; ) q cl /Q l . G G L GZ l GZ l Moreover, as l q /Q l )1, from inequalities (32) we have GZ G L 䊐 ; )cl , ∀l3BF. G GZ l
(32)
4. Generation of the set F In this section, a procedure, called GEN(F ), that generates, for a given median j, a cluster H set F that satis"es conditions (23a) and (23b) is described. GEN(F ) is derived from an H H exact dynamic programming procedure for generating the set P by limiting the size of the state H space graph in such a way that the states generated at the last stage of the recursion correspond to a set F . H Let c(B)" (d !u)!w be the reduced cost of cluster B for median (B) with respect GZ GL G to the dual solution (u, w) produced by procedure H. Let NM "(i , i ,2, i ) be the ordered set L\ of the (n!1) customers obtained from N by removing the median customer j. From a computational point of view, the vector NM is ordered so that (d I !uI ))(d I> !uI> ), k"1,2, n!2. G H G G GH Consider a feasible cluster B- j, i , i ,2, i and let f (B) be the reduced cost of the feasible I I cluster of minimum reduced cost that can be obtained by expanding B with the customer subset i ,i , ,i . The value of f (B) can be computed as follows: I> I> 2 L\ I f (B)"c(B)#g (B), (33) I I> where L\ g (B)"Min (d P !uP )y I> GH G P PI> L\ s.t. q P y ) Q ! q , H G G P PI> GZ y 30, 1, r"k#1,2, n!1. P We assume g (B)"0 so that f (B)"c(B). It is easy to verify that function f (B) has the following L L\ I properties:
f (B))f (B), k"1,2, n!2, (34) I I> (P2) f (B))f (Bi ), k"1,2, n!2. (35) I I> I> Let R be the set of all feasible clusters that can be obtained by using the "rst k customers I i , i ,2, i of NM and such that I f (B)(z(UB)!z(DSP), ∀B3R . I I (P1)
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
377
We assume that j3B and q )Q , ∀B3R . Every cluster B such that (B)Oj is removed GZ G H I from R . It is evident that the resulting set R is a subset of P , that is: L\ L\ H R "B3P : c(B)(z(UB)!z(DSP), and therefore, R can be used instead of P for L\ H L\ H generating F . However, the size of each R , k"1,2, n!1, might be too large even if the gap H I z(UB)!z(DSP) is small. To overcome this problem, we propose a procedure called GEN(F ) H that generates, for a given median j, a sequence of subsets RM , k"1,2, n!1, satisfying the I following two conditions: RM ) , I (36) Max [f (B)]) Min [f (B)]. I I ZRI RM I ZRM I The cluster set F corresponds to RM after having removed any cluster B such that (B)Oj. H L\ Procedure GEN(F ) makes use of a temporary set ¹ representing, at each stage k, a subset of H I R such that Max I [f (B)])Min RI I [f (B)]. The set RM is then extracted from ¹ . At each I Z2 I Z 2 I I I stage k of GEN(F ), RM -¹ -R . H I I I Algorithm GEN (F ) H Step 0. Set ¹ " j and k"0. Step 1. (Dexne the subset RM -¹ ) I I If ¹ ) , then set RM "¹ and de"ne h "z(UB)!z(DSP). I I I I If ¹ ' , then let RM be the largest subset of ¹ such that RM ) and I I I I ". Max RM I [f (B)])Min I RM I [f (B)]. Set h "Max RM I [f (B)] and ¹ I Z2 I I Z I I> Z Step 2. (Generate ¹ ) I> For any B3RM , consider the two clusters B and B"Bi as possible elements of I I> ¹ . L> If f (B))h , then add B to ¹ . I> I I> If q )Q and f (B))h then add B to ¹ . GZ Y G H I> I I> Step 3. Set k"k#1. If k(n!1, then go to Step 1. Step 4. (Dexne F ) H Remove from RM every cluster B such that (B)Oj. L\ Set F be the largest subset of RM such that F ) and H L\ H "Max FH [c(B)]. Max FH [c(B)])Min RM L\ FH [c(B)]. Set h Z L\ Z Z It is easy to note that h *h *2*h since f (B))h , ∀B3¹ while RM -¹ . L\ I> I I> I> I> In order to prove that every RM , k"1,2, n!1, satis"es conditions (36), it is su$cient to show I that f (B)*h , I I
∀B3R ¹ . I I
(37)
In fact, due to Step 1 if condition (37) holds, f (B)*h , ∀B3R RM and f (B)*h , I I I I I I ∀B3¹ RM . I I Assume that R satis"es condition (36). This is true for k"2 if *2 since R "2. By I\ contradiction, assume that there exists a cluster BH3R ¹ such that f (BH)(h . There are two I I I I
CAOR 760 378
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
cases: 1. i 3BH. For property P2, f (BHi ))f (BH). Since f (BH)(h and h *h , f (BH I I\ I I I I I\ I I\ i )(h and (BHi )3RM . Consider Step 2 at stage k!1. The expansion of cluster I I\ I I\ (BHi )3RM produces BH that is added to ¹ since f (BH)(h and this show the I I\ I I I\ contradiction. 2. i , BH. Due to property P1 and h *h , f (BH)(h and, hence, BH3RM . Therefore, I I\ I I\ I\ I\ at Step 2, the cluster BH is added to ¹ and this shows the contradiction. I 5. Methods for solving the CPMP In this section, two methods for solving the CPMP are described. The "rst, called EHP, consists in reducing the number of variables of the integer program SP so that the resulting problem can be solved by an integer programming solver (e.g. CPLEX). The second is the branch and bound method proposed by Pirkul [2] which has been implemented in order to compare the computational performance of algorithm EHP. 5.1. The EHP procedure Whenever the procedure H ends without having found the optimal solution, it is necessary, as described in Section 2, to solve the following problem SP: (SP) cl xl Z Y s.t. xl "1, ∀i3N, l B B Z G Y xl "p, l B Z Y xl 30, 1, ∀l3B,
z(SP)"Min
l B
where the set B is computed according to expression (16). However, the size of B may be still too large. In this case, the set B is replaced with a subset F-B in such a way that the resulting problem becomes solvable by an integer programming code. The solution achieved might not be an optimal CPMP solution, but it is possible to evaluate its distance from optimality. Let (u, w) be the heuristic solution of DSP of cost z(DSP) produced by the bounding procedure HDSP. And let cl "cl ! l u !w be the reduced cost of cluster l3B with GZ G respect to (u, w). Let F be the subset of P , j3N, that satis"es conditions (23), where the H H reduced costs cl are replaced with cl and z(DSP) is substituted with z(DSP). The subsets F , H j3N, are computed by means of procedure GEN(F ) using the dual solution (u, w) instead of H (u, w).
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
379
Let xH be an optimal integer solution, of cost zH(SP), of the set partitioning problem SP where the set B is replaced by the set F (we assume zH(SP)"R if the set F does not contain any feasible CPMP solution). Notice that, if zH(SP)(R, the solution xH is a feasible CPMP solution. The following recognizes if xH is an optimal CPMP solution. Let ¸ be a lower bound to the reduced cost of the least reduced cost cluster in the set P F , H H H j3N and assume that ¸ "R if F ( . Let GAPMIN"Min [¸ ]. The following two cases H H HZ, H exist: (a) zH(SP))z(DSP)#GAPMIN: xH is optimal for the CPMP, since any feasible CPMP solution involving some cluster of the set B/F would have a cost greater or equal than z(DSP)#GAPMIN. (b) zH(SP)'z(DSP)#GAPMIN: xH might not be an optimal CPMP solution, and z(DSP)#GAPMIN is a valid lower bound to the optimal CPMP solution cost.
5.2. Branch and bound method BB In order to evaluate the computational performance of algorithm EHP we implemented the Branch and Bound procedure (BB) proposed by Pirkul [2]. BB makes use of the lower bound LR() obtained from formulation F by relaxing the set partitioning constraints (2) in a Lagrangean fashion (see Section 3.1). The tree-search is a two level binary tree-search in which the "rst level of the tree is formed by "xing variables , j"1,2, n, which de"ne the HH medians of the solution. Whenever a leaf node of this top level tree is reached, it can be treated as the root of a new sub-tree which is explored by "xing variables , i, j"1,2, n, iOj, which GH correspond to assigning the customers to the medians de"ned at the top level of the tree. At each tree node, a lower bound is computed by means of the subgradient procedure applied to the Lagrangean relaxation LR(). The values of the multipliers are initialized with the penalties associated with the lower bound found at the predecessor node. Fifty subgradient iterations are carried out at each tree node, with the exception of the root node where 300 subgradient iterations are performed.
5.3. Dealing with real-world constraints In constrained clustering, as described in Section 1, additional requirements are imposed on the clusters such as bounds on clusters cardinality, incompatibility between customers, distance constraints, etc. Additional constraints can be easily incorporated both in procedure EHP, rejecting any infeasible cluster, and in method BB, changing the branching strategy of the second level of the tree by rejecting the branches that lead to infeasible solutions. However, these simple modi"cations alone could lead to severely degraded performance if no corresponding change is introduced in the bound. The adjustment of the bound in BB can be very involved for some kind of additional constraints, such as customer incompatibility and maximum cluster cardinality. On the other hand, the bounding procedure HDSP can easily be amended to consider any additional
CAOR 760 380
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
cluster constraint by changing Step 2 of algorithm GEN(F ) to avoid the construction of infeasible H clusters.
6. Computational results The algorithms EHP and BB described in Section 5 have been coded in Fortran 77 and run on a Silicon graphics Indy machine (MIPS R4400/200 MHz processor) on six classes of test problems, called A, B, C, D, E, and F, respectively. Classes A and B are drawn from the literature and classes C, D, E and F are four new sets of test problems generated in order to evaluate the performance of EHP on CPMP instances with large dimensions or additional constraints. CPLEX 4.0 is used as the LP-solver in procedure H and as the integer programming solver in EHP. The problems of classes A and B correspond to the CPMP instances used by Osman and Christo"des [7]; set A contains 10 problems of size n"50 and p"5 while set B contains 10 problems of size n"100 and p"10. In these two sets of problems the dissimilarity matrices correspond to Euclidean distance matrices. Problems of classes C and D have the same structure of problems of classes A and B but are of bigger dimension: n"150 and p"15 for class C and n"200 and p"20 for class D. The problems of classes E and F are derived from problems of class A by imposing additional constraints on the clusters. These constraints are bounds on the cluster cardinality and incompatibilities between entities. The incompatibilities are de"ned by an incompatibility matrix [n ] where IH n "1 if entity k cannot be in the same cluster of entity j, 0 otherwise. We impose a maximum IH cluster cardinality of 6 for problems of class E and 11 for problems of class F. The incompatibility matrices are generated by randomly de"ning "ve incompatibilities in such a way that the optimal CPMP solution of the corresponding class A problem becomes infeasible. Table 1 reports the incompatibilities de"ned for problems of classes E and F. For solving problems of classes E and F, algorithm BB has been modi"ed by changing the branching strategy of the second level of the tree, rejecting the branches that lead to unfeasible solutions. Table 1 Problem classes E and F: incompatibilities between customers Problems
Incompatibilities
E1, F1 E2, F2 E3, F3 E4, F4 E5, F5 E6, F6 E7, F7 E8, F8 E9, F9 E10, F10
3, 1, 2, 2, 3, 2, 2, 3, 2, 1,
7 8, 43 11, 41 19, 47 33, 48 7 3, 19 10, 37 11, 32 20, 50 25 5, 49 8, 39 14, 45 27, 50 42 9, 46 10, 17 11, 16 24, 48 12 6, 48 7, 9 18, 31 21, 44 8 4, 14 10, 46 12, 47 18, 45 6 4, 18 10, 36 15, 49 27, 50 8 4, 6 9, 49 22, 44 23, 46 15 8, 42 9, 41 13, 50 17, 40 6 5, 41 10, 46 15, 50 28, 45
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
381
Fig. 1. Example of CPMP solutions with additional constraints.
Fig. 1 shows examples of CPMP solutions with additional constraints. Fig. 1(a) presents the optimal solution of problem CCPX5 of class A. Figs. 1(b) and (c) display CPMP solutions obtained by considering "ve incompatibilities between customers and by imposing a maximum cardinality of 6 and of 11 for problem E5 of class E (Fig. 1(b)) and for problem F5 of class F (Fig. 1(c)), respectively. In the "gures, the customers with the same letter (A, B, C, D or E) are incompatible, i.e. they cannot be in the same cluster. A time limit of 3600 CPU seconds was imposed on both the EHP and BB algorithms for problems of classes A, B, E and F, while a time limit of 7200 CPU seconds was imposed on the same procedures for problems of classes C and D. The parameter , used in GEN(F ), was set to 2000 in H both procedures H and EHP for all test problems. The results obtained are presented in Tables 2}7. The columns in these tables are de"ned as follows: Probl.: a problem instance identi"er. z(UB): cost of the CPMP solution found by the heuristic algorithm of Maniezzo et al. [8], increased of 1. This was done to evaluate if the EHP and BB algorithms are able, within the imposed time limits, to "nd a solution of cost at least equal to the cost of the solution found by the heuristic. zH(SP): cost of the optimal CPMP solution found by algorithm EHP (or cost of the best solution found). z(DSP): lower bound produced by procedure H after 300 subgradient iterations. t : computing time of bounding procedure H. & %E : percentage error of the lower bound z(DSP) (i.e. %E "100z(DSP)/zH(SP)). "1. "1. z(DSP): "nal lower bound produced by procedure HDSP. t : total computing time of procedure HDS. &"1. %E : percentage error of the lower bound z(DSP) (i.e. %E "100z(DSP)/zH(SP)). &"1. &"1. LS: "z(DSP)#GAPMIN, where GAPMIN is the value de"ned in Section 5.1. (The solution of cost zH(SP) produced by EHP is optimal if zH(SP))LS). F : number of clusters generated in EHP.
CAOR 760 382
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
Table 2 Computational results on problems of class A HDSP
EHP
Probl.
z(UB) z(DSP) t
z(DSP) t
CCPX1 CCPX2 CCPX3 CCPX4 CCPX5 CCPX6 CCPX7 CCPX8 CCPX9 CCPX10
714 741 752 652 666 779 788 821 716 830
705.0 740.0 749.0 651.0 664.0 778.0 778.3 772.1 712.6 818.4
Averages
704.2 739.5 748.2 650.2 663.1 777.6 778.1 770.8 712.1 815.1
%E & "1. 0.9 98.8 0.3 99.9 0.7 99.6 0.4 99.9 0.5 99.9 0.4 100.0 1.7 98.9 1.9 94.0 1.2 99.6 2.3 98.3 1.0
98.9
&"1. 1.4 0.4 0.9 0.6 0.7 0.5 2.6 26.7 1.6 4.5 4.0
%E &"1. 98.9 100.0 99.7 100.0 100.0 100.0 98.9 94.2 99.7 98.7 99.0
BB
zH(SP) LS
F
713 740 751 651 664 778 787 821 715 829
570 10 64 23 19 41 1915 30491 355 2511
715 742 753 653 667 780 788 789 717 834
t
#&. 1.7 0.4 1.0 0.6 0.7 0.6 4.7 3600.0 1.8 12.2 362.4
zH(BB)
t
713 740 751 651 664 778 787 820 715 829
5.4 0.8 2.4 0.7 4.1 2.0 29.3 1990.3 14.3 288.5 233.8
Optimal solution obtained by procedure HDSP. No solution found of cost smaller than z(UB).
t : total computing time of procedure EHP including t . #&. &"1 zH(BB): cost of the optimal CPMP solution found by algorithm BB (or cost of the best solution found). t : computing time of algorithm BB. Table 2 shows that most of the problems of class A are relatively easy, with the only exception of problem CCPX8. In fact, bound z(DSP) is already close to the optimum and bound z(DSP) is capable of "nding four optimal solutions. The good quality of the bound enables both EHP and BB to be very e!ective on all problems, except CCPX8. Note that on this set, BB has been able, within the time limit of 3600 s, of proving the optimality of all instances, while EHP was unable to "nd an optimal solution of CCPX8. Table 3 shows the results of problems of class B. In this problems, z(DSP) does not improve much over z(DSP), as z(DSP) is already very close to the optimal solution cost on all instances. Algorithm EHP turns out to be superior to BB in fact it is able to prove the optimality of eight solutions, while BB is able to solve to optimality only four problems. Furthermore EHP "nds an improved solution of problem CCPX18 without proving its optimality, since its cost is higher than the corresponding LS value. Algorithm BB cannot improve the upper bound of 6 out of the 10 instances. Moreover, when both algorithms "nd an optimal solution, EHP does it in a much smaller computing time. Tables 4 and 5 have been included to show the e!ectiveness of our approach on large size instances. Problems of sizes n"150 and 200 are, in fact, beyond the dimensions currently reported in the literature. The quality of the bound is good, being the instances of these classes structured as those of classes A and B. Algorithm BB fails to improve the initial upper bound solution on 8 over 10 problems, even though for these bigger problems we allowed 7200 CPU seconds. Optimal
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
383
Table 3 Computational results on problems of class B HDSP Probl.
z(UB) z(DSP)
CCPX11 CCPX12 CCPX13 CCPX14 CCPX15 CCPX16 CCPX17 CCPX18 CCPX19 CCPX20
1007 967 1027 983 1092 955 1035 1044 1032 1006
EHP t
& 3.6 4.1 2.8 4.2 4.8 3.2 4.7 4.3 4.3 6.7
1000.1 958.2 1021.1 971.1 1079.1 950.2 1024.1 1031.2 1025.1 972.2
Averages
%E "1. 99.4 99.2 99.5 98.9 98.9 99.6 99.0 98.9 99.4 96.6
z(DSP) t 1001.1 958.5 1021.6 971.8 1079.8 951.3 1025.3 1031.7 1026.3 972.8
4.3 99.0
&"1. 6.4 11.6 3.7 12.0 16.6 4.9 11.8 15.9 7.5 77.3
%E &"1. 99.5 99.2 99.6 99.0 99.0 99.7 99.2 98.9 99.5 96.7
16.8
99.0
BB
zH(SP) LS
F
1006 966 1026 982 1091 954 1034 1043 1031 1006
2551 17458 995 17362 19394 1160 10252 23211 3105 62872
1008 968 1028 984 1092 956 1036 1040 1033 984
t
#&. 11.7 832.8 5.1 304.3 558.4 6.4 280.5 1137.8 22.2 3600.0
zH(BB) t 1006 967 1026 983 1092 954 1035 1044 1031 1013
675.9
1099.3 3600.0 88.3 3600.0 3600.0 1911.6 3601.0 3600.0 3115.8 3600.0 2781.8
This solution was not proved to be optimal. No solution found of cost smaller than z(UB). Table 4 Computational results on problems of class C HDSP
EHP
Probl.
z(UB) z(DSP) t
CCPX21 CCPX22 CCPX23 CCPX24 CCPX25 CCPX26 CCPX27 CCPX28 CCPX29 CCPX30
1290 1292 1220 1236 1189 1228 1270 1181 1260 1243
1270.1 1281.2 1192.1 1222.2 1186.1 1226.1 1260.1 1174.2 1236.2 1229.2
& 12.1 9.8 12.3 9.1 5.1 4.4 8.9 6.4 5.9 9.7
%E "1. 98.5 99.2 97.7 98.9 99.8 99.9 99.3 99.5 98.1 98.9
z(DSP) t &"1. 1270.8 119.5 1282.0 24.7 1192.5 180.2 1222.6 36.6 1188.0 1.2 1227.0 0.9 1260.9 23.1 1174.5 4.9 1236.9 119.3 1230.0 46.3
8.4 99.0
55.6
%E &"1. 98.5 99.3 97.7 98.9 100.0 100.0 99.4 99.5 98.2 99.0 99.0
BB
zH(SP) LS
F
1290 1291 1220 1236 1188 1227 1269 1180 1260 1243
111235 36091 124682 55178 203 210 35235 12425 111573 82852
1279 1291 1199 1229 1191 1229 1268 1181 1246 1236
t
#&. 7200 4597 7200 7200 6 5 3641 35 7200 7200 1657
zH(BB) t 1290a 1292 1220 1236 1188 1227 1270 1181 1260 1243
7200.0 7200.0 7200.0 7200.0 928.1 211.0 7200.0 7200.0 7200.0 7200.0 569.2
No solution found of cost smaller than z(UB). This solution was not proved to be optimal.
solutions can be found by BB only when the lower bound is less than 0.2% from the optimal solution value. In no case was BB able to solve a problem with n"200. On the other hand, EHP is able to solve to optimality four instances with n"150 and two instances with n"200. Furthermore, it "nds the optimal solution, without proving its optimality, of one additional instance in both sets.
CAOR 760 384
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
Table 5 Computational results on problems of class D HDSP
EHP
Probl.
z(UB)
z(DSP) t
CCPX31 CCPX32 CCPX33 CCPX34 CCPX35 CCPX36 CCPX37 CCPX38 CCPX39 CCPX40
1447 1352 1391 1395 1401 1384 1399 1462 1427 1393
1438.2 1335.1 1389.1 1362.1 1380.1 1371.1 1370.1 1438.2 1415.1 1390.2
& 15.5 16.7 6.8 18.2 17.7 16.6 18.4 20.1 16.8 8.1
%E "1. 99.5 98.8 99.9 97.6 98.5 99.1 97.9 98.4 99.2 99.9
15.5
98.9
z(DSP) 1438.7 1335.4 1390.0 1363.2 1380.8 1372.0 1370.9 1438.7 1415.5 1391.5
t &"1. 28.5 56.0 1.2 199.0 80.2 51.2 148.2 145.8 44.4 1.8 75.6
%E &"1. 99.5 98.8 100.0 97.7 98.6 99.1 98.0 98.4 99.2 100.0
zH(SP)
LS
F
1446 1352 1390 1395 1401 1384 1399 1462 1427 1392
1444 1340 1392 1369 1386 1376 1376 1443 1420 1395
42529 72627 299 94585 75159 64271 90667 95636 64984 498
98.9
t
#&. 4272.1 7200.0 8.0 7200.0 7200.0 7200.0 7200.0 7200.0 7200.0 10.1 1430.1
This solution was not proved to be optimal. No solution found of cost smaller than z(UB). Table 6 Computational results on problems of class E HDSP Probl.
z(UB) z(DSP)
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10
477 519 504 464 521 564 536 479 497 544
Averages
423.2 472.2 440.1 422.1 474.1 513.2 487.3 436.1 452.3 495.1
EHP t
& 1.5 1.8 1.7 1.9 2.0 1.8 2.6 1.9 2.8 2.5
%E "1. 88.9 97.6 87.5 94.0 96.8 92.5 92.1 93.4 91.9 94.5
2.0
92.9
z(DSP) t 435.0 480.3 468.6 441.8 490.0 544.9 501.0 461.9 489.7 512.5
&"1. 72.9 2.6 35.2 8.7 3.3 20.9 26.0 5.6 21.1 5.2
20.1
%E
&"1. 91.4 99.2 93.2 98.4 100.0 98.2 94.7 98.9 99.5 97.8 97.1
BB
zH(SP) LS
F
476 484 503 449 490 555 529 467 492 524
46392 240 26695 809 91 1312 16714 298 269 709
476 486 499 468 492 572 531 469 513 526
t
#&. 2644.6 2.7 74.8 9.9 3.3 23.5 369.2 5.7 21.4 5.7 316.1
zH(BB) t 477 484 504 464 491 564 536 479 497 544
3600.0 909.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3330.9
This solution was not proved to be optimal. No solution found of cost smaller than z(UB).
Table 6 presents the results for the instances of class E that include additional constraints. Algorithm EHP is able to prove the optimality of nine solutions over 10, while BB can solve only one problem. Notice that the quality of lower bound z(DSP) is much worse than that on the class A problems (from which the class E problems are drawn) and that z(DSP) gives an important contribution in "lling the gap zH(SP)!z(DSP).
CAOR 760 R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
385
Table 7 Computational results on problems of class F HDSP Probl.
z(UB) z(DSP)
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10
819 835 857 734 776 878 878 870 804 922
Averages
704.3 739.3 747.4 650.1 663.2 777.3 777.2 770.4 712.1 816.2
EHP t
& 1.5 1.6 1.7 1.6 2.2 1.6 2.2 1.9 2.0 4.2
%E "1. 86.1 93.8 87.3 90.5 85.6 91.6 91.2 90.3 92.7 91.2
z(DSP) t 738.7 758.7 768.7 682.3 723.9 810.9 818.4 801.6 736.4 845.9
4.2 90.0
&"1. 126.1 68.4 138.6 81.5 108.3 87.2 113.6 108.0 85.1 111.1
%E &"1. 90.3 96.3 89.8 95.0 93.4 95.5 96.1 94.0 95.9 94.5
102.8
94.1
BB
zH(SP) LS
F
818 788 856 718 775 849 852 853 768 895
39487 39837 46250 37240 39930 41040 46663 44253 40311 46395
763 771 787 697 742 831 832 818 746 869
t
#&. 478.2 85.9 249.0 94.7 126.1 172.9 131.4 224.3 170.8 568.5
230.2
zH(BB) t 819 835 857 734 776 878 878 870 804 922
3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0
This solution was not proved to be optimal. No solution found of cost smaller than z(UB).
Table 7 shows the results for problems of class F. The problems of class F are harder than the corresponding ones of class E for both algorithms EHP and BB. This re#ects on a worse quality of lower bound H even though z(DSP) continues to contribute signi"cantly over z(DSP). EHP "nds for all problems a better solution than the initial upper bound, but cannot prove the optimality of any of the solutions obtained. Algorithm BB, on the other hand, is unable to improve over the initial upper bound in any of the 10 instances.
7. Conclusions In this paper, a new method for the capacitated p-median problem (CPMP) based on a set partitioning formulation of the problem has been presented. A valid lower bound to the optimal solution cost is obtained by combining two di!erent heuristic methods for solving the dual of the LP-relaxation of the exact formulation. The dual solution obtained is used for generating a reduced set partitioning problem that can be solved by an integer programming solver. The solution achieved might not be an optimal CPMP solution, however the new method allows to estimate its maximum distance from optimality. Computational tests on problems proposed in the literature and on new sets of test problems with large dimensions or additional constraints show the e!ectiveness of the proposed algorithm.
References [1] Garey MR, Johnson DS. Computers and intractability: a guide to the theory of NP completeness. San Francisco: Freeman, 1979.
CAOR 760 386
R. Baldacci et al. / Computers & Operations Research 29 (2002) 365}386
[2] Pirkul H. E$cient algorithms for the capacitated concentrator location problem. Computers and Operations Research 1987;14(3):197}208. [3] Mirchandani PB, Francis RL. Discrete location theory. Chichester: Wiley, 1990. [4] Neebe AW, Rao MR. An algorithm for the "xed charge assigning users to sources problem. Journal of the Operational Research Society 1983;34(11):1107}13. [5] Hansen P, Jaumard B, Sanlaville E. Weight constrained minimum sum-of-star clustering. Gerad Technical Report G-93-38, 1994. [6] Mulvey JM, Beck MP. Solving capacitated clustering problems. European Journal of Operational Research 1984;18:339}48. [7] Osman IH, Christo"des N. Capacitated clustering problems by hybrid simulated annealing and tabu search. International Transactions in Operational Research 1994;1(3):317}36. [8] Maniezzo V, Mingozzi A, Baldacci R. A bionomic approach to the capacitated p-median problem. Journal of Heuristic 1998;4:263}80. [9] Hansen P, Jaumard B. Cluster analysis and mathematical programming. Mathematical Programming 1997;79:191}215. [10] Moon I, Chaudhry S. An analysis of network location } problems with distance constraints. Management Science 1984;30:290}307. [11] Choi I, Chaudhry SS. The p-median problem with maximum distance constraints: a direct approach. Location Science 1993;1(3):235}43. [12] Moon I, Papayanopoulos L. Facility location on a tree with maximum distance constraints. Computers and Operations Research 1995;22:905}14. [13] Chaudhry SS, Choi I, Smith DK. Facility location with and without maximum distance constraints through the p-median problem. International Journal of Operations and Production Management 1995;15(10):76}82. [14] Christo"des N, Beasley JE. Extensions to a Lagrangean relaxation approach for the capacitated warehouse location problem. European Journal of Operational Research 1983;12:19}28. [15] Leung JMY, Magnanti TL. Valid inequalities and facets of the capacitated plant location problem. Mathematical Programming 1989;44:271}91. [16] Aardal K. Capacitated facility location: separation algorithms and computational experience. Technical Report, Department of Econometrics, Tilburg University, 1994. [17] Van Roy TJ. A cross decomposition algorithm for capacitated facility location. Operations Research 1985;34:145}63. [18] Beasley JE. An algorithm for solving large capacitated warehouse location problems. European Journal of Operational Research 1988;33:314}25. [19] CPLEX Optimization Inc. Using the CPLEX callable library and CPLEX mixed integer library, 1996. [20] Mingozzi A, Christo"des N, Hadjiconstantinou E. An exact algorithm for the vehicle routing problem based on the set partitioning formulation. Internal Report Department of Mathematics, University of Bologna, Italy, 1994. [21] Bianco L, Mingozzi A, Ricciardelli S. A set partitioning approach to the multiple depot vehicle scheduling problem. Optimization Methods and Software 1994;3:163}94. [22] Mingozzi A, Baldacci R, Giorgi S. An exact method for the vehicle routing problem with backhauls. Transportation Science 1998;33:329}415. [23] Golden B, Skiscim C. Using simulated annealing to solve routing and location problems. Naval Research Logistic Quarterly 1986;33:261}79. [24] Mingozzi A, Boschetti MA, Bianco L, Ricciarelli S. A set partitioning approach to the crew scheduling problem. Operations Research 1999;47(6):873}88.