FARD Milad Keshvari

Laurent Alfandari

Milad Keshvari Fard

TRADE-OFFS BETWEEN THE STEPWISE COST FUNCTION AND ITS LINEAR APPROXIMATION FOR THE MODULAR HUB LOCATION PROBLEM RESEARCH CENTER

There exist situations where the transportation cost is better estimated as a function of the number of vehicles required for transporting a load, rather than a linear function of the load. This provides a stepwise cost function, which defines the so-called Modular Hub Location Problem (MHLP, or HLP with modular capacities) that has received increasing attention in the last decade.

In this paper, we consider formulations to be solved by exact methods. We show that by choosing a specific generalized linear cost function with slope and intercept depending on problem data, one minimizes the measurement deviation between the two cost functions and obtains solutions close to those found with the stepwise cost function, while avoiding the higher computational complexity of the latter. As a side contribution, we look at the savings induced by using direct shipments in a hub and spoke network, given the better ability of a stepwise cost function to incorporate direct transportation.

Numerical experiments are conducted over benchmark HLP instances of the OR-library.

Introduction

Hubs are facilities that enable to decrease the number of transportation links between origins and destinations in many-to-many distribution systems. They play an important role for several purposes such as consolidation, switching, connection and sorting points in a large number of businesses and industries. The smaller number of links required, as well as the economy of scale provided by consolidation of the loads are the main incentives for a hub and spoke transportation system. The consolidation happens all over the transportation network: on the routes from origins to hubs, from hubs to destinations, and between the hubs. Considering the higher amounts of load between the hubs, larger transportation modes may be used, which results in lower per-unit costs. The telecommunication industry, airline industry, postal delivery, 1 public transit, maritime industry and freight transportation companies are areas that benefit from hubs.

The Hub Location Problem (HLP) consists of locating hub facilities, as well as allocating supply and demand nodes to the hubs in order to route the traffic between origin-destination (O-D) pairs. The objective is to minimize the total transportation cost of the distribution network system. This problem first appeared in [START_REF] Toh | The impact of hub and spoke network centralization and route monopoly on domestic airline profitability[END_REF], where the authors discussed the application of hub location in the airline industry. It was followed then by O' Kelly (1986) and [START_REF] Campbell | Integer programming formulations of discrete hub location problems[END_REF] who presented the first mathematical formulation and first integer programming formulation for the problem. Later, [START_REF] Skorin-Kapov | Tight linear programming relaxations of uncapacitated p-hub median problems[END_REF] proposed a tighter integer programming formulation, which has been widely used in the literature.

Due to the wide applications of the HLP, the problem has been studied under different assumptions and settings: piecewise linear cost function (O'Kelly andBryan 1998, Ishfaq and[START_REF] Ishfaq | Intermodal logistics: The interplay of financial, operational and service issues[END_REF], capacitated hubs Contreras et al. (2012), single versus multiple allocation of the nodes to the hubs (O'Kelly 2012), continuous p-HLP [START_REF] Aykin | Interacting new facilities and location-allocation problems[END_REF], fixed costs for establishing hubs [START_REF] Tanash | An exact algorithm for the modular hub location problem with single assignments[END_REF], hub location under uncertainty [START_REF] Contreras | Stochastic uncapacitated hub location[END_REF][START_REF] Alumur | Hub location under uncertainty[END_REF][START_REF] Meraklı | A capacitated hub location problem under hose demand uncertainty[END_REF], or competition [START_REF] Lüer-Villagra | A competitive hub location and pricing problem[END_REF]Marianov 2013, Mahmutogullari and[START_REF] Mahmutogullari | Hub location under competition[END_REF], multi-period hub location [START_REF] Gelareh | Multi-period hub location problems in transportation[END_REF], etc. Different solution approaches were proposed for the problem as well. One can refer to [START_REF] Alumur | Network hub location problems: The state of the art[END_REF], [START_REF] Campbell | Twenty-five years of hub location research[END_REF], [START_REF] Farahani | Hub location problems: A review of models, classification, solution techniques, and applications[END_REF], [START_REF] Contreras | Hub location problems[END_REF] for an overview of different models, classifications and solution approaches in the HLP literature. Moreover, [START_REF] Contreras | General network design: A unified view of combined location and network design problems[END_REF] provide an overview of hub location models in the context of network design, where some of the classical HLP assumptions such as the full connectivity of the hubs, or direct connection of all terminal nodes to the hubs are relaxed, to attain more efficient or effective transportation networks.

The majority of hub location research assumes a simple function for the transportation cost, with a cost rate (per unit of flow and unit distance) that is discounted on the arcs between hubs, to reflect economies of scale from consolidation on these arcs [START_REF] O'kelly | The location of interacting hub facilities[END_REF][START_REF] Skorin-Kapov | Tight linear programming relaxations of uncapacitated p-hub median problems[END_REF][START_REF] Ebery | The capacitated multiple allocation hub location problem formulations and algorithms[END_REF][START_REF] Cánovas | Solving the uncapacitated multiple allocation hub location problem by means of a dual-ascent technique[END_REF][START_REF] Gelareh | Multi-period hub location problems in transportation[END_REF]. In practice however, in airline of trucking companies, these transportations are usually characterized by the vehicles traveling between the network's nodes [START_REF] O'kelly | Multiple allocation hub location model with fixed arc costs[END_REF], which can lead to using a stepwise cost function associated with the number of vehicles used on a link. Considering instead simple linear cost functions may not be a good reflection of transportation network operations in some cases [START_REF] O'kelly | Multiple allocation hub location model with fixed arc costs[END_REF], and consequently, result in sub optimal decisions. The Hub Location Problem with a stepwise cost function is often called Modular Hub Location Problem (MHLP), or HLP with modular capacities, and has received increasing attention in the literature. Studying the problem of locating hubs in a telecommunication network, [START_REF] Yaman | Solving the hub location problem with modular link capacities[END_REF] considered a stepwise cost function for the required capacity to be installed on the links that connect hubs together, and proposed a heuristic solution for the problem. [START_REF] Tanash | An exact algorithm for the modular hub location problem with single assignments[END_REF] considered a modular HLP in the context of airline industry, and proposed a Branch-and Bound algorithm for the problem. Other papers such as [START_REF] Takano | A genetic algorithm for the hub-and-spoke problem applied to containerized cargo transport[END_REF], [START_REF] Corberán | Strategic oscillation for the capacitated hub location problem with modular links[END_REF] and [START_REF] Hoff | Heuristics for the capacitated modular hub location problem[END_REF] studied similar problems with a stepwise cost function and developed heuristics to address the complex nature of the problem. Moreover, there exist relevant articles such as [START_REF] Jaillet | Airline network design and hub location problems[END_REF], [START_REF] Janic | Modelling the full costs of an intermodal and road freight transport network[END_REF] and [START_REF] Yaman | Star p-hub median problem with modular arc capacities[END_REF] that consider similar cost functions for hub location in the context of network design.

Figure 1 illustrates the difference between a stepwise and a linear cost function. In a stepwise cost function, the transportation cost does not systematically increase as the load increases, but when the number of required vehicles increases (Figure 1a). On the other hand, the two cost functions imply different cost-per-unit load (Figure 1b). While a simple linear function assumes a fixed unit cost, in a stepwise function the cost per unit of load depends on the overall load to be transported; when all vehicles are full, the per-unit load transportation cost is at the lowest, while having a relatively empty vehicle would lead to a higher per-unit transportation cost. In fact, using a simple linear cost function may provide a rather far estimation of the transportation cost on a link. Such a difference in cost calculation may result in different location/ allocation decisions than those which would be selected if a stepwise cost function were used.

However, the MILP formulations with stepwise cost function contain a larger number of integer variables, which makes it harder to use for larger problems. It is not surprising then that most of the papers which study the MHLP come up with heuristic solutions.

To overcome this issue, we seek to approximate the stepwise cost function by a generalized linear function with a non-zero intercept, which could refer to the fixed (or sunk) cost of vehicle(s) required for the transportation. Although there exist a number of papers in the HLP literature that consider such a linear cost function [START_REF] Adler | Hub-and-spoke network alliances and mergers: Price-location competition in the airline industry[END_REF][START_REF] O'kelly | Fuel burn and environmental implications of airline hub networks[END_REF][START_REF] O'kelly | Multiple allocation hub location model with fixed arc costs[END_REF], none of these research define their linear functions as approximations of stepwise functions, nor provide analyses for that purpose.

We determine the slope and intercept of the linear function so that the expected total transportation cost is the same for both stepwise and linear functions. Our results show that the optimal decisions found under these two cost functions are most often very close.

Indeed, comparing two HLP formulations, one with stepwise and another with a linear cost function, our findings indicate a strong correlation between the difference in the total transportation cost computed by the two formulations on one hand, and their difference in their optimal location/ allocation decisions on the other.

In order to develop our analysis, we focus on one of the most common settings in the HLP literature:

the Uncapacitated Single Allocation HLP, where hubs are fully connected and each node is exclusively allocated to one hub. Fully inter-connected hubs are common in practice and arise in a broad range of applications [START_REF] Contreras | General network design: A unified view of combined location and network design problems[END_REF]. Moreover, restricting the non-hub nodes to be connected to exactly one hub is economically justified in several applications, due to amalgamation of flows into efficient bundles [START_REF] Skorin-Kapov | Tight linear programming relaxations of uncapacitated p-hub median problems[END_REF].

We also consider the possibility of direct shipments from origins to destinations in a hub and spoke network with modular capacities, since a stepwise or a generalized linear cost function provide a better ground to model this feature. Despite its extensive applications in real world problems, only a small fraction of the HLP literature considers the possibility of direct shipments in the network [START_REF] Steadieseifi | Multimodal freight transportation planning: A literature review[END_REF].

Indeed, there are several occasions where it is more economic to send the flow directly from the origin to the destination. Other concerns such as the preference of passengers for direct flights in the airline industry [START_REF] Aykin | Lagrangian relaxation based approaches to capacitated hub-and-spoke network design problem[END_REF], or benefiting from the existing train lines in railroad transportation [START_REF] Racunica | Optimal location of intermodal freight hubs[END_REF], make the direct shipment inevitable for many applications. Our results indicate a significant reduction in transportation cost by incorporating direct shipments in the MHLP, and this cost reduction is decreasing in the shipment frequency throughout the supply chain.

The paper is structured as follows. Section 2 provides the MILP formulations for the MHLP as well as a generalized linear cost function. Section 3 compares the stepwise and the generalized linear cost functions, and provides linear approximations of the stepwise cost function, as well as an approximation of the deviation between the total cost, computed by a stepwise and a linear function. Section 4 extends the formulations to incorporate the possibility of direct shipments. The computational results and managerial insights are reported in Section 5, together with managerial insights. Finally, Section 6 concludes the paper.

MILP formulations

Consider a set V of demand nodes, where each node is a potential location for setting up a hub. We assume an Uncapacitated Single Allocation HLP, where each node is allocated to only one hub, and all hubs are connected. For simplicity, hereinafter we refer to links that connect the terminal nodes to their hubs as c/d links (collection/distribution), and the links that connect two hubs as h/h links. We assume c/d links use small vehicles of capacity Q 2 , and a unit cost per load and distance of C 2 , while h/h links use larger

vehicles of capacity Q 1 > Q 2 and unit cost C 1 > C 2 .
Economies of scale of inter-hub transportation imply that C1 Q1 < C2 Q2 . Assuming modular capacities over the links, the transportation cost does not increase linearly in the load, but rather as a stepwise function of the number of vehicles, as in [START_REF] Jaillet | Airline network design and hub location problems[END_REF], [START_REF] Hoff | Heuristics for the capacitated modular hub location problem[END_REF]. It calculates the transportation cost of a load L over a distance d as

Cost stp (L) = C L Q d, (1)
where Q is the capacity of the vehicle, C represents the cost per distance unit of that vehicle (possibly including fuel, maintenance, insurance, driver costs, etc.) and x denotes the smallest integer larger than or equal to x. An HLP formulation with a stepwise cost function encourages solutions with higher load aggregation and is sparing in the use of linkages and vehicles.

Let L ij denote the load to be transported from node i to node j, and d ij the distance between nodes i and j. We note

V c i = j Lij Q2
the number of vehicles required to collect the load from node i to its related hub, and

V d j = i Lij Q2
the number of vehicles required to distribute the loads to node j from its hub.

The decision variables for the Modular Hub Location Problem (MHLP) are:

• Y ik = 1 if node i is allocated to the hub at location k, and 0 otherwise;

• Z ijkm = 1 if the load shipped from node i to j transits through hubs k and m, and 0 otherwise.

• V

h/h km : the number of large vehicles used to transport the load from hub k to hub m.

The Modular Hub Location Problem (MHLP) under the uncapacitated single-allocation assumption can be written as follows.

min i k C 2 d ik V c i Y ik + k m C 1 d km V h/h km + j m C 2 d mj V d j Y mj (2) s.t. k Y kk = p (3) k Y ik = 1 ∀i (4) Y ik ≤ Y kk ∀i, k (5)
m

Z ijkm = Y ik ∀i, j = i, k (6) k Z ijkm = Y jm ∀i, j = i, m (7) i j L ij Z ijkm ≤ Q 1 V h/h km ∀k, m = k (8) V h/h km ∈ Z + ∀k, m (9)
Z ijkm ≥ 0, Y ik ∈ {0, 1} ∀i, j, k, m (10)
The terms in Equation ( 2) calculate the cost of vehicles used for collection, inter-hub transportation, and distribution respectively. Constraint (3) implies that p hubs have to be located. Constraints (4) state that each node is allocated to only one hub. Constraints (5) imply that to allocate a node i to node k, a hub has to be constructed at location k. Constraints ( 6) and ( 7) relate variables Y to variables Z, and imply that if Z ijkm = 1 then Y ik = Y mj = 1. Constraints (8) compute the number of vehicles on each h/h link as

V h/h km = i j L ij Z ijkm /Q 1 .
Finally, Constraints (10) and (9) determine the type of each variable. Note that a similar formulation was studied in [START_REF] Tanash | An exact algorithm for the modular hub location problem with single assignments[END_REF]. There exist other formulations for the Modular HLP in the literature, such as the formulation with quadratic constraints in [START_REF] Yaman | Solving the hub location problem with modular link capacities[END_REF], [START_REF] Hoff | Heuristics for the capacitated modular hub location problem[END_REF], or the flow-based formulation of [START_REF] Tanash | An exact algorithm for the modular hub location problem with single assignments[END_REF]. However, since the focus of the current research is on the tradeoff between different cost functions, and all formulations provide the same integer solution, we use the current path-based MILP formulation to study the choice between stepwise and linear functions.

Due to the large number of integer variables, the above MHLP formulation is likely to be intractable for medium-size problems. For example, for a similar formulation, [START_REF] Tanash | An exact algorithm for the modular hub location problem with single assignments[END_REF] developed a branch and bound algorithm that solves instances up to size n = 40 whereas, for the classical HLP, instances of size 200 and up to 500 could be solved, using Benders-Decomposition techniques (Contreras et al. 2011a). Since the formulations with linear cost functions allow to solve larger problems to optimality, we approximate the stepwise cost with a linear function. Such a function should allow for a non-zero intercept, to consider the fixed cost of using vehicles over a transportation link, and has been used in the literature [START_REF] Adler | Hub-and-spoke network alliances and mergers: Price-location competition in the airline industry[END_REF][START_REF] O'kelly | Fuel burn and environmental implications of airline hub networks[END_REF][START_REF] O'kelly | Multiple allocation hub location model with fixed arc costs[END_REF]. We refer to this function as Generalized Linear Cost (GLC) function, and define it as:

Cost lin (L) = (aL + b)d (11)
where a and b represent the slope (e.g. cost rate) and the intercept (e.g. fixed cost) of transportation cost function respectively. Based on Equation ( 11) we define the HLP with Generalized Linear Cost function (HLP-GLC) formulation as follows:

min i k (b 2 Y ik + a 2 j m Z ijkm L ij )d ik + j m (b 2 Y jm + a 2 i k Z ijkm L ij )d jm k m (b 1 w km + a 1 i j Z ijkm L ij )d km (12) s.t. w km ≥ Y kk + Y mm -1 ∀k, m (13)
Constraints (3) to ( 7), ( 10)

w km ≥ 0 ∀i, j (14)
where a 1 and b 1 represent the slope and intercept of the linear cost function on h/h links, while a 2 and b 2 identify the cost function on c/d links, and auxiliary variable w km is equal to 1 if nodes k and m are both hubs, and 0 otherwise. The terms in the objective function (12) refer to the collection cost, inter-hub transportation cost and distribution cost, respectively. Constraints (13) determine if the link between nodes k and m is a h/h link. Note that similarly to Z ijkm , variables w km only take binary values in the optimal solution. The HLP-GLC formulation reduces to the [START_REF] Skorin-Kapov | Tight linear programming relaxations of uncapacitated p-hub median problems[END_REF] formulation by setting b 1 = b 2 = 0. Thus, regardless of the linear cost function's intercept, hereinafter we refer to the MILP formulations that use a linear cost function as HLP-GLC.

Both formulations above reduce to the classical HLP formulation of [START_REF] Skorin-Kapov | Tight linear programming relaxations of uncapacitated p-hub median problems[END_REF], by removing Constraints ( 8) and ( 9) from the MHLP formulation, and Constraints ( 13) and ( 14) from the HLP-GLC formulation, as well as replacing the objective function with

min i j k m L ij (a 2 d ik + a 1 d km + a 2 d mj )Z ijkm ( 15
)
with a 1 = αa 2 , where α ≤ 1 is a discount factor to address economies of scale over h/h links.

In the next section, we introduce a way to define the slope a and intercept b of the linear cost function of Equation ( 11), so that the expected measurement deviation between the linear and stepwise cost functions is minimized. Furthermore, we estimate the expected deviation between the two cost functions for any given parameters a and b.

3 Comparison of linear vs stepwise formulations

Notations and probabilistic assumptions

For a given link l with a total flow of L l , we define the measurement deviation between the transportation costs computed by both the stepwise and the linear cost functions as l = Cost stp (l) -Cost lin (l). Now, let Y be an arbitrary solution of the HLP (i.e., location of the hubs and allocation of the terminal nodes to these hubs are given). Using the linear cost function of Equation ( 11 ). We refer to this gap as the solution gap δ between the optimal solutions of the two formulations:

δ = T C stp (Y * lin ) -T C stp (Y * stp ) ≥ 0. ( 17
)
We use δ as an indicator of the sub-optimality of the HLP-GLC's solutions, when compared to those of the MHLP. Due to the higher tractability of the HLP-GLC formulations, one may prefer HLP-GLC over MHLP as long as the solution gap δ is small. Although it is hard to infer some bound for δ in general, our numerical results in Section 5 show that a linear cost function that in average, results in a smaller measurement deviation for any given solution, is more likely to find optimal solutions with smaller solution gaps δ as well. Therefore, for a given stepwise cost we need to find a linear approximation with the minimum measurement deviation.

In case we know the exact load on a link l, we can define a linear cost function specific to l that results in l = 0. Assume a load L to be transported using vehicle(s) of capacity Q and cost per unit distance of C; in this case, such a linear approximation of form (aL + b)d is given by:

a = L Q C L , b = 0. ( 18
)
In a modular HLP, the above cost function is applicable for the c/d links of the HLP-GLC model presented in Section 2, since the total load to be transported from a node to its associated hub is known. However, this approach is not applicable to the h/h links, as the total flow on a link depends on the allocation of terminal nodes to the hubs. Moreover in some applications, a unique cost function for all links of the same type (e.g. all c/d links) may be preferred.

Therefore, we need to define a generalized linear function for a family of links that share the same stepwise function, e.g. c/d links or h/h links. Since the flow on each link depends on factors such as the load matrix and location/allocation decisions, and is not known before solving the problem, we design a probabilistic study to find a linear approximation of the stepwise function that in average, has a zero measurement deviation for that family of links, i.e. E l ( l ) = 0 .

We assume that for all O-D pairs (i, j), the load L ij is a random variable with expected value of μ = E[L ij ] and variance σ2 . For the numerical purposes, these parameters can be calculated from the load matrix.

Let random variable L l be the total flow passing over link l. This implies that L l is the sum of a number n l of loads L ij that pass through link l. According to the Central Limit Theorem (CLT), as n l increases the distribution of L l tends to a Normal distribution; the speed of convergence depends on the distribution of L ij s. For example, according to [START_REF] Renyi | Probability theory[END_REF], the sum of n random variables with identical uniform distributions follows an Irvin-Hall distribution, that is very close to a Normal distribution even for small values of n, typically n ≥ 4 [START_REF] Killmann | A note on the convolution of the uniform and related distributions and their use in quality control[END_REF]. Indeed, even for relatively small or medium networks with tens of demand points, the load on each c/d or h/h link would consist of several O-D loads, which is enough to satisfy the CLT, regardless the distribution of L ij . Hence in the sequel, we approximate the distribution of L l , f (L l ) by a Normal distribution, i.e. f (L l ) ∼ N (µ l , σ 2 l ).

Best approximation of the stepwise cost function by a generalized linear function

Theorem 1 specifies the best linear approximation of a stepwise cost function. The resulting linear function could be used for a family/group of similar links such as h/h or c/d links.

Theorem 1. The best linear approximation of a stepwise cost function with parameters C and

Q that satisfies E[Cost stp (L l )] -E[Cost lin (L l )] = 0 is given by Cost lin (L l ) = C Q L l + C 1 2 + δ l 2Q -∆ l , (19)
where δ l = Q 2µ l Q -2µ l , with µ l as the expected load on link l and

∆ l = µ l Q t=1 tQ tQ-δ l f (L l )dL l + µ l min{µ l , µ l Q Q-δ l } f (L l )dL l .
All proofs can be found in Appendix A. Note that for a generic link l, the slope C Q does not depend on the load's density function f ; however, this holds true for the intercept, only if µ l is a multiple of Q 2 .

Corollary 1. If µ l is a multiple of Q/2, then Cost lin (L l ) = C Q L l + C 2 .
In case µ l is not a multiple of Q 2 , we need to estimate the load's distribution, f (L l ). Decomposing the analysis into c/d and h/h links, we have:

(i) For c/d links, the load distribution is f (L l ) ∼ N (n -1)μ, (n -1)σ 2 , ( 20
)
since the load carried on a c/d link from terminal node i to its hub is the sum of loads from/to node i to/from n -1 other nodes. Thus using the CLT, the load distribution for a randomly selected c/d link follows Equation ( 20), which enables to compute b 2 using Equation ( 19).

(ii) For h/h links, we can estimate µ l and σ l assuming cluster symmetry, i.e., each hub serves n/p customers (if n/p is fractional, we keep the fractional value for approximation), which enables to consider the same expected load on all h/h links. Therefore according to the CLT, the load distribution on h/h links is

f (L l ) ∼ N n 2 p 2 μ, n 2 p 2 σ2 . ( 21
)
The numerical results is Section 5 confirm the robustness of this assumption.

Nonetheless, even without any knowledge of the load distribution, one can simply approximate the intercept of the linear cost function, only based on the cost and capacity of the vehicle being used as follows.

Corollary 2. When 2µ l Q is large, the GLC function proposed by Theorem 1 can be approximated by a GLC with slope a = C Q and intercept b = C 2 .

Predicting the measurement deviation when using an arbitrary linear cost function

Using the linear cost functions identified by Theorem 1, results in an expected measurement deviation equal to zero. However, this may not hold true if other linear functions are being used. Since our results in Section 5.1 indicate a positive correlation between the measurement deviation (in absolute value) and the solution gap, it would be of our interest to approximate the expected measurement deviation between the stepwise cost and any given linear function. To have a robust measure that is independent of the problem size or load magnitude, we define the relative measurement deviation ˜ (Y ) as

˜ (Y ) = (Y ) T C lin (Y )
Proposition 1 provides an approximation of the expected relative measurement deviation, ˆ , for a hub and spoke network.

Proposition 1. An approximation ˆ of the expected relative measurement deviation E[˜ (Y )] when using a linear cost function with slope a t and intercept b t as an approximation of a stepwise function with cost C t and capacity Q t , assuming symmetry in allocation of terminal nodes to the hubs, is given by

ˆ = A 1 ρp(p -1) + 2A 2 (n -p) b 1 ρp(p -1) + a 1 ρn 2 p-1 p μ + 2b 2 (n -p) + 2a 2 (n -p)(n -1)μ ( 22
)
where

A t = C t E L l Qt -a t µ l -b t , t = 1, 2
, and ρ is the ratio of the average length of a h/h link by that of a c/d link.

For Equation ( 22), all parameters are known a priori except ρ. In Appendix B we propose a heuristic to estimate ρ based on the number of hubs p. Moreover, A 1 and A 2 can be computed by numerical methods, using the approximations of the load distributions given by Equations ( 20) and ( 21). Although we assume a symmetry in allocation of nodes to the hubs to obtain this result, our results in Section 5.2 show that the prediction of Proposition 1 is robust, regardless of the allocation of nodes to the hubs or distribution of the load.

Note that any stepwise cost function of a more general form, such as

Cost stp = V + (R + C L Q )d
, where V and R are constants (see [START_REF] Swan | Aircraft trip cost parameters: A function of stage length and seat capacity[END_REF]), can be approximated by a linear function of the form Cost lin = v + (aL + b)d, using Theorem 1 with small modifications. Similarly, Proposition 1 is applicable to such cost functions with few changes. However, since such cost functions are less common in the HLP literature, we do not develop this variant in this paper.

Extension of formulations to incorporate direct shipments

As mentioned earlier, using direct connections in parallel with through-hub shipments is common in practice. In this section we incorporate direct shipments into our MILP formulations, MHLP and HLP-GLC.

The incentive is the better ability of a stepwise or a generalized linear function, due to enforcing fixed transportation costs, to model direct connections in a hub and spoke network.

With a simple linear cost function (with zero intercept), the transportation cost for direct connections can be modeled with a higher cost-per-unit rate, e.g. c for direct links vs. c for c/d links, with c ≥ c. This implies that regardless of the load magnitude, for any O-D pair (i, j) that c d ij < c(d ik +αd km +d mj ), direct shipment is preferred over transportation through hubs k and m. Consequently, if c is not considerably larger than c, one could expect a significant proportion of shipments -especially between two non-hub nodes that belong to the same cluster-to be transported directly. Such a solution would require a larger number of vehicles, with a lower capacity utilization rate and a higher total transportation cost in reality. This issue, however, is easily solved using a stepwise cost function, since for having a direct shipment the full capacity cost of required vehicle(s) for that link has to be paid.

To incorporate direct shipments in the MHLP formulation presented in Section 2, we define binary variables q ij which is equal to 1 if j receives its demand from i directly, and 0 otherwise. Without loss of generality, we assume C 3 and Q 3 represent the cost and capacity of the vehicles used for direct shipments.

Finally, we define

V s ij = Lij Q3
as the number of vehicles required for transporting the load directly from i to j. Note that unlike Section 2, the vehicles required for collection V c ik and distribution V d mj become decision variables, since part of the flow on the c/d links can be transported by direct shipments and consequently, less vehicles may be required. The MHLP with direct shipments can be written as 3), ( 4), ( 5), ( 8)

min i k C 2 d ik V c ik + k m C 1 d km V h/h km + j m C 2 d mj V d mj + i j C 3 d ij V s ij q ij (23) s.t. Equations (
k m Z ijkm + q ij = 1 ∀i, j = i (24) j m Z ijkm ≤ M Y ik ∀i, k (25)
i k

Z ijkm ≤ M Y jm ∀j, m (26)
l

L il Y ik - j L ij q ij ≤ Q 2 V c ik ∀k, i = k (27) l L lj Y mj - i L ij q ij ≤ Q 2 V d mj ∀j, m = j (28) Z ijkm ≥ 0, Y ij , q ij ∈ {0, 1} , V c ik , V h/h km , V d mj ∈ Z + ∀i, j, k, m (29)
Constraints (24) state that i sends the load to j either through hubs or by direct shipment. Constraints ( 6) and ( 7) in the original model were replaced by constraints ( 25) and ( 26) to allow for direct shipments.

In these two constraints, M is a large constant and can be set to n. Finally, constraints ( 27) and ( 28) calculate the number of vehicles required for collection and distribution, respectively.

We can add two additional constraints ( 30) and ( 31) to tighten the above formulation by limiting the definition of direct links to the links that connect two non-hub nodes. Without changing solutions, these constraints enhance the computational speed:

q ij ≤ 1 -(Y ij + Y ji )/2, ∀(i, j) (30)
q ij ≤ 2 -(Y ii + Y jj ), ∀(i, j) (31)
Similarly to Section 2, the large number of integer and binary variables makes the above formulation computationally demanding. Therefore, we need to extend the HLP-GLC formulation to allow solving larger problems optimally. On the other hand, a positive intercept in the linear cost function could resemble a sunk cost for assigning a vehicle to direct shipments, which allows for solutions of higher quality.

To incorporate direct shipments in the HLP-GLC formulation, we change the objective function to

min i k (b 2 Y ik + a 2 j m Z ijkm L ij )d ik + j m (b 2 Y jm + a 2 i k Z ijkm L ij )d jm k m (b 1 w km + a 1 i j Z ijkm L ij )d km + i j q ij (b 3 + a 3 L ij ) d ij (32)
Moreover, similarly to the extended version of MHLP, we need Equation ( 24) to relate q ij to Z ijkm .

Likewise, we can also include Equations ( 30) and (31) to tighten the formulation.

In order to find a proper linear approximation for the stepwise function on direct links (i.e. linear parameters a 3 and b 3 ), one can consider cost functions specific to each direct link as Equation ( 18), which result in zero measurement deviation on these links. Using Corollary 2 on the other hand, provides a unique linear cost function for all direct connections, which depends only on the vehicle's characteristic.

Theorem 1 can also be used if the distribution of L ij is known. However, assuming less-than-truck load for all O-D pairs (i.e. 0 < L ij ≤ Q 3 , ∀ i, j), we can define a 3 and b 3 in an effective way.

Assuming a 3 = C3 Q3 , in order to have a zero expected measurement deviation on the direct links we get:

E C 3 L ij Q 3 = E(a 3 L ij + b 3 ) ⇔ C 3 = C 3 Q 3 E(L ij ) + b 3 ⇔ b 3 = C 3 (1 - μ Q 3 ).
Therefore, for the direct links, one can use a unique linear cost function (a 3 L ij + b 3 )d ij with:

a 3 = C 3 Q 3 , b 3 = C 3 (1 - μ Q 3 ). ( 33
)

13

This section on numerical results is organized as follows. In Section 5.1, we compare the solution gaps and measurement deviations for HLP-GLC under different linear cost functions. We discuss about the choice of proper cost function for a modular HLP, considering factors such as solution gap and computational time.

In Section 5.2 we conduct simulations to test the accuracy and robustness of Proposition 1 regarding the predicted measurement deviation, and the factors impacting this deviation. In Section 5.3, we compare the ability of the stepwise and GLC functions to model the impact of shipment frequency in hub and spoke networks. Finally in Section 5.4, we study the advantages of allowing for direct shipments in hub distribution systems. To run these tests, we use a computer of 16 GB RAM and 3.60 GHz processor, with OS Windows 7 and CPLEX 12.6, and a computational time limit of three hours.

Measuring solution gaps and choosing the proper cost function

As defined in Section 3, the solution gap δ determines how close the optimal decisions of HLP-GLC are to those of MHLP. This indicator will help to choose an appropriate GLC function to replace with the stepwise function when the MHLP formulation is computationally intractable. Similarly to Section 3.3, we define the relative solution gap, δ as

δ = T C stp (Y * lin ) T C stp (Y * stp ) -1
To have reliable insights on the solution gap δ, we use all three main benchmark datasets in the HLP literature, namely the CAB dataset of O'Kelly (1987), the AP dataset of [START_REF] Ernst | Efficient algorithms for the uncapacitated single allocation p-hub median problem[END_REF] and the Turkish dataset of [START_REF] Tan | A hub covering model for cargo delivery systems[END_REF].

Regarding the stepwise cost function used in the MHLP formulation, we consider two transportation modes: large and small trucks with capacities of Q 1 = 4.5 and Q 2 = 1 tonnes, with costs of C 1 = 1.93 and C 2 = 1.45 dollars per kilometer respectively 1 . These parameters result in an α close to the one recommended by AP data set (0.296 against 0.25 2 ). Furthermore, similarly to [START_REF] Contreras | Hub location as the minimization of a supermodular set function[END_REF], we consider identical collection and distribution costs, as we assume that both are performed by the same type of vehicle. Finally for simplicity, we assume Q 3 = Q 2 and C 3 = C 2 , i.e. c/d links and direct links use the same type of vehicle.

We solve HLP-GLC using four different linear cost functions, and under two scenarios of allowing or forbidding direct shipments. In the absence of direct shipments, we use the MHLP and HLP-GLC

1 The data come from section 14 of the Owner Drivers and Forestry Contractors Act 2005 of Australia 2 Value of α according to AP dataset is equal to transfer cost collection cost = 0.75 3 = 0.25 19),(21) b: from Eq (19),( 21)

c/d links a = C2 Q2 a = C2 Q2 a = C2 Q2 a i,c = j Lij Q2 C2 j Lij , a j,d = i Lij Q2 C2 i Lij b = 0 b = C2 2 b: from Eq (19),(20) b = 0 h/h links a = C1 Q1 a = C1 Q1 a = C1 Q1 a = C1 Q1 b = 0 b = C1 2 b: from Eq (
direct links a = C2 Q2 a = C2 Q2 a = C2 Q2 a ij = Lij Q2 b = 0 b = C2 2 b : C 2 (1 -μ Q2 ) b = 0
Table 1: Calibration of the four GLC functions formulations presented in Section 2, while in the presence of direct shipments, the extended versions of the formulations presented in Section 4 are used. The four GLC functions are:

• GLC 0 with a slope of a t = Ct Qt for t = 1, 2, 3, and intercept b t = 0; we use this cost function to represent the HLP formulations with a zero intercept in the literature;

• GLC 1 with a t = Ct Qt and b t = Ct 2 recommended by Corollary 2;

• GLC 2 with a t = Ct Qt and intercepts based on the overall load of the network, i.e. using Theorem 1 for c/d and h/h links, and Equation (33) for direct links.

• GLC 3 where we determine specific cost functions for each c/d1 and direct links, but use Theorem 1 for h/h links.

Table 1 presents the parameters of the four aforementioned GLC functions. To summarize, GLC 0 represents a simple linear cost function with zero intercept, GLC 1 uses a generalized linear cost function which only uses vehicle's cost and capacity, and GLC 2 considers the load information as well, to determine the cost functions. These three cost functions identify a single linear cost function for each family of links.

Finally, GLC 3 provides a taylored cost function for each c/d and direct link, although it uses the same cost function for all h/h links.

We test instances of size 10, 15, 20 and 25, with 2 and 3 hubs. Since the AP data set does not have an instance of size 15, we use the first 15×15 matrices from the load and distance matrices of the problem of size n = 40 of this dataset. Similarly, for the Turkish data set we use the first n × n matrices of load and distance to generate datasets of different size n.

We introduce a new parameter to have a better understanding of the problem. For the HLP formulations with a linear cost function of zero intercept, multiplying the demand matrix by a constant does not change the location/allocation decisions; whereas it may change the solutions of the MHLP or HLP-GLC (with a cost function of non-zero intercept). To study the impact of the load's magnitude being transported over the network, we define β as the average number of vehicles required to carry loads on a c/d link.

β = (n -1)μ Q 2 . ( 34
)
This parameter can be used further to study the shipment frequency in the supply network (as we will see in Section 5.3), or to choose the proper transportation mode given the loads to be distributed over the network. We consider scenarios of β equal to 0.25, 0.5, 0.75, 1 and 1.5. For each case, the load matrix is multiplied by a factor in order to get the favorite value of β.

For each set of n, p and β, we solve the MHLP and HLP-GLC (under each of the aforementioned GLCs), for each of the three data sets. We then calculate the relative solution gap δ, as well as the absolute value of relative measurement deviation |˜ |, and compute the average of these two indicators over the three datasets,for each problem. Table 2 reports the average relative solution gap, and average relative measurement deviation (in parentheses) for each HLP-GLC formulation. Moreover, Table 3 provides the average of these two indicators over all test problems. For simplicity, henceforth we refer to different HLP-GLC formulations, only with the name of their cost functions; for example, in Table 2 GLC 0 refers to the HLP-GLC formulation with the cost function GLC 0 . Note that due to the high computational requirements of the MHLP formulation with direct shipments, we only consider samples of size 10 and 15 when direct shipment is allowed.

In Table 2, * implies that at least for one of the three datasets, we could not find the optimal solution of MHLP within three hours; for these cases, we consider the best integer solution found. Indeed for these cases, δ is a lower bound for its true value. Moreover † implies that the solution of HLP-GLC could not be found within three hours.

Without direct shipments, the highest relative solution gap is for GLC 0 , with an average δ of 6.03%, while it exceeds 20% solution gaps for some problems. However, as β increases, GLC 0 results in smaller δ. Nonetheless in GLC 0 , allowing for direct shipments, while charging the same transportation rate as for c/d links, results in solutions of poor quality, as our results show an average δ of 449.25%.

GLC 1 and GLC 2 result in smaller gaps: in the absence of direct shipments, the average δ is 1.34% for GLC 1 and 0.94% for GLC 2 . Allowing for direct shipments, increases the average solution gaps of these two cost functions to 2.74% and 2.13%, respectively. For both formulations δ is usually small, especially in the absence of direct shipments. Although statistically speaking, GLC 2 outperforms GLC 1 (using the paired t-test, the p-value=0.0452 when direct shipment is not allowed, and p-value=0.0493 when it is allowed), the difference between the two linear approximations is not large. Finally, GLC 3 outperforms all other HLP-GLC formulations, with an average δ of 0.34% in the absence of direct shipments and 1.80% in the presence of direct shipments.

p = 2 p = 3 β GLC0 GLC1 GLC2 GLC3 GLC0 GLC1 GLC2 GLC3
A partial correlation test with control variables of n, p, β and the type of GLC, indicates a significant positive correlation between δ and |˜ |, at 0.001 significance level. The results are similar for correlation between δ and | |. This finding implies that the magnitude of measurement deviation is a reliable predictor of the size of solution gap, regardless of the choice of linear function, problem size, etc. Furthermore, a similar analysis indicates a negative correlation between δ and β at 0.001 significance level, with n, p and the type of cost function as control variables. In other words, as β increases, the gap between the solutions of MHLP and HLP-GLC decreases.

Moreover, we found that there is no significant relation between the problem size n and the relative solution gap δ, when approximating stepwise cost by a GLC function. Conducting a partial correlation Direct shipment is NOT allowed Direct shipment is allowed test between δ and n, with the type of linear function as control variable, we have a p-value of 0.9882 when direct shipment is not allowed, and p-value of 0.8875 when it is allowed. Note that for some of the larger problems we could not find the optimal solution of MHLP within the time limit, implying that the solution gap could be even larger, which strengthens our findings.

GLC 0 GLC 1 GLC 2 GLC 3 GLC 0 GLC 1 GLC 2 GLC 3 average δ 6
Table 4 reports the average computational time over the three datasets. Note that * implies that at least for one of the datasets, the optimal solution could not be found within the three hours time limit. As we can see in Table 4, in the absence of direct shipments, GLC 0 provides the lowest computational time, while there is no significant difference between GLC 1 to GLC 3 , that use linear functions with non-zero intercept. Note that the negative impact of increasing p on the efficiency of GLC 1 to GLC 3 is stronger than that on GLC 0 . Not surprisingly, all formulations with a linear cost function outperform the MHLP in terms of computational time.

p = 2 p = 3 β MHLP GLC 0 GLC 1 GLC 2 GLC 3 MHLP GLC 0 GLC 1 GLC 2 GLC 3 Direct shipment is NOT allowed n = 10 0.
Allowing for direct shipments, however increases the computational time for all formulations, although its negative impact is outstanding for GLC 01 . The same holds true for the MHLP formulation, where in most cases, the optimal solutions could not be found within the 3 hours time limit for n ≥ 20.

Based on our discussions on Tables 2 to 4, we can come to a conclusion regarding the choice of the cost function. First and foremost, the solving approach is a key determinant. A compact formulation with less integer variables is usually preferred for exact solving, although the use of exact methods is limited to small networks (for MHLP, the compact formulation tested solved problems with size up to n = 25, whereas problems of size n = 40 could be solved by the taylored branch-and-bound algorithm of [START_REF] Tanash | An exact algorithm for the modular hub location problem with single assignments[END_REF]). As for methods that explore the feasible space with neighborhood search (e.g. metaheuristics) -which constitute a large part of the HLP literature [START_REF] Farahani | Hub location problems: A review of models, classification, solution techniques, and applications[END_REF])-the choice of the cost function does not usually make a difference in computational time; therefore, with these solution approaches, a stepwise function for the modular HLP may be preferred.

If one wishes to directly solve the problem to optimality by a standard solver such as Cplex, the MHLP formulation is not recommended even for medium size problems, given the large computational time required. In such cases, GLC 1 , GLC 2 and GLC 3 cost functions are preferred, as they usually offer a low solution gap with a reasonable computational time. The GLC 0 formulation could be used either if β is large enough (resulting in small solution gaps), or when n and p are too large to be addressed by other GLC functions. However, in general, GLC 0 is outperformed by other formulations.

If charging a unique transportation cost on every c/d (or direct) link is not necessary, GLC 3 is preferred over GLC 2 and GLC 1 as it provides the lowest solution gap with similar computational time. Otherwise, between GLC 1 and GLC 2 , which both offer a single cost function for each family of links, the latter provides smaller solution gaps in average, while the former is simpler to use, as it only depends on the vehicles type.

Figure 2 shows the tradeoff between different cost functions for a specific problem, which is somehow a representative of the above discussions. Similarly to Section 3, we assume a MHLP with no possibility of direct shipments, and two transportation modes, with the smaller vehicle for c/d links and larger vehicle for h/h links. We generate networks of n demand points, with coordinates x and y randomly distributed following a uniform distribution. Loads L ij between each O-D pair i and j also follow a uniform distribution from 0 to Q 2 , although in Figure 5, we try different load distributions in order to test the robustness of the predictions.

As we do not need the optimal solution to study the measurement deviation, we use a fast heuristic to generate solutions so as to make conclusions based on a large number of instances. We use the K-clustering heuristic of [START_REF] Lloyd | Least squares quantization in pcm[END_REF], where n nodes are allocated among p clusters, and the closest node to the weightcenter of each cluster is chosen as the hub. Given the location/allocation decisions, we compute the total transportation cost of the resulting solution using both linear and stepwise cost functions, and compute the relative measurement deviation ˜ for the problem. On the other hand, we use Equation ( 22) to predict the measurement deviation ˆ for the same problem. We estimate the value of ρ in Equation ( 22

a) at = C t Q t , bt = 0 -=1 -=4 -=7 -=10 (
Average measurement error (%) Figure 4 shows that a larger β (which could be associated with a lower shipment frequency or smaller vehicles for the c/d links) results in smaller measurement deviations. This result is intuitive; as an example, if the total load on a c/d link is L l = 0.5Q 2 , we observe a 100% measurement deviation on that link (since the stepwise function considers the cost of a full vehicle for this load, while the linear function only pays for half of the vehicle's capacity), however for L l = 10.5Q 2 , this deviation decreases to 5%. Although in both cases, one truck is dispatched half empty, the difference comes from the relative emptiness of the vehicles.

This explains why in Table 2, the relative measurement deviations of the GLC 0 are small for larger values of β. There are other factors that can impact the measurement deviation. For example, a higher economy of scale of the large vehicle comparing to the small one, leads to larger relative measurement deviations.

However for conciseness purposes, we skip analyzing all effective factors on the measurement deviation.

Impact of shipment frequency

Determining the optimal location of distribution centers generally requires a trade-off among inventory costs, transportation costs, customer responsiveness, etc. [START_REF] Nozick | Inventory, transportation, service quality and the location of distribution centers[END_REF]. [START_REF] Daskin | Facility location in supply chain design[END_REF] and Shen (2007) provide overviews of the research that study facility location in an integrated supply chain design. Similarly, when locating hubs in a supply chain network, we may consider other factors such as inventory policies and shipment frequency. In this section, we analyze the ability of MHLP and the various GLC formulations to capture the impact of shipment frequency on the optimal cost over a given time period.

To that purpose, we solve a number of problems from the CAB dataset, under different values of n and β, using three formulations: MHLP, HLP-GLC 0 and HLP-GLC 3 . We remind from Subsection 5.1 that parameter β reflects the shipment frequency; a larger β (i.e. higher average load per vehicle) can be associated with a lower shipment frequency throughout the supply chain, and therefore a smaller overall transportation cost is expected due to consolidation of loads. Our goal is to see how the overall transportation cost in a sufficiently long period of time is impacted by the shipment frequency, and how the different formulations for the modular HLP can take it into account.

To that purpose, we consider different values of β as different shipment frequencies, with β = 1.5 representing the lowest shipment frequency. Smaller βs refer to more shipments in the same unit of time.

For example, if β = 1.5 refer to the shipment frequency of once every 6 days, β = 0.25 would refer to a daily shipment pattern. To compare these scenarios, we need to compute the overall aggregated transportation cost in a longer period, e.g. 6 days, for all β. This is simply done by multiplying the transportation cost under each β i by 1.5/β i . For our computations, we use the assumptions presented in Subsection 5.1 regarding the stepwise and linear functions.

We can see in Figure 6 that both stepwise cost (Figure 6a) and GLC with non-zero intercept (Figure 6b) reflect the variation of the transportation cost when shipment frequency changes. The reason is the fixed costs incurred by these cost functions for every single transportation. However GLC 0 , having no fixed cost or intercept, is unable to model this consolidation effect (Figure 6c). This could help to enlighten future research on the HLP in an integrated supply chain design.

Impact of direct shipment

Next, we study the advantages of allowing direct shipments in a hub and spoke network under the modular framework. To measure the cost reduction due to the possibility of direct shipments, we compare the total transportation cost calculated by the MHLP formulation, in the presence and absence of direct average cost reductions.

Figure 7 shows that the advantage of using direct shipments increases in n and β. Moreover, it implies that by using specific vehicles for direct shipments, the resulting cost saving could be even higher. Obviously, by tailoring the transportation mode used for direct shipments, to the loads' distribution, one can achieve further cost reductions.

Figure 8 demonstrates share of direct shipments in the total transportation cost. As we can see, under larger values of β, direct shipments account for a considerable fraction of the total transportation costs, highlighting the importance of having direct shipments beside the through-hub transportations.

Conclusion

In this paper, we studied the modular hub location problem, where the transportation cost is a stepwise function of the number of vehicles. Since the MILP formulations with stepwise cost functions are computationally demanding, we approximated the stepwise cost with a linear function with non-zero intercept, based on a probabilistic analysis to minimize the expected measurement deviation between the two cost functions. We showed empirically that this measurement deviation is positively correlated with the gap between the optimal solutions of the two formulations. We also provided a prediction of the average measurement deviation between a stepwise cost function and any given linear approximation different from the best-approximation one. The accuracy and ronbustness of this prediction is confirmed by numerical experiments. Finally, our results show that allowing for direct shipments within a modular HLP can result in significant cost reduction, which could not be properly measured by using a simple linear cost function.

We finish by highlighting avenues for future research. In this paper we solved small size problems due to the computational limitation of solving the MILP formulations to optimality by a standard solver.

However, efficient solution algorithms are required to solve medium and large size MHLP or their HLP-GLC approximations. While we calibrated the GLC function to approximate the stepwise function for a given problem, a similar approach might be used to approximate nonlinear or piecewise linear cost functions.

Furthermore, using a stepwise cost function allows for better modeling of some extensions of the problem, such as allocation of multiple hubs to a terminal node, simultaneous intermodal transportation, limited fleet size and transportation capacity, HLP with routing in a dynamic setting, etc. Finally, the link with fleet management in a supply chain could be explored further. The type and number of required vehicles, or limitations in using them, could impact the location/allocation decisions. Having a holistic approach to handle these issues, though complicated, could result in significant cost savings for companies.

Figure 1 :

 1 Figure 1: Stepwise vs linear cost functions of the load, when C and Q represent the cost and capacity of the vehicle

  ), we calculate the total linear cost of Y , T C lin (Y ) = l∈Y Cost lin (l), where l ∈ Y refers to the subset of links used in solution Y . Similarly, we can calculate the stepwise cost of Y as T C stp (Y ) = l∈Y Cost stp (l). The two values are not necessarily the same, and we refer to the difference as the measurement deviation, denoted by : (Y ) = T C stp (Y ) -T C lin (* stp and Y * lin are optimal solutions of the MHLP and the HLP-GLC formulations, respectively. If the location/allocation decisions of Y * stp and Y * lin are not identical, they are likely to result in different stepwise costs, i.e. T C stp (Y * lin ) = T C stp (Y * lin

Figure 2 :

 2 Figure 2: Stepwise cost of optimal GLC solutions (a) and computational time (b) for n = 25, p = 2, β = 0.5, for CAB data set in the absence of direct shipments

Figure 3 :

 3 Figure 3: The effect of GLC slope and intercept on the relative measurement deviation ˜ (n = 100, p = 2, and β = 2.5)

Figure 3

 3 Figure 3 reports the measurement deviation under different GLC parameters. As expected, the linear cost function with a t = CtQt and b t = Ct 2 , with t ∈ {1, 2} proposed by Corollary 2, results in a very small relative measurement deviation. However, other linear cost functions with different slopes and intercepts may offer small measurement deviations as well, if their parameters are chosen properly.

Figure 4 :

 4 Figure 4: The effect of β on ˜ (n = 200 and p = 8)

Figure 5 :

 5 Figure 5: The effect of distributions of L ij on the measurement deviation , with n = 100, p = 3, and β = 1.5

Figure 6 :

 6 Figure6: Impact of β on the overall aggregated transportation cost, with p = 2 for the CAB dataset when direct shipment is not allowed.

Figure 7 :

 7 Figure 7: Average cost reduction under MHLP due to using direct shipments, with p = 2.

Figure 8 :

 8 Figure 8: Share of direct shipments in the overall transportation cost, with p = 2 (Results for MHLP).

Table 2 :

 2 Average relative solution gap δ (measurement deviation ˜ ) in %

Table 3 :

 3 Average relative solution gap δ and measurement deviation ˜ in % over all test problems

	.03	1.34	0.94	0.34	449.25	2.74	2.13	1.80
	average ˜ 189.86	16.15	6.97	2.13	1793.53	17.46	6.28	2.95

Table 4 :

 4 Average computational time in seconds

			25 1.7	0.3	0.5	0.5	0.5	16.1	0.2	0.5	0.7	0.7
			0.5	2.1	0.3	0.5	0.5	0.5	17.6	0.3	0.5	0.6	0.6
			0.75 1.8	0.2	0.4	0.4	0.4	50.3	0.3	0.5	0.6	0.6
			1	1.7	0.3	0.4	0.4	0.4	42.7	0.2	0.5	0.5	0.6
			1.5	0.7	0.3	0.4	0.3	0.4	6.2	0.2	0.5	0.5	0.5
			0.25 32.9	0.8	2.6	2.5	3.5	576.8	0.9	7.8	8.1	10.0
			0.5	23.6	0.8	2.2	2.1	2.7	9803.0 * 0.9	5.9	8.0	7.5
		n = 15	0.75 28.2	0.9	1.9	1.7	1.7	5471.0	0.9	4.6	6.7	7.4
			1	10.8	0.8	1.9	1.8	1.4	1724.4	0.9	4.4	4.7	6.4
			1.5	13.7	0.9	1.8	2.0	2.1	173.9	0.9	4.7	3.4	2.9
			0.25 507.2	4.2	12.2	10.8	15.5	10800 *	3.7	46.3	46.5	50.1
			0.5	225.1	4.2	11.8	10.4	10.1	10800 *	3.8	43.3	54.0	50.6
		n = 20	0.75 164.9	4.1	9.8	6.7	8.2	10800 *	4.0	30.1	50.2	37.4
			1	87.2	4.0	8.9	12.4	11.8	10513.5 * 3.8	23.4	19.5	31.4
			1.5	111.7	4.0	7.5	8.4	6.7	3805.0 * 3.7	19.1	18.1	15.5
			0.25 6462.6 * 17.0	51.1	49.8	81.2	10800 *	21.8	253.1	340.5	260.3
			0.5	2743.9	17.4	51.2	30.0	47.1	10800 *	22.6	146.3	210.1	185.0
		n = 25	0.75 2531.1	15.6	48.4	52.2	57.6	10800 *	20.4	135.4	137.4	146.0
			1	2558.3	16.0	39.3	45.8	45.6	10800 *	22.2	100.7	99.6	100.0
			1.5	1639.4	14.9	34.7	34.4	43.6	7763.1 * 19.4	102.7	87.1	77.2
			0.25 3.2	24.4	0.7	0.7	0.8	22.0	1.7	1.0	1.0	0.9
	Direct shipment is allowed	n = 10 n = 15	0.5 0.75 53.3 5.4 1 83.8 1.5 430.6 0.25 54.3 0.5 1276.7 0.75 4647.2 * 3636.1 * 4.8 24.7 0.6 25.5 0.6 24.6 0.6 25.2 0.6 3647.4 * 5.1 3644.3 * 5.2 1 7643.4 * 3627.2 * 5.2	0.7 0.6 0.6 0.6 5.5 5.4 5.2 4.9	0.7 0.6 0.8 0.7 6.3 5.7 5.3 4.6	28.1 224.3 221.3 264.0 755.9 10800 * 10800 * 10800 *	1.6 1.7 1.7 1.8 824.7 817.8 844.5 820.0	0.9 0.9 1.0 1.1 12.9 12.3 17.0 13.5	0.8 1.5 1.3 1.2 13.3 17.4 20.6 16.7	0.8 0.8 1.1 2.2 16.9 13.9 32.6 17.3
			1.5	10800.0 * 3629.7 * 5.4	5.6	5.8	10800 *	816.3	15.3	10.7	12.3

Note that in presence of direct shipments, in order to identify the a specific cost function for each link, we simply consider the whole flow that would pass over the link if the terminal node had no direct connections.

The long average computational time of GLC 0 with direct shipment, is mainly due to the large computational time for the Turkish dataset, however, even the CAB and AP datasets reported a larger computation time under GLC 0 than GLC 1 to GLC 3 .

Appendix A Proofs

Proof of Theorem 1: Since the load on a c/d or a h/h link includes the load of several O-D pairs, and can be approximated by a Normal distribution with mean µ l , which implies P(0 ≤ L l ≤ 2µ l ) 1, as L l > 0.

Moreover, since the Normal distribution is symmetric around its mean value µ l , we have for any u < v:

We have:

Proof of Corollary 1: If µ l is a multiple of Q 2 we have

Proof of Corollary 2: Using the approximation

the approximation holding when 2µ l /Q is large, which results in b = C 2 . This ends the proof.

Proof of Proposition 1: From Equations ( 1) and ( 11) we can compute the expected measurement deviation on a given link l as:

where t = 1 refers to h/h links and t = 2 refers to c/d links.

Let d denote the average length of a c/d link in the hub and spoke network, while ρd represents the average length of a h/h link. We first show that the expected measurement deviation under cluster symmetry is:

Indeed, the number of c/d links in the network is 2(n -p), while there exist p(p -1) h/h links, given the cluster symmetry assumption, as the network among the p hubs is a complete graph. The expected measurement deviation can be written as:

Moreover we have:

Notice that in the last equation, p-1 p refers to the average fraction of the load originated from one demand point that have destinations outside of the source cluster. Dividing Equation ( 36) by ( 37) provides an approximation of the relative measurement deviation, ˆ 1 .

1 Note that we also calculated E( T C lin ) directly, using numerical methods. However, Equation ( 22) is preferred, as it revealed to provide empirically the same level of precision of the prediction, with a closed-loop formula that is analytically explicit and much faster to compute.

Appendix B Estimation of ρ

In order to estimate ρ in Equation ( 22), we use the K-means clustering approach [START_REF] Lloyd | Least squares quantization in pcm[END_REF]. This heuristic divides n data points into K clusters so as to minimize the sum squared error between the empirical mean of a cluster and the points in that cluster. Similarly to K-means, the HLP partitions the n demand points into K = p clusters, with a hub that resembles the centroid in the K-means clustering. To estimate ρ, we generate n = 400 nodes spread over the x-y space according to a Uniform distribution, and construct p clusters using K-means clustering with K = p. We consider the centroid of each cluster as the location of the hub, and calculate ρ as the ratio of the average distance between the centroids of different clusters, divided by the average distance from each node to its cluster's centroid.

We solve this problem for different values of p ranging from 2 to 20, with 30 replications for each p.

Figure 9 shows the 95% confidence interval for the values of ρ. We can see that by increasing p, ρ increases as well. We tried different values for n as well, however the problem size does not have a significant effect on ρ.


Stefano Valtolina

Barbara Rita Barricelli

Alessandro Rizzi

Sabrina Menghini

Ascanio Ciriaci

Socio-Technical Design of an App for Migrants Rescue Operations

Keywords: Participatory design, Human Work Interaction Design, international development, domain experts, migration, emergency medical services

This position paper illustrates the research and development work done in the last years for understanding how to support domain experts in the rescue operations of migrants who attempt to reach Italian coasts via sea journeys on Mediterranean routes. The context, characterized by humanitarian, social, and organizational issues, presents complex challenges that can only be tackled with a multidisciplinary, participatory, and internationalized approach. The final system takes inspiration from other projects developed for different purposes but acting in similar context and according to similar socio-technical dimensions.

HWID for Emergency Medical Services

In recent years, Italy is handling the difficult situation of migratory flows ending with landing on the southern coast of the Country. According to UNHCR (United Nations High Commissioner for Refugees) report [START_REF]The Mediterranean Refugees/Migrants Data Portal UNHCR[END_REF], since the beginning of 2017, 50,275 immigrants have entered Italy from the Mediterranean. In [START_REF] Who | Mass casualty management systems strategies and guidelines for building health sector capacity[END_REF] World Health Organization defines a mass casualty incident as "an event which generates more patients at one time than locally available resources can manage using routine procedures. It requires exceptional emergency arrangements and additional or extraordinary assistance". This definition is well suited to describe what happens during the rescue operations for managing immigrants 'landing when it is necessary to offer medical assistance to a number of people who often exceed what the relief structures can accommodate. These conditions make it essential to study and implement specific strategies and work plans observed by all actors involved. Assistance operations are carried out by staff specialized in various disciplines such as doctors, nurses, and paramedics. Experts from different domains are called upon to perform their profession in a complex environment and with very challenging timing and intervention modes. This means having to perform sensitive tasks in a short time, though maintaining a high level of security, efficiency, and reliability of performance. This leads to the need of designing and developing IT applications to support the whole rescue operations. In particular, our research and development work is framed into studying how to design the interaction of IT solutions for enabling Emergency Medical Services (EMS). EMS are defined in [START_REF]Emergency medical services systems in the European Union[END_REF] as "[…] the ambulance services component that responds to the scene of a medical or surgical emergency, stabilizes the victim of a sudden illness or injury by providing emergency medical treatment at the scene and transports the patient to a medical facility for definitive treatment". However, the issues in this context do not relate only with medical assistance: operation workflows have to be put in place and leadership and organizational aspects have to be faced. To manage rescue operation in an efficient way means to nominate one or more managers who can successfully lead and coordinate all team members. One of the most critical actions to be taken in rescue operations is the triage, i.e. efficiently determining severities of injuries and prioritizing treatments; this action constitutes one of the most important tasks for basic life support. Furthermore, there are other crucial activities, like organizing and running specifics areas of operation for triage, treatment and transportation. Therefore, it is mandatory to take care of the specific flow of information between the operation managers and the team members. Such a research context can be clearly seen as framed into Human Work Interaction Design (HWID) [START_REF]Human Work Interaction Design: Designing for Human Work[END_REF][START_REF]Human Work Interaction Design: Usability in Social, Cultural and Organizational Contexts[END_REF][START_REF]Human Work Interaction Design. Work Analysis and HCI[END_REF][START_REF]Human Work Interaction Design[END_REF], a lightweight version of Cognitive Work Analysis, addressing the concept of Work in Human-Computer Interaction. The background and experience we bring in the field is twofold. On one hand, INMM -In Manibus Meis -is a registered supplier to NATO and is responsible for providing medical information support systems for first responders and military rescuers, medical control systems for first aid rescue teams, supporting systematic collaborative to emergencies management. On the other hand, the members of MIPS (Multimedia Interaction Perception Society) Laboratory of Università degli Studi di Milano bring into play their experience in interaction design for domain experts in several application domains [START_REF] Barricelli | Cultures of participation in the digital age: Coping with information, participation, and collaboration overload[END_REF][START_REF] Barricelli | Culture of participation in the digital age: Social computing for learning, working, and living[END_REF][START_REF] Barricelli | MANTRA: Mobile anticoagulant therapy management[END_REF][START_REF] Zhu | A meta-reflective wiki for collaborative design[END_REF][START_REF] Gianni | An ecosystem of tools and methods for archeological research[END_REF].

Motivations

The project aims at investigating how mobile digital tools used during critical and emergency activities not only facilitate the intervention on the field but also guarantee a benefit to whole related organizational structure.

The idea takes inspiration from other projects carried out in different contexts but designed around similar socio-technical principles. One example can be mKrishi@fisheries (see Fig. 1), an app developed by Tata Consultancy Services (TCS) in collaboration with the Innovation Lab Mumbai, Central Marine Fisheries Research Institute (CMFRI), Indian National Center for Ocean Information Services (INCOIS) and Indian Council of Agricultural Research (ICAR). This app aims at providing fishers in the Maharashtra state of India with a geo-referenced map indicating to where there is likely to be many fish. Through an algorithm able to analyze satellite images, weather conditions, data of the past catch, this app can help fishers to find locations which higher concentration of fish can be reached in safety condition.

Besides a real help to fishers, the app offers great benefits to the network of local fishery societies. For example, it was possible to register a reduction of consumption of fuel, whose distribution is under the responsibility of the fishery societies. The equal access to the information for each fisher has favored their cooperation and collaboration. Moreover, the growth in the use of mKrishi has seen and is still seeing the birth of a platform that includes services for equal fish marketing and suggestions for a more sustainable fishing procedure.

Examples like this leads us to believe that a single app used to support rescue and assistance operations in Italy can improve the whole organizational structure favoring better cooperation between all operative units and government agencies in order to enable more equal immigrants' treatment and distribution among all European countries.

Challenges

In designing and developing IT solutions to be used in such a critical context, we identified six main challenges:

1. Time and resources management: the applications need to support the rescuers in gathering medical data and in managing the operation as a whole as quickly as possible but keeping a high-level quality of the actions. 2. Clinical risk reduction: electronic guides are provided to avoid incomplete and incorrect medical data collection in stressful situations, which could impair the final outcome exposing patients to possibility of mistakes along the rescue chain. 3. Human rights preservation: the entire workflow needs to be addressed efficiently and in reasonable time but always paying attention to not overcome human rights and dignity in the process. Specifically, discriminations on any ground have to be avoided. 4. Privacy preservation: medical data have to be managed in compliance with law requirements. This means to collect, store, protect and use all gathered data in conformance with the requirements of legislation and regulations, both on a National and EU level [13]. 5. Internationalization: there are two different aspects of the context that require an internationalized approach. Firstly, the migratory wave is characterized by a multiplicity of different nationalities. To enable the collection of medical data and informed consent, and to efficiently and effectively deploy medical care, any IT application has to be designed and developed in more than one language. Secondly, it is desirable to trigger an information exchange process in order to facilitate the transmission and analysis of the data between European countries. 6. Security: information security controls had to be implemented to protect databases against compromises of their confidentiality, integrity and availability.

Our Contribution

INMM in collaboration with researchers of Università degli Studi di Milano have designed and developed ITHEALTH (International Traveller Health Surveillance System), a digital tool that through a tablet device provides rescuers with a set of functionalities for gathering patients' medical data and for managing and coordinating rescue operations. ITHEALTH allows rescuers to assign a bracelet to the patient tagged with a unique alphanumeric code (manually, by reading an NFC chip or scanning a QR Code), and then screening is assisted by the system in use at entry points (seaports and on board of ships) and along with transfer and relocation of migrants and refuges. Each authorised caregiver is assigned a portable device (tablet or smartphone); each migrant is assigned a medical tag. All personal data is securely encrypted in it, the medical tag is worn by the migrant as his/her right to access to health care. The resulting electronic health records (EHRs) are automatically and securely stored locally both in the ITHEALTH storage system and on the medical tag. Only authorised personnel who have access credentials to the system can see the data, thus physically the data is sent and visible only on authorized devices: laptop or desktop computer. When any connection is available data can be transferred to a server, installed on a laptop computer, normally placed at an operating centre and/or at the Hospital. There are additional features included:

In case of arrival from an epidemiological 'area of risk', the system, updated with previous preloaded information, matches data and instruct caregivers to deepen screening while alerting of a possible threat. In case a disease is confirmed, the caregiver is assisted to command the prompt evacuation by adopting the EVAC protocol, hospitalisation ad or isolation of the case, sending an alarm to the main Institution that there may be the risk of a public health threat to monitor.

The workflow implemented in ITHEALTH follows official protocols and standard procedures, so that screening is guided through unified protocols, throughout the whole chain of care.

The digitalization of such protocols allows to face the first two challenges mentioned in the previous Sectioni.e. time and resource management and clinical risk reduction. The quality of the workflow is guaranteed by the implementation of standard protocols, whereas the digital processes enable a quick data collection, management, and delivery. When a rescuer collects data, an informed consent is shown to the migrant to be signed. This page is translated in a set of languages and offers information about the reasons behind the data acquisition process. Fig. 2 depicts a set of screenshots that present some steps of the medical acquisition process of ITHEALTH. Fig. 2(a) requires to insert or a new immigrant or to load data previously registered in an electronic bracelet assigned to an immigrant. The forms in Fig. 2 (c) and (d) are used for collecting immigrant's status and biographical data and medical data respectively. In this latter case, data concern heart rate, oxygen saturation and body temperature. Finally, Fig. 2 (d) requires immigrant signs the informed consent about the gathered data. This page is translated into the immigrant's language if it is possible.

Fig. 2. Screenshots of some pages of the ITHEALTH medical data acquisition process

When an intervention is completed, the resulting EHRs are automatically and securely stored locally both in the ITHEALTH storage system and on the medical TAG. When communications are available, the rescuer can send the EHRs and the additional information through any communication channel available (recently radio communications have been exploited too). The server keeps a database that securely stores all that has been done during the operations, enabling the creation of reporting to use for coordinating different rescues and for enabling the cooperation between various intervention agencies. About human rights and privacy presentation, EHR systems need to manage new and additional safeguards to address the fundamental conflicts and dangers of exchanging information in an electronic environment [13]. As said in the previous section, any IT application has to be designed and developed in more than one language. This is to mitigate misunderstandings between the rescuer and the immigrant about data processing and collection purposes could block or slow down the information transmission process. To face off this possible communication gaps ITHEALTH tries to report medial information by using appropriate images and very simple interfaces that can be understood even if the migrant does not speak one of the languages known by the rescuers. Fig. 3 presents another screenshot used for reporting possible pain, fractures, burns, wounds, skin manifestations or amputation. By sketching on the figure, the rescuers can annotate where the problem is located by using a color code to specify the type of problem. This visual strategy is very useful and simple to use for the rescuer and very understandable for the immigrant too. Finally, to deal with the second aspect concerning international issues we designed an information exchange protocol among the parties in order to: (i) set forth the information to be exchanged, the operational procedures to be followed, and the security mechanisms and other safeguards to be maintained; (ii) and set out the ways that such exchange of the particular information would be consistent with the purposes. To this aim, ITHEALTH provides modules that allow coordinators of the involved teams to follow remote rescue operations, giving orders, guiding the actions of the individual rescuer, recording data about injured, and setting up coordination tasks. In order to evaluate ITHEALTH we carried out several tests in different scenarios for testing how the tool can support rescuers during their actions and other rescue operations in order to decrease their workload while accomplishing several unusual tasks in parallel and under time pressure.

Fig. 1 .

 1 Fig. 1. A screenshot of the mKrishti app with indications to where finding fish (left) and Fisheries landing center in Alibaug village (right). Screenshots courtesy of TCS, Innovation Lab Mumbai, ICAR, CMFRI Mumbai, and INCOIS.

Fig. 3 .

 3 Fig. 3. A screenshot that depicts how medical data are gathered in a visual and intuitive way.