Optimal location problems with routing cost

In the paper a model problem for the location of a given number $N$ of points in a given region $\Omega$ and with a given resources density $\rho(x)$ is considered. The main difference between the usual location problems and the present one is that in addition to the location cost an extra {\it routing cost} is considered, that takes into account the fact that the resources have to travel between the locations on a point-to-point basis. The limit problem as $N\to\infty$ is characterized and some applications to airfreight systems are shown.


The optimal location problem
The location problem consists in determining, given a region Ω ⊂ R d and a nonnegative function ρ : Ω → R describing the density of resources in Ω, the position of a given number N of points x 1 , . . . , x N in Ω in order to minimize the work necessary to concentrate the resources in the points x i . Assuming that the work to move a unit mass from x to y is proportional to |x − y| p , the optimal location problem can be written as the minimization problem (1.1).
When the number N of points is large, as explained in the Introduction, the numerical computation of the optimal points x i is heavy. Therefore an asymptotic analysis as N → +∞, which provides, instead of the precise location of the points x i , only their asymptotic density in Ω, can be very helpful and with much lighter computational costs.
In the case under consideration Ω is a geographic region in R 2 on which a density of resources ρ is distributed; we assume that the density ρ is known. We want to locate in Ω a given number N of airports in the most efficient way according to a global cost that we are going to define. In this section we take into account only the location cost, while the routing cost will be considered in Section 3. The location cost consists in evaluating the work necessary to concentrate the resources distributed on Ω into the airports; we denote by x i the positions of the airports and we assume that the work to move a unit mass from a point x to a point y is proportional to |x − y| p .
Let us denote by m i the quantity of resources that will be concentrated at the point x i and by Ω i the so-called Voronoi cell corresponding to x i , that is the subregion of Ω that sends its resources to the point x i . In other words, we have Then we have that the total cost to concentrate the resources spread on Ω i into the airport 2 x i is given by where A is a proportionality constant. Summing up over all the N airports we have that the total cost is given by that can be also written in the form where Σ is the unknown set of N points to be determined. The most efficient choice of the positions of the airports, considering only the location cost, is then obtained by solving the minimization problem (1.1). When the number N tends to +∞, instead of looking at the precise positions x i in Ω of the airports, one will simply target at determining the limit density µ of the points x i . In order to do it, we identify each set Σ ⊂ Ω of N points with the measure If we assume that, up to a normalization, the density ρ has a unitary total mass, the location cost (2.1) is proportional to the p-th power of the Wasserstein distance between the probabilities ρ dx and µ N . The asymptotic analysis of the cost above has been performed (see for instance [1,5] and references therein) and we summarize here below the available results that, to be correctly stated, require the use of the Γ-convergence, a variational theory developed by De Giorgi and his school starting from the seventies.
When N → +∞ the cost (2.1) is asymptotically equivalent to the limit cost expressed in terms of the limit density µ of points, where C p,d is a constant depending on the exponent p and on the dimension d. It has to be noticed that in the integral above only the absolutely continuous part of µ has to be taken into account, neglecting the singular part. In the case d = 2 the constant C p,2 can be explicitly computed and we have where E is the regular hexagon of unitary area centered at the origin. For instance one has C 1,2 ∼ 0.377 and C 2,2 ∼ 0.16. A plot of the value of C p,2 for p ∈ [0, 2] is given in Figure 1. On the other hand, if we are interested not only in the location x i of the i-th airport but also in the mass m i that is there concentrated, instead of the measures µ N above we have to consider the measures and the optimization problem is written in terms of the p-Wasserstein distance as We notice that, without the normalization ρ dx = 1, passing to the probability ρ(x) ρ dx, the optimization problems (2.2) and (2.3) remain of the same form.

Routing costs
In this section we consider some extra cost that occurs in the model of airports location. Once the resources, spread on Ω with density ρ, have been concentrated at the points x i , with mass m i each, the goods will travel among airports that are connected on a point-to-point basis. We consider two situations: in the first one the routing cost depends linearly on the transported mass, while in the second the connection cost between x i and x j only depends on the distance |x i − x j |. In the following we denote by V (x) the function |x| q .

Mass dependent routing costs
In this subsection we assume that the mass m i concentrated at the point x i is dispatched to the remaining points x j proportionally to the masses m j ; moreover, we assume that the cost to move a unit mass from a point x to a point y is proportional to |x − y| q for a suitable power q. Therefore, the cost to move the entire mass m i is where B is a proportionality constant and m = j m j = ρ dx. Finally, the total routing cost is If we write the routing cost in terms of the measure ν N we obtain and the total cost taking into account location and routing terms gives the optimization problem The characterization of the limit problem as N → ∞ in this case is easy and we can write it as where the minimization above is intended in the class of all measures ν having the same total mass as ρ.
The necessary conditions of optimality for the optimization problem (3.2) can be obtained by differentiating the Wasserstein distance term (see [4]) and the routing cost; we obtain where φ is the Kantorovich potential for the transport from ρ to ν and c is a constant playing the role of the Lagrange multiplier of the mass constraint on ν. In (3.3) the measure ν appears in a very implicit way and can be determined only numerically. One connection between the Kantorovich pontential φ and the transport map T from ρ to ν is given by the Monge-Ampère equation ρ = ν(T ) det(∇T ).

Differentiating in (3.3) we obtain
A∇φ + 2B m ∇V * ν = 0 and T (x) = x − ∇φ(x). Therefore we have the system In dimension 1 we can proceed by an iterative scheme, fixing an initial ν 0 and obtaining T 0 from the first equation in (3.4). Then we can recover ν 1 by the second equation and, assuming T 0 invertible, .
We can now proceed by iterating the scheme above.
Example 3.1. In this particular example we can find an explicit solution taking n = 1, p = 2, and V (s) = |s| 2 . If we suppose that the barycenter of ν is in the origin, we obtain: Putting this in the 1-dimensional Monge-Ampère equation and indicating by v the density of ν, we obtain and changing variables,

Mass independent routing costs
In this subsection we assume that the cost to connect the airports located at the points x i and x j does not depend on the transported mass and amounts simply to K|x i − x j | q where now the constant K is the cost of flying along a unit distance. In this case it is more convenient to use the probability measures µ N introduced in Section 2 which provide the routing cost in the form Taking into account the asymptotic expression of the location cost given in (2.2), we obtain the optimization problem where now µ runs in the class of all probabilities on Ω. Setting ε = A C p,d N −2−p/d /K we are now faced with the problem The necessary conditions of optimality for problem (3.5) simply follow by differentiation of the cost functional and give: where * denotes the convolution operator and c is a constant coming from the mass constraint on µ.
When ε → 0 the optimal densities µ ε of problem (3.5) tend to a Dirac mass δ x 0 for a suitable point x 0 . In order to identify the limit problem as ε → 0, and so to identify the point x 0 around which the optimal densities µ ε concentrate (it can be seen as the main hub of the airports system), it is convenient to rescale the cost above dividing it by its minimum value. Considering the measures a simple calculation provides for the minimal cost of problem (3.5) for a suitable constant C. Optimizing with respect to δ the quantity above we obtain Therefore the rescaled functionals become Note that the optimal measures for F ε and for G ε are the same. In order to characterize the asymptotic behavior of the minimizing sequences (µ ε ) we will compute in the next section the Γ-limit of the sequence of functionals (G ε ). The general theory of Γ-convergence (see for instance [8]) then provides the identification of the main hub x 0 around which the measures µ ε tend to concentrate.

The Γ-convergence result
First of all we notice that, due to the presence of the coefficient ε −1/(1+p/d) in front of the routing term, the Γ-limit on a measure µ will be +∞ whenever Ω×Ω V (x − y) d(µ ⊗ µ) = 0. Therefore, we may limit ourselves to analyze only the measures for which the routing term vanishes, i.e. the Dirac masses µ = δ x 0 .
It is convenient to set notice that α + β = 1 and that α = βp/d. We will show that the Γ-limit of the sequence of functionals G ε , computed on the dirac mass δ x 0 and with respect to the weak* convergence of measures, coincides with the functional

The Γ-limsup inequality
In order to obtain a Γ-limsup inequality we have to choose a suitable sequence µ ε δ x 0 and compute the limit of G ε (µ ε ). We take where the function φ will be chosen later. Then µ ε is a probability measure and we have We choose now φ in order to minimize the quantity at the right-hand side. An easy computation gives

The Γ-liminf inequality
In order to conclude that the Γ-limit of the functionals G ε is the functional H it remains to show the Γ-liminf inequality, which amounts to prove that for every sequence µ ε The following lemma will be useful.
Proof. Since µ is singular, it is concentrated on a measurable set S with |S| = 0 and there is a sequence of open sets (A k ) containing S and such that |A k | → 0. Since Ω \ A k are closed sets, we have for every k A diagonal argument achieves the proof.
Take now a generic sequence µ ε δ x 0 , denote by u ε (x) the density of the absolutely continuous part of µ ε with respect to the Lebesgue measure, and let A ε be the open sets provided by Lemma 4.1. Define We have where we used the fact that Ω×Ω V (x − y) d(ν ⊗ ν) ≥ 0 for every measure ν and that µ 2 ε ≥ u ε 1 A c ε dx. Using the Young inequality we finally obtain (4.1). As a conclusion, the Γ-limit computation is achieved and the optimal main hub for the limit location-routing problem of (3.5) is located at the point x 0 which minimizes the quantity Note that this minimization problem for x 0 is of the form of a Torricelli optimal location problem with suitable exponents.

9
In this section we perform some numerical simulations, based on the results of the previous section, that can be applied to real cases. In the first subsection, some 1-dimension and 2-dimension examples will be presented for different routing costs and density functions and varying ε. The last subsection introduces an application of the model to the USA airfreight system in order to compare the results with the current location of US airfreight hubs.
We are interested to find minimizer of the functional representing the sum of location and routing costs and the x 0 minimizing the Γ-limit functional H(δ x 0 ). In fact, one of the properties of Γ-convergence is the convergence of minima, so if H is the Γ-limit of G ε , the limit of minimizers of G ε is a minimizer of H. We will find numerically the optimal µ ε and we observe that for small ε they are close to a Dirac mass at a suitable point x 0 , according to the previous section's result.
In Section 3 we found an optimality condition for the minimizers of F ε . Unfortunately, condition (3.6) does not admit an explicit solution, so we approximate it numerically. More specifically, we approximate the minimizer via an iterative scheme: we start from the uniform distribution with total mass 1 and then we define the iteration term according to the necessary condition: Here the Lagrange multiplier c, according to condition (3.6), is proportional to ε 1/(1+p/d) .

1-D and 2-D examples
In the one-dimensional case, the domain Ω is the interval [−1, +1] that is discretized in order to solve numerically the problem. Consequently, both the functions ρ and µ are expressed through an array of values in correspondence of the discretization points. At each step of the convergence procedure shown in Figure 2, µ n+1 is obtained from the relationship (5.1) and then normalized to a probability measure. The first simulation, reported in Figure 3, is related to a non-symmetric distribution of population density ρ and a quadratic routing cost function: We assume that the convergence is reached when the maximum error between the values of µ n and µ n+1 is less than 2%. In this conditions, about 10 iterations are requested to solve the problem and the computational time results to be proportional to the number of point used to discretize the domain; when it is divided into 200 steps, the calculation time is about 100 sec.
The results for different values of ε coefficient show that, as the ε decreases, the limit density µ tends to have a concentration centred on a single point. At the limit as ε → 0, the density µ becomes a Dirac mass located at the point x 0 that minimizes the functional (4.2). 10 Figure 2: numeric procedure for the determination of probability distribution µ The convergence towards the limit conditions of ε → 0 is slow and it cannot be reached numerically because the onset of numerical errors below the value of ε 10 −4 . Therefore, the routine can be completed by calculating also the value of the functional H(δ x 0 ) reported in (4.2), and find the point x 0 of minimum. In this case, the minimum of the functional (4.2) can be found explicitly: which gives The analytical solution equals to the value determined by the numerical procedure that is also reported in Figure 3.
The population ρ often can have an uneven distribution among the domain, and therefore an adequate function is requested in order to model correctly this aspect. A first solution can be provided by treating this distribution as a sum of M Gaussian functions: where the coefficients A j ,B j ,X j are used to set respectively the height, the width and the position of the j-th peak. The results reported in Figure 4 refer to the case of a population ρ with 8 peaks of different position, height, and area of influence (width); also in this case the simulations have been conducted with two different values of the coefficient ε.
The density of probability µ is highly dependent by the values of the coefficients; the point x 0 is in this case numerically determined: x 0 +0.044.
As the ε decrease, the influence of the routing cost becomes larger on despite of the location ones so that the system tends to minimize the airport distance. The large differences between the solutions remarks the importance in choosing a value of the coefficient as realistic as possible. Moreover, one can note that the computational time is not affected by the "complication level" of the population function but only by the used discretization step. The effective decisional process related to the facilities (the airports) location, can be done in a post-processing phase: in this way, we can decide how many airports can be located in a given region, proportionally to the area limited by the density distribution; for example, the numbers on the X-axis of Figure 4 equals to the airport on each step (each step length is 0.2).
The routine has been applied also in the 2-D case, considering a correspondent peaks distribution shown in Figure 5. We remember the Gaussian function in the case of two variables and also both the routing cost function and the exponent of location: V (x) = |x| 0.5 , p = 1. Although the calculation procedure does not change, the computational time results much higher than in the previous cases because of the great number of points requested to discretize properly the domain. Neverthless, it remains notably lower than the ones of the common Operating Research models (about 2200 sec. when the domain is divided into 1600 cells).
The result shows that the probability density follows the shape of the initial population ρ(x, y) (we can note that the exponent q=0.5 determines the minor importance of the routing costs on despite of location ones) the point of maxima can be observed near the central peaks where the effect of both the location and routing terms are summed.
When the population distribution becomes irregular the position x 0 of the "main hub" cannot be estimated immediately but the functional (4.2) can be easily computed and its minimum can be founded. For the same 2-D case, the values of the functional are shown in Figure 7 while its minimum point is depicted in Figure 7 together with the level curves of the population ρ.

Application to the US airfreight system
Two main problems will be faced in order to apply location-routing models to real cases: • Location and Routing terms are related with ground and air transportation costs respectively. Preliminarily, we can suppose a linear dependence of ground cost with the transport distance but the same assumption becomes not valid in the case of air transportation.
• The distribution of population ρ identifies the airfreight demand among the domain; data are directly available only for some areas (occidental countries) while in most cases an extrapolation from some socio-economic data is needed. 13 Therefore, the aim of the present section is to set the coefficients and exponents appearing in (3.5) in such the way the terms of cost functional reflect as realistic as possible the dynamics of the real world. The numerical routine will be finally applied on the US domain that will be represented as a polygon on a Cartesian system.

Routing and airfreight cost
In the functional (4.2), the routing costs can be modelled as a function of the transport distance, through the general power relationship V (x) = K|x| q by simply setting the coefficient K and the exponent q. In the case of airfreight, most of cost terms depend strongly by the economies of scale in which the company operates (countries connected, commercial accordances, kind of service done, aircraft used) so that it often is not possible to determine a general function that could be valid in every case. Nevertheless, if we cannot determine an explicit function, we can determine its "shape" by supposing that the part of cost variable with the transport distance, is mainly related to the fuel consumption during flight. In a first approximation, the amount of fuel required for a given mission can be determined (for example in the case of constant power aircraft) by using the so-called Breguet relations: where the Range equals to the transport distance, E is the aerodynamic efficiency of the considered aircraft, and k c and η p are respectively the Specific Fuel Consumption and the propeller efficiency: since these parameters are all known for each aircraft and engine, the amount of fuel and its cost can be calculated in dependence of the Range flown. The operating costs are usually reported in terms of Costs per Unit of carried mass and flown distance, Cost/(T on · Km)(simply by dividing by the total payload and the transport distance) and the results of this procedure deriving from Breguet, has been compared in Figure 8 with some statistical models ( [11] and [12]) that use a regression of both historical data about existing freighter and data collection of the financial report of transport companies. The Curves in Figure 8 have a similar shape and they differ only by a translating coefficient that can be related to the different economies of scale which data are extrapolated 15 from. Moreover, the Cost/T on can be determined by integrating the curves in Figure 8 so that finally a suitable value of the exponent q is determined: The value results lower than 1 because, conceptually, the air cannot be considered as a constant mass transport: in the case of existing aircraft in fact the Weight of embarked fuel represents until 30-40% of the total mass so that the flight condition and, consequently, also the fuel vary notably during cruise (as the aircraft is lightening, the burn fuel decreases).

Modeling the airfreight demand
The airfreight demand is often not directly measurable so that also the initial population ρ has to be properly modelled. The identification of the socio-economic parameters affecting the airfreight demand, results very difficult; also in this case the models are highly affected by the economies of scale and the geographical region on which the air transport is operated. In the present study, we assume that the airfreight depends by some socio-economic parameters in a likely linear regression as proposed in [13]. In the model proposed, the airfreigth demand in some point of the domain is obtained by the following expression: • C 0 , .., C 5 : coefficient coming from linear regression of economic data • AF : volume of airfreight demand (T ON ) • P C: per capita personal income ($1,000) • T SE: traffic shadow effect. In first approximation, this parameter will not considered in order to avoid any iterative process also for the input data.
• T SL: transportation-shipping-logistics employment market share (%) • M D:# of medical diagnostic establishments • HT : average high-tech employee wage ($1,000) An airport has a relatively small catchment region (cities or districts), so that on the "ground side", the airfreight demand has influence on a very small area on despite of the worldwide dimension of the air transport; for this reason the airfreight function must to be refined also if the domain is very large and, consequently, the required discretization step is small. In this contest,the socio-economic data in (5.5) are extrapolated from common statistical reports (National Bureau for USA or Eurostat for E.U.) for each metropolitan districts, so that a detailed function can be easily determined.
The points used to define the spatial distribution of airfreight demand, are reported in Figure 9 in blue dots, and their position refers to the centroids of the so called metropolitan statistical area of the US territory. 16  Since data are know in correspondence of these points, the ρ function is then obtained through a cubic interpolation with a matlab routine. In Figure 10 is reported the ρ obtained by this procedure: it has an uneven distribution, peaks are concentrated in very rich or very populated regions and their area of influence is relatively small.

Results
The Figure 11 shows the level curves of the µ function as result of the real case study in which the coefficient ε is set ε 10 −1 .
The black dots of Figure 11 indicate the positions of the 10 major cargo airports in US. The global maxima of the µ density results very close to the Memphis airport which is the busiest center of airfreight transport and the hub of the FedEx: its airfreight volumes are doubled respect to the other airports. Also the other local maxima are located near the effective position of the other airports: largest errors are appreciated along the boundary areas where also the used projection method presents the largest errors.
The level curves of the functional H(δ x 0 ) are plotted in Figure 12. Although the exponent q results lower than the unity, the effects of the routing costs tend to predominate the location 17 ones and consequently the importance of initial population ρ is reduced on despite of the distance power relationship. The position of the minimum point (the "main hub") differs from the global maxima of the limit density µ displayed in Figure 11. In this case, the difference is due to the value of the coefficient ε used to determine the limit density µ, which is relatively high so that the results of the functional (3.5) and (4.2) are not coincident.

Conclusions
Some considerations can be done. Our initial problem was to determine the optimal position of a certain number N of airports into a domain with location and routing cost condition. This problem is hard if we try to solve it with a direct approach because of his intrinsic complexity. After the modelling phase we concentrate to mass independent routing cost and we caracterized the asymptotic behaviour of the total cost problem. This means that instead of finding the exact position of the N airports, we compute a probability density that represent the "importance" of a certain point in the area taken as domain. Moreover, gamma limit result allow us to find the position of the optimal main hub minimizing functional H(δ x 0 ).
Supported by the examples of the one and two dimensional cases, we observe that these two limit problems are very "easy", in terms of computational costs. So instead of looking at the initial problem is more convenient to be reduced at the other two. This makes possible to apply the procedure described to real cases, as in the USA airfreight system.