Solving the median problem with continuous demand on a network

Solving the median problem with continuous demand on a network∗ Rafael Blanquero and Emilio Carrizosa, Universidad de Sevilla, Spain {rblanquero,ecarrizosa}@us.es November 5, 2012 Abstract Where to locate one or several facilities on a network so as to minimize the expected users-closest facility transportation cost is a problem well studied in the OR literature under the name of median problem. In the median problem users are usually identified with nodes of the network. In many situations, however, such assumption is unrealistic, since users should be better considered to be distributed also along the edges of the transportation network. In this paper we address the median problem with demand distributed along edges and nodes. This leads to a globaloptimization problem, which can be solved to optimality by means of a branch-and-bound with DC bounds. Our computational experience shows that the problem is solved in short time even for large instances. Keywords: Network Location, Median Problem, Continuous Demand, DC Functions, Global Optimization 1 Introduction Location problems on networks have attracted the interest of researchers and practitioners since the 60s of last century. Under the usual assumption, the demand is concentrated at the nodes of the network, and the facilities can be located either at the nodes or along the edges of the transportation network. Assuming that the demand is concentrated at the nodes is not realistic when modeling socalled networks spatial phenomena, e.g. [18], that is, phenomena which do happen in points along edges of the network, such as, for instance, traffic accidents, or close to such edges, as happens in urban settings, where edges model the city streets, close to the buildings where demand happens. See [18, 19] for further discussion on the advantages on continuous network models against traditional approaches (discrete or planar location models). ∗Research supported by grants from the Spanish Ministry of Education, Culture and Sport MTM2009-14039- C06-06, Junta de Andaluc´ıa TIC-6064, FQM-329, in part financed by the European Regional Development Fund (ERDF). 1 For this reason, several researchers have addressed location problems on networks under the assumption that demand is not only concentrated on nodes, but also is (continuously) distributed along the edges of the network. Most papers consider the case in which, for each edge of the network, demand is uniformly distributed. [17] addresses the problem of minimizing the expected distance from the users to one facility, [6, 8] consider the two-facility case on very particular network topologies (trees), whereas [15] addresses the minimization of the variance of distances from the users to the facility. Assuming uniform demands on edges should be seen as a first step towards gaining realism of the model, while maintaining tractability. Indeed, the resulting objective functions admits a rather simple form, as a piecewise polynomial function in one variable, whose optimization is reduced to inspecting all critical points, namely, extreme points and points at which the derivative of the polynomial function vanishes. Assuming more general distributions for the demand has been advocated by several authors. [18] suggests the use of general density distributions, from which random samples are generated, yielding a discrete approximation to the problem, which is the one which is later analyzed. Statistical kernel methods have also been recently proposed to model the demand, [20, 23], though, as far as the authors know, no optimization has been carried out, excepting, as said above, discretization via simulation. In this paper we consider a single-facility location problem, namely, the 1-median problem: the point minimizing the expected distance from the users is sought. With respect to the state-of-the-art, we give a further step towards realism, by assuming arbitrary distributions for the demand along edges. Contrary to the planar 1-median problem, known to be convex, [7], the 1-median problem on networks with continuous demand poses nontrivial challenges: we have no longer a simple expression for the objective, and its optimization calls for the use of global-optimization techniques. In particular, it is shown that the objective function is DC, [10, 11, 21, 22] i.e., it can be written as the difference of two convex functions, and thus its optimization can be addressed via branch-and-bound methods customized for DC functions, [1, 4, 5]. The remainder of the paper is organized as follows. In Section 2 the 1-median problem with demand distributed on edges and on nodes of the network is formally introduced. Properties of the objective function are discussed in Section 3, where it is shown, in particular, how the function can be expressed as the difference of two convex functions on each edge of the network. These properties will be the cornerstone for a branch-and-bound algorithm, as described in Section 4. Our numerical experience is reported in Section 5, showing that our algorithm enables us to solve problems on large networks in reasonable time. 2 Problem formulation Let N= (A, E) be a connected and undirected network, with node set A={a1, . . . , an}and edge set E, with |E|=m. Let us denote by lij the length of each edge eij = [ai, aj]∈Eand let d(x, y) be the distance between two points x, y ∈N, obtained as the shortest path from xto y. In particular, the distance dij =d(ai, aj) between each pair of nodes {ai, aj}can be worked out by using standard algorithms, [2]. Note that for every pair of nodes ai, ajwith [ai, aj]∈E, it follows that dij ≤lij ,and equality holds if and only if the edge [ai, aj] is the shortest path joining aiand aj. 2 Given a node ak∈Aand a point x∈[ai, aj], obtained after covering a distance lxon the edge [ai, aj], the distance d(x, ak) from akto xis, as a function of x, a piecewise linear concave function given by: d(x, ak) = min{rk ij(x), sk ij(x)}(1) where rk ij(x) = d(ai, ak) + lxsk ij(x) = d(aj, ak)+(lij −lx) (2) We assume that the demand not only occurs at nodes but also along the edges of the network. More precisely, the demand of a node a∈Awill be denoted by ωa≥0,the total demand of a given edge e∈Eis pe≥0,and it is distributed along eaccording to a random variable with cumulative distribution function (cdf) Fe. Under the previous assumptions, the median problem with continuous demand can be written as follows: min x∈NH(x) := X a∈A wad(x, a) + X e∈E peZy∈e d(x, y)dFe(y) (3) Observe that we are making no assumption on the type of distribution followed for the demand. In case the demand on the edges is continuously distributed along e, i.e., when the cdf Fehas a pdf fe,(3) can be rewritten as min x∈NH(x) := X a∈A wad(x, a) + X e∈E peZy∈e d(x, y)fe(y)dy. (4) 3 Properties It has been noticed in [12] that the objective function of (3) is neither convex nor concave as a rule. Indeed, the objective function can exhibit local optima which are not globally optimal, as the following example shows. So the use of global optimization techniques is required if one seeks the optimal solution. Example 1 Let us consider the network N= (A, E)with A={a1, a2, a3}and E={[a1, a2],[a1, a3],[a2, a3]}. Arc lengths, node demands and arc demands are the following: l12 = 1 l13 = 1 l23 = 1 w1= 0 w2= 0 w3= 0 p12 = 0.35 p13 = 0.30 p23 = 0.35 We also assume that the demand along each arc eis distributed according to a beta distribution, i.e., the probability density function fehas the form fe(x) = Γ(αe+βe) Γ(αe)Γ(βe)xαe−1(1 −x)βe−1x∈[0,1]. The parameters αe, βeof these probability distributions in the three edges are as follows: Arc e αeβe [a1, a2] 0.6 0.5 [a1, a3] 0.4 0.8 [a2, a3] 0.5 0.5 3 Under these assumptions, the objective function of Problem (3) restricted to the edge [a1, a2] takes the form H12(x)=0.35 Z1 0 |x−y|f12(y)dy+0.30 Z1 0 min{x+y, 3−x−y}f13(y)dy+0.35 Z1 0 min{2+x−y, 1−x+y}f23(y)dy. (5) Multimodality of the function is clearly seen in Figure 1. 0.685 0.69 0.695 0.7 0.705 0.71 0 0.2 0.4 0.6 0.8 1 Figure 1: Objective function H12(x) in Example 1 The special structure of the objective function of Problem (3) will be exploited in order to design a deterministic global optimization algorithm that allows us to find an optimal solution to the problem. More precisely, we will show that H(x) in (3) belongs to the broad class of DC functions, [11, 10, 21]. This key property will allow us to solve Problem (3) by branch-and- bound algorithms, as the one described in Section 4, since lower and upper bounds can easily be obtained for DC functions as soon as a DC decomposition is available. Definition 2 Let Ω⊂Rnbe a convex set. A function h: Ω →Ris called DC in Ωif there exist two convex functions h+: Ω →R,h−: Ω →Rsuch that h(x) = h+(x)−h−(x)∀x∈Ω (6) A pair (h+, h−)satisfying (6) is called a DC decomposition of hin Ω. An interesting property of the class of DC functions is that it is closed under the most common operations in optimization, [3, 10, 11, 21, 22]. In particular, if h1, . . . , hrare DC functions and λi∈R, i = 1, . . . , r, then Pr i=1 λihi, maxi=1,...,r hi, mini=1,...,r hiand k(h1, . . . , hr)kare also DC and their DC decompositions can be easily obtained from the DC decompositions of hi, i = 1, . . . , r. The following result shows that the objective function H(x) in (3) is DC on each arc of the network. 4 Proposition 3 Given ¯e∈E, the function H¯e: ¯e7→ Rdefined as H¯e(x) := X a∈A wad(x, a) + X e∈E peZy∈e d(x, y)dFe(y) is DC. A DC decomposition of H¯eon ¯eis given by the pair (H+ ¯e, H− ¯e),with H+ ¯e(x) = p¯eZy∈¯e |x−y|dF¯e(y) H− ¯e(x) = H+ ¯e(x)−H¯e(x) (7) Proof. Let us rewrite H¯e(x) as H¯e(x) = h1(x) + h2(x) + h3(x) where h1(x) = X a∈A wad(x, a) (8) h2(x) = X e∈E,e6=¯e peZy∈e d(x, y)dFe(y) (9) h3(x) = p¯eZy∈¯e d(x, y)dF¯e(y) (10) The function h1is piecewise linear and concave (see [13] for instance), and h2is also concave, see [12]. In what follows it is shown that h3is DC, and a DC decomposition is given. Given x, y ∈¯e= [ai, aj],the distance d(x, y) between xand yis obtained as the minimum of the lengths of the following paths: 1. the subedge of ¯ewith x, y as end points, 2. the subedge of ¯ejoining xand ai, the shortest path joining aiand aj,and then the subedge of ¯ejoining ajand y, 3. the subedge of ¯ejoining xand aj, the shortest path joining ajand ai,and then the subedge of ¯ejoining aiand y. In other words d(x, y) can be expressed as d(x, y) = min {|x−y|, x +dij + (lij −y),(lij −x) + dij +y}(11) = min {|x−y|, lij +dij − |x−y|} (12) =|x−y| − max {0,2|x−y| − (lij +dij)}.(13) Observe that the last expression gives a DC decomposition of don ¯e. Hence, h3can be written as h3(x) = p¯eZy∈¯e |x−y|dF¯e(y)−p¯eZy∈¯e max {0,2|x−y| − (lij +dij)}dF¯e(y),(14) which yields a DC decomposition of h3.Taking into account that H¯e=h1+h2+h3,with h1, h2 concave and h3decomposed as a difference of convex functions in (14), it follows that (7) gives a DC decomposition of H¯e,as asserted. 2 We end this section with further properties of the objective function Heunder some assumption on the edge e. These results extend previous well-known results for particular topologies, e.g. for networks which are chains or trees. 5 Corollary 4 Let ¯e∈Ebe an edge such that p¯e= 0.Then H¯eis concave on ¯e. Proof. If p¯e= 0,then, by Proposition 3, (0,−H¯e) is a valid DC decomposition for H¯e.2 Proposition 5 Let ¯e∈Ebe an edge such that E\ {¯e}is a disconnected network. Then H¯eis convex on ¯e. 4 The algorithm Problem (3) will be solved by using a standard branch and bound method [14, 16] which will find out the optimal solution within a relative accuracy of ε > 0. The bounds of the objective function required to applying the algorithm will be worked out by taking into account its DC structure. A brief description of such an algorithm is shown next. •Phase 1: Initialization 1. Fix the required accuracy ε > 0. 2. Set UB = +∞(upper bound initialization). 3. Compute the all-pairs distance matrix. 4. Set the list Λ of remaining segments as empty. 5. For each edge e∈Edo: (a) Consider eas a segment with its nodes as the segment vertices. (b) Evaluate the objective function at the segment midpoint. If this value is lower than UB, then update UB and store in xUB the midpoint as incumbent. (c) Calculate a lower bound for the segment e,LB(e). (d) If LB(e)< UB/(1 + ε), then insert einto Λ. •Phase 2: Branch and Bound process Repeat as long as no stop was reached: 1. Select from Λ the minimum lower bound segment emin and remove if from Λ. 2. If LB(emin)≥UB/(1 + ε), then stop the algorithm with xUB as optimal solution and UB as optimal objective value. 3. Split emin by its midpoint into two smaller segments, e1 min and e2 min. 4. Evaluate the objective function at the midpoint of the two small segments. If any of these values is lower than UB, then update UB. 5. Compute a lower bound of the objective function on each small segment. 6. If LB(ei min)< UB/(1 + ε) for i= 1 or i= 2, then insert ei min into Λ. 7. If UB has been updated in this iteration, then discard all segments from Λ whose lower bound is greater than UB. 6 The algorithm uses a data structure Λ where all the segments (bits of edge) that can contain an optimal solution are stored. The loop in Phase 1 establishes the initial composition of the data structure by selecting the edges whose lower bound is not greater that the global upper bound of the optimal objective value. At the same time, that upper bound is improved by evaluating each edge at its middle point and replacing the bound with the objective value when this is smaller. Phase 2 of the algorithm consists of an undefined loop where the segment with the worst lower bound is processed; the algorithm finishes when the difference between that lower bound and the global upper bound is smaller than the tolerance εchosen in Phase 1. If the stopping rule is not fulfilled, the selected segment is split into two equal segments which are processed in the same way that the edges in Phase 1. Eventually, if a change in the upper bound took place during an iteration, all the segments in the data structure that cannot contain an optimal solution are removed. The computation of the objective function’s lower bound on each segment requires more attention and is going to be detailed next. 4.1 Constructing lower bounds Given an edge ¯e= [ai, aj]∈E, Proposition 3 provides a DC decomposition for H¯e–the restriction of the objective Hto e– and this fact can be exploited in order to obtain the lower bounds required in the previous algorithm. Starting from the DC representation (7), a concave underestimate L¯e(x) of H¯eis obtained by replacing H+ ¯ewith an affine underestimate built in the usual way, L¯e(x) = H+ ¯e(x0) + ξ(x−x0)−H− ¯e(x) where ξis any point in ∂H+ ¯e(x0),the subdifferential of H+ ¯eat x0∈¯e. Since ∂H+ ¯e(x0) = Z(0,x0) dF¯e(y)+[−1,1] Z{x0} dF¯e(y)−Z(x0,l¯e) dF¯e(y), one has 2F¯e(x0)−1∈∂H+ ¯e(x0), and thus, a concave underestimate L¯e(x) is given by L¯e(x) = H+ ¯e(x0) + (2F¯e(x0)−1)(x−x0)−H− ¯e(x) Due to the concavity of L¯e, it is enough to evaluate this function at the extreme points of ¯eto obtain its minimum on the segment, which is also a lower bound for H+ ¯e. Hence, LB(¯e) = min{L¯e(ai), L¯e(aj)} 5 Computational results The effectiveness of the proposed algorithm was investigated with the aid of numerical cases. The algorithm described in Section 4 was coded in Fortran and compiled using Intel c Fortran Compiler XE 12.0. Executions were carried out on an Intel Core i7 computer with 8.00 Gb of 7 RAM memory at 2.8 Ghz, running Windows 7. The solutions were found to a relative accuracy of 10−3and the integrals were calculated by means of the functions qdags and qdagp available at the IMSL Fortran Numerical Library. We experimented with a set of 43 test networks obtained from [9, 24]. The number of nodes of these test problems ranges from 150 to 1000, and the number of edges from 296 to 3083. Each problem was solved 10 times over each network using randomly generated parameters: the demands of nodes were obtained from a Uniform distribution on [0,1], as it is also the case of the overall demand of each edge. Regarding the demand along each edge, it was assumed to be distributed following a Beta distribution with parameters randomly generated on the interval [0.1,5], which provides a wide range of density functions with very different shapes. After the resolution of each set of 10 instances, some statistical measures (minimum, maximum, average and standard deviation) were calculated for the following indicators of the algorithm performance: •Number of iterations of the Phase 2 of the algorithm. •Maximum size of the data structure used for storage reached during the algorithm execution. •CPU time. Table 1 shows the computational results, where the number of nodes |A|and edges |E|of the graphs are reported as well as the above-mentioned computational measures. The number of iterations and the maximum size of the branch-and-bound list remain low in all the executions. However, the CPU time results are quite high mainly due to the computation of the integrals in Phase 1 (steps 5-b and 5-c) and Phase 2 (steps 4 and 5). One can see that the CPU time required to solve different instances of the same problem shows a great stability. We summarize the findings of this paper. We have addressed the problem of locating one facility on a network with demand distributed on nodes and edges, following arbitrary distributions. The problem is shown to be multimodal, calling for the use of global optimization techniques. The objective function on each edge has been shown to be DC, and a DC decomposition is given. This enables us to obtain concave underestimates of the objective, which are used in a branch and bound procedure. Our numerical tests show that problems of large size are solved rather quickly. Extensions of our techniques to the multifacility case are now under study. References [1] L. Bello, R. Blanquero, E. Carrizosa, “On minimax-regret Huff location models”. Computers & Operations Research 38 (2011) 90–97. [2] D.P. Bertsekas, Network Optimization: Continuous and Discrete Models. Athenas Scientific, Belmont, Mass. (1998). [3] R. Blanquero, E. Carrizosa, “Optimization of the norm of a vector-valued DC function and applications”. Journal of Optimization Theory and Applications 107 (2000) 245–260. 8 Network |A| |E|iterations B&B list time min max mean±std min max mean±std min max mean±std KROB200G 200 386 6 18 11.60±4.03 11 18 13.60±2.46 27.89 32.76 30.14±1.63 KROB150G 150 296 13 21 17.20±2.70 4 4 4.00±0.00 17.61 20.89 19.19±1.17 KROA200G 200 392 11 17 13.70±1.89 8 13 9.60±1.58 30.58 34.60 32.15±1.37 UR137 980 1744 10 16 12.20±2.04 11 11 11.00±0.00 494.23 543.32 511.20±13.80 KROA150G 150 297 11 22 16.90±4.09 6 16 13.20±3.16 17.46 20.01 18.55±0.80 PR152G 152 296 2 7 4.00±1.56 5 6 5.10±0.32 14.77 18.47 16.34±1.03 RAT195G 195 336 4 11 8.10±2.18 3 6 3.30±0.95 17.89 20.31 18.67±0.89 TS225G 225 306 9 15 12.10±1.73 9 11 9.40±0.70 8.94 10.30 9.85±0.46 UR532 298 597 4 44 12.00±13.26 3 28 10.10±7.78 69.69 82.91 73.39±4.14 UR542 343 862 9 17 13.90±2.77 8 18 11.90±3.48 169.28 182.12 175.00±3.98 UR552 388 1135 15 30 22.90±4.72 14 14 14.00±0.00 333.14 347.49 339.86±5.33 UR562 416 1403 23 37 28.20±3.97 27 30 28.60±0.84 532.68 569.67 551.77±12.75 UR732 452 915 1 15 7.90±5.00 6 11 8.00±2.21 164.72 179.26 171.37±4.72 UR535 458 812 0 4 1.60±1.78 3 3 3.00±0.00 102.79 121.29 113.57±5.78 UR545 476 1104 7 20 14.00±3.80 15 24 19.60±2.67 264.00 285.97 275.21±9.58 UR555 490 1305 8 22 16.20±4.59 25 30 26.90±1.79 401.20 429.67 415.26±10.39 UR537 493 868 8 14 11.00±2.05 10 13 12.10±1.45 118.01 134.64 124.53±5.49 UR565 496 1513 15 27 21.90±3.60 23 28 26.30±1.64 573.02 631.32 604.81±19.32 UR547 498 1112 11 19 15.00±2.62 9 14 11.40±1.35 259.46 278.79 269.54±6.59 UR557 498 1310 9 21 16.30±4.16 10 17 11.50±2.42 404.76 425.90 413.57±7.60 UR567 499 1426 13 20 16.10±2.02 5 9 6.80±1.99 510.09 540.01 519.53±8.27 UR742 538 1325 9 20 15.10±4.09 15 18 16.90±0.74 393.67 425.37 406.78±10.22 UR752 580 1735 30 54 43.20±6.68 25 42 29.50±5.60 787.07 825.64 799.86±12.27 UR762 593 2089 31 47 39.40±5.13 21 22 21.80±0.42 1179.76 1271.95 1229.23±26.58 UR132 605 1122 12 27 21.40±3.92 8 8 8.00±0.00 224.33 239.21 233.02±5.58 UR735 662 1200 10 15 12.00±1.56 8 11 10.50±1.08 234.36 255.16 247.08±6.57 UR142 709 1815 5 12 8.40±2.32 4 4 4.00±0.00 748.35 805.64 779.01±17.74 UR745 713 1616 6 16 10.80±3.74 13 14 13.60±0.52 536.82 591.37 561.00±15.58 UR755 724 1966 12 21 15.70±3.27 19 23 20.50±1.84 926.85 972.39 945.26±13.39 UR765 741 2278 6 18 11.90±3.93 13 18 14.20±1.40 1316.04 1366.79 1341.44±16.73 UR737 744 1315 6 14 10.90±2.42 4 4 4.00±0.00 273.45 320.83 290.30±12.75 UR747 745 1659 9 20 14.10±3.31 11 13 11.90±0.57 574.32 609.25 586.98±12.07 UR757 748 1969 12 23 16.70±3.30 19 21 19.90±0.74 891.87 967.30 929.94±21.33 UR767 749 2314 14 22 17.00±2.58 24 53 36.60±9.81 1369.49 1427.39 1390.21±17.92 UR152 766 2390 17 33 24.70±4.74 16 18 17.10±0.88 1455.80 1564.55 1491.42±30.88 UR162 802 2897 30 47 37.80±5.33 23 34 31.00±2.98 2321.47 2419.84 2366.32±32.81 UR135 892 1619 0 7 3.40±2.50 4 6 5.50±0.85 430.73 465.54 447.52±11.17 UR145 929 2117 12 25 20.40±4.03 16 36 28.90±5.72 953.42 1011.92 975.20±18.78 UR155 975 2680 8 22 14.90±3.78 10 10 10.00±0.00 1701.16 1816.37 1755.81±32.00 UR165 980 3068 12 18 15.10±1.97 20 29 24.10±2.47 2368.22 2453.96 2425.75±27.64 UR147 996 2254 10 21 14.60±3.20 5 9 6.90±1.66 1037.24 1113.22 1081.61±23.26 UR157 1000 2690 11 26 18.80±5.35 34 38 36.70±1.42 1671.54 1779.63 1726.31±30.29 UR167 1000 3083 14 25 18.20±3.26 17 34 25.70±7.15 2408.45 2495.42 2451.24±31.00 Table 1: Computational results [4] R. Blanquero, E. Carrizosa, “Continuous location problems and Big Triangle Small Triangle: Constructing better bounds”. Journal of Global Optimization 45 (2009) 389–402. [5] R. Blanquero, E. Carrizosa, P. Hansen, “Locating objects in the plane using Global Optimization techniques”. Mathematics of Operations Research 34 (2009) 837–858. [6] M.L. Brandeau, S.S. Chiu, R. Batta, “Locating the two-median of a tree network with continuous link demands”. Annals of Operations Research 6(1986) 223–253. [7] E. Carrizosa, E. Conde. M. Mu˜noz-M´arquez, J. Puerto, “The generalized Weber problem with expected distances”. RAIRO. Recherche op´erationnelle 29 (1995) 35–57. [8] T.M. Cavalier, T.M., H.D. Sherali, “Network location problem with continuous link demands: p-medians on a chain and 2-medians on a tree”. European Journal of Operational Research 23 (1986) 246–255. [9] A. Corber´an, J.M. Sanchis, “A branch & cut algorithm for the windy general routing problem and special cases”. Networks,49 (2007) 245–257. [10] R. Horst, N.V. Thoai, “DC programming: Overview”. Journal of Optimization Theory and Applications 103 (1999) 1–43. 9