Diego Delle Donne
email: delledonne@essec.edu

Laurent Alfandari
email: alfandari@essec.edu

Claudia Archetti
email: archetti@essec.edu

Ivana Ljubić
email: ljubic@essec.edu

Freight-on-Transit for urban last-mile deliveries: A Strategic Planning Approach

Keywords: freight-on-transit, network flow, column generation, mixed integer linear programming

We study a delivery strategy for last-mile deliveries in urban areas which combines freight transportation with mass mobility systems with the goal of creating synergies contrasting negative externalities caused by transportation. The idea is to use the residual capacity on public transport means for moving freights within the city. In particular, the system is such that parcels are first transported from origins (central distribution centers) to drop-in stations, which are stop locations on public vehicles itineraries. Then, they are transported through public vehicles to drop-out stations, from where they are delivered to destination by freighters using green vehicles (such as bikes, drones, porters, etc.). The system is known as Freight-On-Transit (FOT). In this paper, we focus on the strategic decisions related to defining the public transportation network that will take part to the delivery system, i.e., which public vehicle lines and stop locations will be included (and thus equipped for the service). We propose different formulations for the problem and effective heuristic solution approaches based on column generation. We perform exhaustive tests aimed at providing managerial insights on the performance and the efficiency of the system.

1 Introduction E-commerce is experiencing an explosion of revenues and market share. According to [START_REF] Monday | E-commerce share of total global retail sales from 2015 to 2024[END_REF], e-commerce market share of retail sales is expected to reach almost 22% in 2024. This phenomenon has unprecedented impacts on people's daily life. Nowadays, and especially as a consequence of the COVID-19 pandemic, consumers are more and more used to buy on the internet instead of going to physical stores. In addition to the daily life alteration, another consequence of e-commerce explosion is related to the new logistic and distribution model: goods are now directly delivered to customers places. The effect is the huge increase of commercial vehicles travelling around urban and rural areas for home deliveries. Given the granularity of customers orders (e-buyers typically place small orders with high frequency, see [START_REF] Ipc | Online shopping frequency according to online shoppers worldwide as of October 2018[END_REF]) and the importance given to fast shipping (see [START_REF] Mhl | Delivery Time Top Priority for Online Shoppers[END_REF]), existing delivery operations related to e-commerce are becoming unsustainable. This is especially true when considering last-mile deliveries in urban areas, where most of online buyers concentrate. In China for example, more than 40 billions of packages were delivered in urban areas in 2017 (+ 28% in one year, see [START_REF] Of | Statistics of China[END_REF]). The main drawbacks associated with the steep increase of the number of (half-empty) vans travelling around urban roads are related to pollution, congestion and risk of accidents. According to [START_REF]Handbook on External Costs of Transport (European Commission)[END_REF], negative externalities generated by road freight transportation (mainly related to road accidents, congestion, pollution and climate change) represent a huge loss for society, equivalent to around 6.6% of EU GDP. Most of the economic and social losses are concentrated in the city and are around 2% to 4% of cities' GDP.

Various strategies have been recently developed to face the issue. Autonomous vehicles are seen as the future of last-mile deliveries (see [START_REF]Projected market size of the autonomous last mile delivery worldwide from 2021 to 2030 (in billion U.S. dollars)[END_REF]) as they are environmentally friendly and some of them, like drones and robots, avoid congestion. However, there are issues related to their use. For example, drones raise security issues given the density of population of urban areas. Amazon Prime air program, which uses drones for performing Amazon's deliveries, has been launched in 2016 but is not operative yet. Delivery robots with wheels (e.g., Yape, Starship, Amazon Scout, Fedex's SameDay Bot), using sideways at pedestrian speed, now effectively operate in a number of cities around the world (see [START_REF] Alfandari | A tailored Benders decomposition approach for last-mile delivery with autonomous robots[END_REF][START_REF] Mourad | Integrating autonomous delivery service into a passenger transportation system[END_REF] for related delivery optimization problems). Other non-autonomous solutions involve different green transportation means, like bicycles. Cargo-bikes are nowadays an established reality in many cities.

Combination of public transport systems with freight transportation

One opportunity that has not yet been fully analysed and exploited is combining mass mobility systems with freight transportation, in particular, using Public Transport Systems (PTS) for performing last-mile deliveries. Indeed, PTS suffer from a chronic inefficiency: they are overloaded on peak-hours and largely underused out of the peak (see [START_REF] Chapuis | Restoring public transit amid COVID-19: What European cities can learn from one another[END_REF]). There are thus high potential benefits in exploiting the residual capacity of PTS for last-mile deliveries, to face delivery sustainability and cope with PTS inefficiency. This novel logistics concept of integrating goods and passenger flows to promote higher utilization rates for the public transport network is known as Freight on Transit (FOT) (see, e.g., [START_REF] Ozturk | An optimization model for freight transport using urban rail transit[END_REF] and [START_REF] Cochrane | Moving freight on public transit: Best practices, challenges, and opportunities[END_REF]). Passenger-and-Package Sharing systems have been recommended by the European Commission as soon as 2007 ( [START_REF]Green paper, towards a new culture for urban mobility[END_REF]).

The interest in FOT is not only scientific, as shown in the literature review section, but also practical. Indeed, involved stakeholders, mainly public authorities, have started the discussion a few years ago (see for example, [START_REF] Edrington | Guidebook: Using Public Transportation to Facilitate Last Mile Package Delivery[END_REF]) and some applications already exist. In a promising FOT pilot project launched by Monoprix in Paris (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017), the company used the RER line D in Paris to transport goods between their distribution center in Combs-la-Ville and Bercy station in the city boundary of Paris [START_REF] Boudet | Monoprix choisit Fret SNCF pour approvisionner ses magasins parisiens[END_REF]. In addition, Amazon has recently won a patent for a mobile package pickup system using public transport [START_REF]Amazon packages to ride public transport[END_REF]. The system calls for installing storage compartments on buses or other vehicles, and letting riders with the authorised codes unlock a designated compartment and pick up their item. Other recent interesting initiatives include the municipality of Vienna, where the delivery service is in place combining FOT with cargo bikes. Packages are delivered to the bus garage where they are sorted on bike containers and delivered through cargo bikes (see [START_REF]Climate-Friendly Package Delivery with Cargo Bikes[END_REF]). A similar initiative is taking place in Barcelona [START_REF]El reparto de última milla se sube a la bici en el área metropolitana[END_REF].

In this paper, we investigate a system which is close to the one proposed by Amazon. Specifically, we consider a last-mile delivery service where parcels move from origin points, called Central Distribution Centers (CDC), to destinations (customers places) through a delivery service making use of PTS. In particular, parcels are first transported from CDC to PTS stations, called drop-in stations, where they are loaded on PTS buses (or any other PTS transportation mean, like tram). They are then moved by buses to drop-out stations, where they are unloaded and picked-up by freighters (which can be bikers, porters, drones, droids...) that perform the last leg of the delivery, i.e., up to the customer's place. Figure 1 illustrates the described transport system. In this work, we focus on the strategic decision related to the design of the PTS network used in the delivery service, i.e., the choice of the bus lines and stations to be used in the delivery service.

Strategic decisions

FOT systems involve two kinds of stakeholders: public authorities managing the PTS and logistic companies performing last-mile delivery operations. FOT systems rely on agreements so that logistic companies make use of PTS for performing their deliveries. These agreements might result from a win-win perspective (public authorities are interested in reducing the societal and environmental impact of delivery operations while logistic companies are interested in reducing their costs) or from regulations imposed by public authorities (like limitations on commercial vehicles that are allowed to travel in urban areas).

As mentioned above, in this paper we focus on the strategic decision associated with the design of a FOT system. In particular, we focus on the public authority decision about the setting of the system, i.e.:

1. which part of the PTS network should be involved in the FOT system, i.e., the bus/metro/tram lines that will take part in the service and that, thus, have to be equipped accordingly; Figure 1: Illustrative example of a FOT system.

2. which stops have to be equipped for acting as drop-in or drop-out stations.

Both decisions involve a set-up cost in terms of equipment for accommodating freights that will be transported. There will be storage place for drop-in and drop-out stations and storage compartments in vehicles.

Note that examples of such equipment already exist in practice, like lockers as storage space for stations or the storage equipment on public vehicles proposed by Amazon [START_REF]Amazon packages to ride public transport[END_REF]. Many different factors have to be taken into account in the strategic decision. For the system to be as effective as possible in terms of resource utilization and congestion reduction, the main objective is to maximize the demand served by existing PTS resources. However, budget limitations related to the PTS equipment mentioned above have to be taken into account, since public authorities may not allow to use all bus lines or all stations. Also, some other important factors have to be considered, which are related to the inconvenience caused by such a service on passengers mobility. Indeed, storage equipment, both on vehicles and stations, would reduce the capacity formerly allocated to mobility. In addition, operations related to loading and unloading freights might interfere with passengers activities. Thus, the public authority could introduce restrictions related, for example, to a subset of existing public stops or public vehicles which cannot be involved in the service (for example, because they are typically too crowded), or time windows in which the FOT system cannot operate (peak hours).

In the following, we assume that the latter restrictions are known, meaning that the public authority has already deliberated on this kind of limitations. Thus, the problem is to decide the subset of lines and stops to be involved in FOT in such a way that the flow of goods transported is maximized. We assume that an estimated origin-destination (O-D) matrix is available, where for each pair of origin and destination the corresponding amount of freight to be delivered is known. Thus, the problem becomes designing the public transportation leg of , i.e., deciding which bus lines and stops should be included, with the aim of maximizing the flow of freight that can be moved through FOT while satisfying the budget constraints. We call the problem Strategic Freight-on-Public-Transport problem (SFOT).

Our contribution

Our contribution can be summarized as follows.

• We propose a new combinatorial optimization problem to help public authorities in analyzing, developing and implementing this novel logistic concept, in which freight and passenger transportation are combined. The network design problem maximizes the covered demand, while respecting capacity restrictions reflecting the available budget and the PTS infrastructure. Major strategic decisions to be made are: which PTS lines and which stations should be used, and with which capacities.

• We propose three mixed integer programming (MIP) formulations for finding exact solutions of the problem. They are based on paths, flows and aggregated flows, respectively. We use instances of small and medium size to compare the computational performance of these models.

• We design an effective MIP heuristic based on a column generation approach associated with the path formulation, which is computationally the best-performing model among the three ones proposed.

These MIP heuristic provides high-quality solutions for large instances of realistic size.

• We provide useful managerial insights. In particular, we analyze the impact of the available budget and the available stop-and line-capacity on the demand that can be covered. In addition, we simulate scenarios in which the demand is increased, or the number of available lines is enlarged, to understand their impacts and limitations for achieving the best resource exploitation and absorption of the existing demand.

The paper is organized as follows. In Section 2 we give an overview of the related literature. In Section 3.2 we provide a formal problem definition and we introduce the three MIP formulations. A MIP heuristic derived from the column generation procedure is presented in Section 4. Computational results are provided in Section 5 and final conclusions are drawn in Section 6. Additional computational details are provided in the Appendix.

Literature review

We first start this review of literature with general papers about integrating freight in PTS, outside the field of Operations Research and without any optimization model. In [START_REF] Cochrane | Moving freight on public transit: Best practices, challenges, and opportunities[END_REF] first defines the term Freight on Transit (FOT) was introduced first. The authors conduct a survey from transportation experts to determine the greatest opportunities, challenges and difficulties of FOT. They depict very specific local solutions for the city of Toronto. One of the interesting ideas in this paper is to attach cargo trailers to buses. In [START_REF] Trentini | Toward a shared urban transport system ensuring passengers & goods cohabitation[END_REF] explore innovative solutions to integrate a global urban mobility strategy are explored, presenting new concepts for the TPS. In particular, they introduce two FOT projects, one in Germany and one in London, in which public buses carry small packages are presented. The paper shows some quantitative gains with the new systems. Different stages of analysis necessary to include freight into PTS are studied in [START_REF] Galkina | Investigating using Urban Public Transport For Freight Deliveries[END_REF]. The authors provide a case study based on the city of Bratislava, and propose a basic system. In [START_REF] Gatta | Sustainable urban freight transport adopting public transport-based crowdshipping for B2C deliveries[END_REF], the authors depict crowdshipping systems based on public transport services, with individuals carrying parcels during their daily commute using the metro of Rome. In [START_REF] Bruzzone | The integration of passenger and freight transport for first-last mile operations[END_REF], the integration of passenger and freight transport is studied, and various key performance indicators (KPIs) of such integration are proposed. A case study using vessels of the local PTS in the Lagoon of Venice is presented. In [START_REF] Zhou | Freight Transport Mode Based on Public Transport: Taking Parcel Delivery by Subway as an Example[END_REF], the authors study parcel delivery by subway. Their results show that using the subway to transport goods can efficiently reduce carbon emissions by 20-50 %. The interested reader will find a comprehensive list of papers in the recent survey of [START_REF] Elbert | Freight on urban public transportation: A systematic literature review[END_REF], reviewing both qualitative and quantitative research on the integration of freight in public transport. Note that, while in [START_REF] Trentini | Toward a shared urban transport system ensuring passengers & goods cohabitation[END_REF] no specific infrastructure is foreseen for the delivery service, the remaining contributions listed above all considers costs related to the equipment of the means of transportation and locations involved in the system. We now review contributions which employ the methodology of Operations Research and Optimization. To the best of our knowledge, very few among the existing articles study the strategic network design questions, i.e., how to select a subset of public transport buses and stations for freight transportation. We found three papers that propose MIP formulations: two of them model hub location problems ( [START_REF] Ji | A Multimodal Passenger-and-Package Sharing Network for Urban Logistics[END_REF] and [START_REF] Zhao | Location selection of intra-city distribution hubs in the metro-integrated logistics system[END_REF]) and one models a location-routing problem ( [START_REF] Gianessi | The multicommodity-ring location routing problem[END_REF]). In [START_REF] Ji | A Multimodal Passenger-and-Package Sharing Network for Urban Logistics[END_REF], packages are transported by metro between the hubs, and by taxis or trucks between stores and hubs. In [START_REF] Zhao | Location selection of intra-city distribution hubs in the metro-integrated logistics system[END_REF], same-day delivery during off-peak hours is performed using subway, and the model aims again at locating the metro hubs with a p-median facility location model, with a case study in Shanghai. In the location-routing problem of [START_REF] Gianessi | The multicommodity-ring location routing problem[END_REF], a subset of UDCs is to be selected around the city to be connected by a ring using public trains or shuttles, and electric vans perform pick-up and delivery inside the city. None of the above papers considers the selection of lines and maximum coverage of demand due to line capacities.

In a recent comprehensive survey dealing with freight on urban public transportation [START_REF] Elbert | Freight on urban public transportation: A systematic literature review[END_REF], 15 articles have been classified as quantitative. Among them, 11 deal with operational decisions, most of them involving vehicle routing. In the following, we provide a brief overview of the works that are most closely related to our contribution. In [START_REF] Cheng | When packages ride a bus: Towards efficient city-wide package distribution[END_REF], packages are distributed by buses over a pre-designed network of bus lines and stops. Decisions concern the routing of packages pre-assigned to a bus line. Transshipment is allowed and a multi-commodity flow formulation minimizing the time-span of deliveries, and a heuristic, are provided. In [START_REF] Cheng | Planning city-wide package distribution schemes using crowdsourced public transportation systems[END_REF], the authors consider a variant of the former problem in which a single path is given and transshipment is not allowed. In [START_REF] Masson | Optimization of a city logistics transportation system with mixed passengers and goods[END_REF] a single Consolidation and Distribution Center is given together with a single bus line with multiple buses, where bus stops are fixed and time windows are associated with customers. The authors provide a model based on the classical arc-based 2-echelon VRP, minimizing the number of vehicles used, breaking ties concerning the distance traveled by the freighters. The authors also propose an Adaptive Large Neighborhood Search (ANLS) metaheuristic. In addition, they derive lower bounds by solving the linear relaxation of a set-partitioning formulation by column generation. In [START_REF] Behiri | Urban freight transport using passenger rail network: Scientific issues and quantitative analysis[END_REF], the authors consider transporting freight using passenger trains to Paris, with a single rail line with fixed schedules and stations' capacities. The decisions are related to the assignment of packages to trains. A MIP formulation is proposed, however only heuristic methods are provided for computing solutions. The example of the French store Monoprix is given, transporting textiles, cosmetics, household goods and leisure products from warehouses in the outskirts of Paris to the Paris-Bercy station via RER D train line. In [START_REF] Ozturk | An optimization model for freight transport using urban rail transit[END_REF], an optimization model is proposed for city logistics sharing passenger trains and cargo trains. In [START_REF] Ghilas | Integrating passenger and freight transportation: Model formulation and insights[END_REF], a MIP formulation is provided for a problem that fits our configuration with depot, drop-in and drop-out stations and destinations of demand requests. However, only operational decisions concerning the routing of vehicles carrying both parcels and passengers are optimized. Finally, the approach of [START_REF] Mourad | Integrating autonomous delivery service into a passenger transportation system[END_REF] integrates autonomous pick-up and delivery robots sharing PTS capacity.

Note that the review of [START_REF] Elbert | Freight on urban public transportation: A systematic literature review[END_REF] classifies FOT integration systems in three categories, namely: Shared track where the public infrastructure like rails is used and the vehicles are generally light railways with only freight in the wagons but no passengers (such systems have been implemented in Dresden, Frankfurt and Zurich, see [START_REF] Marinov | Urban freight movement by rail[END_REF]); Shared vehicle where freight is transported inside a separate wagon attached to a light railway with passengers, or a trailer attached to a public bus, with dependent loading/unloading and transshipment operations (see [START_REF] Shen | A new urban logistics transport system based on a public transit service[END_REF] and [START_REF] Behiri | Urban freight transport using passenger rail network: Scientific issues and quantitative analysis[END_REF]); and finally, Shared Wagon, where passengers and freight travel inside the same wagon or bus and peak periods of passenger flow have to be taken into account to share space in a proper and secure way (see [START_REF] Kelly | Innovative interior designs for urban freight distribution using light rail systems[END_REF]).

To summarize, to the best of our knowledge, this article is the first attempt to optimize strategic decisions in the context of designing a FOT system with selection of lines. Moreover, most of the above mentioned articles minimize (operational) cost for satisfying a given demand, without considering the maximization of the covered demand.

Problem statement and formulations

In this section we first provide a detailed problem statement, followed by three MIP formulations.

• S in k = S in ∩ S in k • S out k = S out ∩ S out k • L in s /L out s ⊆ L:
Set of lines having stop s as a drop-in/drop-out stop

• L s := L in s ∪ L out s :
Set of lines passing through stop s Table 1: Notation summary of (public transport) lines L, each of which is associated with a fleet serving the line (i.e., each fleet serves a subset of lines). We let F be the set of fleets. A line ∈ L is associated with a sequence of potential drop-in stops S in followed by a sequence of potential drop-out stops S out . We define S = (S in ∪ S out ) and S = ∈L S as shortcuts. Note that we do not consider the case of 'circular lines', i.e., lines that make a go-and-back' path. This lines would not satisfy the assumption that drop-in stops precede drop-out stops. However, these lines can be easily modelled as two lines going in opposite directions. Each line ∈ L has a daily capacity C ∈ N of packages which can be transported. In turn, each stop s ∈ S has a daily capacity of C s ∈ N packages which can be stored at s (either as drops-in or drops-out). A set of commodities K is given, and each commodity k ∈ K has a given origin o k ∈ O, a given destination d k ∈ D and a number P k ∈ N of packages to be delivered from o k to d k . In addition, only a subset S in k of the drop-in stops (those nearby o k ) and a subset S out k of the drop-out stops (those nearby d k ) can be used to deliver commodity k. We define S in k = S in ∩ S in k and S out k = S out ∩ S out k as shortcuts. Finally, strategic budgets maxfleets, maxlines and maxstops (being non-negative integers) are given to limit the number of fleets, lines and stops, respectively, selected in the solution. We summarize in Table 1 the notation used in this paper.

In this work, we do not consider transshipment to be an option for delivering packages, this is, a package cannot be transported by more than one line. Therefore, the delivery of packages of a given commodity k ∈ K with the FOT system consists on three legs:

1. First leg: From origin o k to a drop-in stop s ∈ S in k .
2. Public-transportation leg: From stop s ∈ S in k , through a line ∈ L in s to a drop-out stop t ∈ S out k .

3. Last leg: From stop t to its final destination d k .

Note that the public-transportation leg is mandatory in the FOT system, i.e., a package cannot be taken from o k to a stop s and then from directly s to d k , even if s ∈ S in k ∩ S out k . We consider that in this case, the package could have been delivered directly, without using the FOT system.

As we are dealing with FOT at a strategical level, we consider each commodity to represent the (estimated) aggregated demand from one depot (origin) to a given neighborhood (destination). Hence, the set of packages corresponding to a commodity can be split and sent by different paths to their final destination.

Definition 1 (SFOT).

With the above given notation, the SFOT asks to select a subset of fleets, lines and stops, satisfying budget constraints and to route packages through these elements, maximizing the number of delivered packages.

MIP formulations

In the remainder of this section we propose three MIP formulations for the problem introduced in Section 3.1. The first is a path formulation where variables represent commodities' paths from origins to destinations. The remaining are two flow formulations. All presented formulations use the same set of binary variables to model discrete decisions concerning network design issues, i.e., the selection of fleets, lines and drop-in/drop-out stops to use:

• For each stop s ∈ S, a binary variable x s is used to decide if s is selected as part of the solution or not.

• We use binary variable y to decide if line ∈ L is selected in the solution.

• For each fleet f ∈ F , we use binary variable z f to state if at least one line ∈ L f is selected.

Each of the three formulations uses a different set of additional variables to model the flow of packages in the resulting network. We first present these formulations, and at the end of this section we introduce additional inequalities which can be used to strengthen their linear programming (LP) relaxations.

Path-based formulation

In addition to variables (x, y, z) introduced above, the path-based formulation uses a continuous variable w k st ∈ R + , for each demand k ∈ K, each line ∈ L and each pair of stops (s, t) ∈ S in k × S out k , to determine the flow of packages of commodity k using line , loaded at drop-in stop s and unloaded at drop-out stop t. Note that the choice of line along with the stops s and t determines a complete path from o k to d k and so each of the variables w k st is associated with a possible path from an origin to a destination. The path-based formulation for the SFOT can then be defined as follows:

max k∈K ∈L s∈S in k t∈S out k w k st (1) ∈L s∈S in k t∈S out k w k st ≤ P k ∀k ∈ K (2) k∈K s∈S in k t∈S out k w k st ≤ C y ∀ ∈ L (3) k∈K ∈L   t∈S out k w k st + t∈S in k w k ts   ≤ C s x s ∀s ∈ S (4) s∈S x s ≤ maxstops (5) ∈L y ≤ maxlines (6) f ∈F z f ≤ maxfleets (7)
s 2 s 4 o k1 d k1 s 3 o k2 d k2 s 1 s 5 s 6 s 7 Figure 2:
The graph G for a line with 3 drop-in stops and 3 drop-out stops by which two commodities, namely k 1 and k 2 , can be transported.

y ≤ z f ∀f ∈ F, ∈ L f (8) x s , y , z f ∈ {0, 1} ∀s ∈ S, ∈ L, f ∈ F (9)
w k st ∈ Z + ∀k ∈ K, ∈ L, s ∈ S in k , t ∈ S out k ( 10
)
Constraints ( 2) limit the aggregated flow of each demand along all possible paths. Constraints ( 3)-( 4) ensure that lines and stops capacities are not violated. Constraints ( 5), ( 6) and ( 7) impose strategical restrictions on the number of stops, lines and fleets, while constraints (8) ensure the logical link between variables z and y, i.e., a line can be selected only if the corresponding fleet is selected too. We point out that constraints (3) are valid under the assumption that there are no drop-in stops after a drop-out stop. Nevertheless, if this assumption does not hold, our formulation can be easily adapted to deal with this case.

A potential disadvantage of this model is its size: even though the formulation is compact, the number of variables is in O(|K||L||S| 2 ) and the number of constraints is in O(|K| + |L| + |S|). To make the model computationally more tractable, we propose in Section 4 to solve the LP-relaxation of this model by dynamically pricing its columns. In the remainder of this section, we propose two alternative formulations, with significantly less decision variables.

Flow-based formulation

For the flow-based formulation, we start by focusing on the flow of packages for a given line ∈ L. This can be done by considering each line separately, due to the fact that transshipment is not allowed. We construct a network, say G = (V , A ∪ A in ∪ A out ), representing a directed simple graph, in which V is the set of stops visited by line along with all origins o k ∈ O and destinations d k ∈ D for each commodity k which can be routed through , i.e.,

V = S ∪ {o k ∈ O : k ∈ K ∧ S in k = ∅} ∪ {d k ∈ D : k ∈ K ∧ S out k = ∅}.
We define three types of arcs:

• Set A are the arcs (s i , s i+1 ), i = 1, . . . , |S | -1,
connecting two consecutive stops of the line, where s i is the i-th stop on the path of line . We slightly abuse the notation here and omit index for stations s i , keeping in mind that we are referring to a single line .

• Set A in are the arcs (o k , s), for each k ∈ K and s ∈ S in k , and • Set A out are the arcs (t, d k ), for each k ∈ K and t ∈ S out k . Figure 2 shows the graph G for a line with 3 drop-in stops and 3 drop-out stops by which two commodities, namely k 1 and k 2 , can be transported.

For each commodity k ∈ K and each arc However, due to the structure of our network (line is a simple path), these variables can be safely projected out, as we shall see below in Proposition 1. Figure 3 illustrates the proposed network and the role of each flow variable in the model. We note that the set of arcs on this network is the union of the arcs of the graphs G , for every line ∈ L (characterized with different colors in the figure). Moreover, due to the fact that transhipment is not allowed, the flow on each graph G is independent from the flow on the other graphs. The set of variables of this formulation is completed by using the binary variables x, y and z introduced above. Thus, the flow-based formulation for the SFOT can be written as:

(o k , s) ∈ A in (resp. each arc (t, d k ) ∈ A out ), variable f in k s ∈ Z + (resp. f out k t ∈ Z + )
2 s 1 s 2 o k t 1 t 2 t 3 d k f in k 1s1 f in k 1s2 f in k 2s2 f out k 1t2 f out k 2t1 f out k 2s2 f out k 1 s2
max k∈K ∈L t∈S out k f out k t (11) s∈S in k ∈L in s f in k s ≤ P k ∀k ∈ K (12) k∈K s∈S in k f in k s ≤ C y ∀ ∈ L ( 13
)
∈L in s k∈K s∈S in k f in k s + ∈L out s k∈K s∈S out k f out k s ≤ C s x s ∀s ∈ S ( 14
)
s∈S in k f in k s = t∈S out k f out k t ∀k ∈ K, ∈ L (15)
f in k s ∈ Z + ∀k ∈ K, ∈ L, s ∈ S in k ( 16
)
f out k t ∈ Z + ∀k ∈ K, ∈ L, t ∈ S out k (17)
x, y, z satisfy ( 5)-( 9)

Constraints ( 12)-( 14) link flow-based and binary variables and limit the flow according to demands and capacities (in a similar way as constraints (2)-( 4) for the path-based formulation). In addition, flow conservation constraints [START_REF]Green paper, towards a new culture for urban mobility[END_REF] ensure that incoming and outgoing flows are equal for each line and each commodity k.

Proposition 1 shows that our flow-based formulation is valid -it can be derived from the multi-commodity formulation imposed on graphs G for each line ∈ L, together with the capacity constraints [START_REF] Elbert | Freight on urban public transportation: A systematic literature review[END_REF]. Indeed, the amount of flow for commodity k traversing the arc (s i , s i+1 ) ∈ A (denoted as φ k si below) is implied by the incoming and outgoing flows (through arcs from A in and A out ). Proposition 1. In any feasible solution for the flow-based formulation (5)-( 9), ( 12)-( 17), the flow of a given commodity k ∈ K along an arc (s i , s i+1 ) ∈ A , denoted by φ k si , is determined by

φ k si = i j=1 (f in k sj -f out k sj ) (18)
where f in k sj and f out k sj are considered 0, whenever the corresponding arcs do not exist in the network G .

Proof. The proof is given by induction on i and it follows from the flow conservation constraints imposed for the multi-commodity flows on G . Clearly,

φ k s1 = f in k s1 -f out k s1
, due to the flow conservation for commodity k on line at this stop and the fact that s 1 has no entering arc from A . For i > 1, the flow conservation for commodity k on line at stop s i states that

φ k si-1 + f in k si = f out k si + φ k si . Therefore, φ k si = φ k si-1 + f in k si -f out k si = i-1 j=1 (f in k sj -f out k sj ) + f in k si -f out k si
(by inductive hypothesis)

= i j=1 (f in k sj -f out k sj ).
As in the path-based formulation, we point out that constraints [START_REF]El reparto de última milla se sube a la bici en el área metropolitana[END_REF] are valid under the assumption that there are no drop-in stops after a drop-out stop. Nevertheless, if this assumption does not hold, the formulation can be easily adapted to deal with this situation (e.g., by using Proposition 1 to determine the flow at each arc of the path of each line).

Aggregated flow-based formulation

The idea behind the aggregated flow-based formulation is to merge the flows of different lines entering (or leaving) the network from the same stop. To this end, for a commodity k ∈ K and a potential drop-in stop s ∈ S in k (resp. drop-out stop t ∈ S out k ), variable g in ks ∈ Z + (resp. g out kt ∈ Z + ) represents the flow of packages from commodity k sent from origin o k to stop s (resp. from stop t to destination d k ) to be delivered by using any of the lines passing through stop s. We note that variable g in ks represents the sum over all lines passing though stop s of variables f in k s from the disaggregated flow-based formulation introduced in Section 3.2.2, i.e., g in ks =

∈L in s f in k s . The same relation exists between variable g out kt and variables f out k t . On the one hand, this aggregation reduces the number of variables in the formulation. On the other hand, the flow passing through arcs connecting two consecutive stops cannot be deduced anymore from variables f (as this was the case for the disaggregated formulation). We therefore need to introduce arc variables φ k s used in Proposition 1 representing the flow of packages of demand k transported through arc (s, s ) for each line ∈ L and each stop s ∈ S , where s is the stop following s on line . Figure 4 gives an example illustrating the role of each variable in the model. As in the previous formulations, we additionally use binary variables x, y and z. We further define prev( , s) as the stop preceding s in line and L start s (resp. L end s ) as the set of lines having s as the first (resp. the last) stop in their path. Finally, let S be the set of stops associated with line except the last stop. The aggregated flow-based formulation for the SFOT can then be written as follows :

max k∈K t∈S out k g out kt (19) 1 2 s 1 s 2 o k t 1 t 2 t 3 d k g in ks1 g in ks2 g out kt3 g out kt2 g out kt1 s 3 s 4 s 5 φ k 1s1 φ k 2s3 φ k 2s2 φ k 1s2 φ k 1s5 φ k 2s4 φ k 1 t1 φ k 1 t2 φ k 2t1 φ k 2t3
s∈S in k g in ks ≤ P k ∀k ∈ K (20) k∈K φ k s ≤ C y ∀ ∈ L, s ∈ S (21) k∈K (g in ks + g out ks ) ≤ C s x s ∀s ∈ S (22)
g in ks + ∈Ls φ k prev( ,s) = ∈Ls φ k s + g out ks ∀k ∈ K, s ∈ S (23)
g in ks ≤ ∈L in s φ k s ∀k ∈ K, s ∈ S ( 24
)
g out ks ≤ ∈L out s φ k prev( ,s) ∀k ∈ K, s ∈ S ( 25
) ∈L φ k s ≤ g in ks + ∈L φ k prev( ,s) ∀k ∈ K, s ∈ S, L ⊆ L s \ L end s ( 26
) ∈L φ k prev( ,s) ≤ g out ks + ∈L φ k s ∀k ∈ K, s ∈ S, L ⊆ L s \ L start s ( 27
)
g in ks ∈ Z + ∀k ∈ K, s ∈ S in k ( 28
)
φ k s ∈ Z + ∀k ∈ K, ∈ L, s ∈ S ( 29
)
g out kt ∈ Z + ∀k ∈ K, t ∈ S out k (30)
x, y, z satisfy ( 5)-( 9)

Constraints ( 20)-( 22) link aggregated flow-based and binary variables and limit the flow according to demands and capacities (in a similar way as constraints ( 2)-( 4) for the path-based formulation). In addition, flow conservation constraints [START_REF] Kelly | Innovative interior designs for urban freight distribution using light rail systems[END_REF] ensure that the aggregated incoming and outgoing flows at each stop s are equal. Finally, constraints ( 24)-( 27) ensure that the resulting flow satisfies the flow conservation constraints for each individual commodity:

• Constraint (24) forbids flow from an origin to a stop s unless it is picked up by lines having s as a potential drop-in stop. In turn, Constraint [START_REF] Masson | Optimization of a city logistics transportation system with mixed passengers and goods[END_REF] states that the flow leaving stop s to its final destination must come from previous stops in the lines having s as a potential drop-out stop. In addition, the union of ( 24) and ( 25) avoids having flow not using any line (i.e., going from o k to s and from s to d k ).

We note that this was not an issue in the former models, as the flow of packages was partitioned by lines (i.e., each variable is associated with a line ), and a stop s cannot be in S in ∩ S out , for any line . However, here we use aggregated flows over lines, thus we need constraints ( 24) and ( 25) to avoid such a case.

• Constraint (26) (resp. ( 27)) states that the flow leaving (resp. entering) s by a subset of lines L cannot be greater than the flow entering (resp. leaving) s from this same set of lines plus the flow from the origin (resp. to the destination) of the corresponding demand k. Note that these constraints are not necessarily implied by [START_REF] Kelly | Innovative interior designs for urban freight distribution using light rail systems[END_REF], as the latter takes into account all flow entering (resp. leaving) s, even the flow of the lines ending (resp. starting) their paths at s (which are not counted in [START_REF] Mhl | Delivery Time Top Priority for Online Shoppers[END_REF] nor in [START_REF] Mourad | Integrating autonomous delivery service into a passenger transportation system[END_REF]). These constraints avoid the presence of transshipment at any feasible solution.

We note that line capacity constraints [START_REF] Ipc | Online shopping frequency according to online shoppers worldwide as of October 2018[END_REF] do not rely on the assumption that there are no drop-in stops after drop-out stops. Furthermore, in order to consider different capacities for each arc of a path, it suffices to replace C by the corresponding arc capacity in the right-hand-side of the constraint.

We remark that the aggregation conducted in this section comes at the cost of introducing an exponential number of constraints of type ( 26) and [START_REF] Mourad | Integrating autonomous delivery service into a passenger transportation system[END_REF]. However, for implementation purposes, these constraints can be omitted in an initial formulation, and dynamically added (in a branch-and-cut fashion) whenever an integer solution is found. Given an integer solution satisfying constraints ( 5)-( 9), ( 20)-( 30) except [START_REF] Mhl | Delivery Time Top Priority for Online Shoppers[END_REF], the separation procedure for constraints [START_REF] Mhl | Delivery Time Top Priority for Online Shoppers[END_REF] asks to find a commodity k ∈ K, a stop s ∈ S, and a subset of lines

L ⊆ L s \ L end s such that ∈L (φ k s -φ k prev( ,s) ) > g in ks . (31)
The left-hand-side of ( 31) is obviously maximized for

L = {∈ L s \ L end s : (φ ks -φ kprev(,s) ) ≥ 0}
Therefore, the separation procedure for constraints (26) can be solved in polynomial time for each k and s by calculating these values. The separation for constraints [START_REF] Mourad | Integrating autonomous delivery service into a passenger transportation system[END_REF] is analogous.

In our implementation we start the solution process by adding, for each k ∈ K and each s ∈ S, a constraint (26) (resp. a constraint [START_REF] Mourad | Integrating autonomous delivery service into a passenger transportation system[END_REF]) for each individual line in L s \ L end s (resp. in L s \ L start s ), i.e., with L = { }. We add the rest of these inequalities in a dynamic fashion whenever they are violated.

Additional inequalities

The following inequalities can be added to all proposed formulation (as they involve only the common binary variables) with the goal of improving the continuous relaxation of the models.

x s ≤ ∈Ls y ∀s ∈ S (32) y ≤ s∈S in x s ∀ ∈ L (33) y ≤ s∈S out x s ∀ ∈ L (34)
z f ≤ ∈L f y ∀f ∈ F (35)
Constraints [START_REF] Trentini | Toward a shared urban transport system ensuring passengers & goods cohabitation[END_REF] assert that a stop s is selected only if at least one of the lines passing through s is also selected (i.e., there are no isolated stops in the solution). Analogously, constraints [START_REF] Monday | E-commerce share of total global retail sales from 2015 to 2024[END_REF] and [START_REF] Zhao | Location selection of intra-city distribution hubs in the metro-integrated logistics system[END_REF] ensure that a line is selected only if at least one of its drop-in stops and one of its drop-out stops is selected. Finally, constraints [START_REF] Zhou | Freight Transport Mode Based on Public Transport: Taking Parcel Delivery by Subway as an Example[END_REF] forbid to select a fleet f if none of its associated lines L f is also selected. We remark that these constraints are not strictly speaking valid inequalities as they could be violated by feasible solutions. However, these constraints exclude solutions with redundant elements (e.g., stops for which no associated line is selected, etc.), hence we can safely add the above constraints to the formulations to improve their continuous relaxation. We observed in preliminary tests that adding these constraints actually improves the solution times of the proposed approaches.

A column generation based MIP heuristic for the path-based formulation

Although being polynomial, the number of variables of the path-based formulation may considerably grow as the size of the instance increases in the number of commodities, lines and stops. This is a strong drawback when attempting to solve real-size instances, in which even solving the continuous relaxation of the model may be challenging. In this section, we propose an alternative heuristic solution method based on the well-known column generation technique. In order to solve the linear relaxation of the model, we start with a restricted master problem (RMP) including only the binary variables, i.e., omitting the path variables w.

In this case, a solution of value 0 is obtained, i.e., no package is delivered (note that this solution is feasible).

After an optimal solution for this RMP is found, we solve the pricing problem, which asks to find a path variable w k st with a positive reduced cost. If such a variable is found, we add it to the RMP and iterate the process. If no such variable exists, we conclude that an optimal LP-solution is found. We refer the reader to [START_REF] Desrosiers | A primer in column generation[END_REF] for more details on the column generation technique. After solving the linear relaxation of the formulation, our heuristic continues by solving the full path-based MIP formulation from scratch but restricting the set of path variables only to those generated during the column generation phase. The obtained solution is therefore a heuristic one, as not all the path variables are present in the MIP. However, the optimal value obtained from the linear relaxation may serve as an upper bound, thus giving a certificate of the optimality gap of the heuristic solution.

Throughout this section we will frequently refer to some elements of the dual problem of the linear relaxation of the path-based formulation, so we now provide the complete dual problem formulation. We use variables λ ( ) I to represent the dual variable associated with the primal constraint ( ) for the corresponding set of indices I. Additionally, we use σ z , σ y and σ x to represent the dual variables associated with the upperbounding constraints of the binary variables z, y and x (appearing in the continuous formulation when the integrality constraints are relaxed). In addition to Constraints (2)-( 8), we consider the primal formulation including the additional valid inequalities (32)- [START_REF] Zhou | Freight Transport Mode Based on Public Transport: Taking Parcel Delivery by Subway as an Example[END_REF] introduced in Section 3.3. Its dual is:

min k∈K P k λ (2)
k + maxstops λ (5) + maxlines λ (6) + maxfleets λ (7)

+ f ∈F σ z f + ∈L σ y + s∈S σ x s (36) λ (2) k + λ (3) + λ (4) s + λ (4) t ≥ 1 ∀k ∈ K, ∈ L, s ∈ S in k , t ∈ S out k (37) λ (7) - ∈L f λ (8) f + λ (35) f + σ z f ≥ 0 ∀f ∈ F (38) -C λ (3) + λ (6) + λ (8)
f( ), -s∈S λ (32) s + λ (33) + λ (34) -λ

(35) f( ) + σ y ≥ 0 ∀ ∈ L (39)
-Csλ (4) s + λ (5) + λ (32) s -

∈L in s λ (33) - ∈L out s λ (34) + σ x s ≥ 0 ∀s ∈ S (40)
λ ( * ) , σ z , σ y , σ x ≥ 0 (41)

The pricing problem

Given a solution (λ, σ) of the dual of the RMP, the reduced cost of a variable w k st is given by

ck st = 1 -λ (2) k + λ (3) + λ (4) s + λ (4) t . (42)
Hence, the pricing problem asks to find values for k ∈ K, ∈ L, s ∈ S in k , t ∈ S out k such that ck st is positive and maximized, which turns to minimize the sum between brackets. Instead of performing a complete enumeration of all potential variables, we propose to tackle the minimization problem as a shortest path problem in a particular directed acyclic graph (DAG). We first develop the proposed solution approach and then further discuss its potential benefits in practice (versus a complete enumeration). Given an instance of the SFOT, we define the pricing network as the directed graph G = (V, A), where

V = {v 0 } ∪ K ∪ K S ∪ K L ∪ K T ∪ {v f }, with
• v 0 and v f are the source and sink nodes of G

• K S = {(k, s) ∈ K × S : s ∈ S in k } • K L = {(k, ) ∈ K × L : S in k = ∅ ∧ S out k = ∅} • K T = {(k, t) ∈ K × S : t ∈ S out k } and A = ({v 0 } × K) ∪ X 1 ∪ X 2 ∪ X 3 ∪ (K T × {v f }), with • X 1 = {(k, (k, s)) ∈ K × K S } • X 2 = {((k, s), (k, )) ∈ K S × K L : s ∈ S in k } • X 3 = {((k, ), (k, t)) ∈ K L × K T : t ∈ S out k } Figure 5
shows an example of the pricing network described above, illustrating all parts of the structure. For a solution (λ, σ) of the dual of the RMP, we define the cost c i of a node i ∈ V as follows:

c v0 = c v f = 0, c k = λ (2) k ∀k ∈ K c (k,s) = λ (4) s ∀(k, s) ∈ K S , c (k, ) = λ (3) ∀(k, ) ∈ K L , c (k,t) = λ (4) t ∀(k, t) ∈ K T .
We observe that each variable w k st is associated with a path from v 0 to v f in the given pricing network, namely P k st := {v 0 , k, (k, s), (k, ), (k, t), v f }, and the sum of the costs on the arcs of this path is equal to the sum in brackets in the corresponding reduced cost (42). Therefore, the pricing problem can be solved by finding a shortest path in the pricing network, using arc lengths cij := c j , for each arc (i, j) ∈ A, where costs c j are defined above.

Reducing the pricing problem to a shortest path problem instead of performing a full enumeration allows to exploit effective solution algorithms, e.g., Dijkstra's shortest path algorithm, to avoid analysing every possible path. In addition, the pricing network can be created once at the beginning (and stored in an efficient data structure) and then used for each pricing step of the algorithm. Moreover, the structure of the network is such that it can also be used implicitly (i.e., not stored in memory but traversed as needed during the search).

(k 1 , 1 ) (k 1 , 2 ) (k 2 , 3 ) (kn, 1 ) (kn, 3 ) (k 1 , s 1 ) (k 1 , s 2 ) (k 1 , s 3 ) (k 1 , s 4 ) (k 2 , s 3 ) (k 2 , s 5 ) (k 2 , s 6 ) (kn, s 1 ) (kn, s 7 ) (k 1 , t 1 ) (k 1 , t 2 ) (k 1 , t 3 ) (k 2 , t 4 ) (k 2 , t 5 ) (kn, t 1 ) (kn, t 6 ) (kn, t 7 ) k 1 k 2 kn v 0 v f

Implementation details

During the solution of the RMP, several pricing problems are solved. Each of them consists in searching for a path variable with positive reduced cost to be added to the RMP. Even if the inclusion of a single variable is sufficient to continue with the process, it is widely known that adding several variables (with positive reduced costs) may significantly improve solution times. Specifically in our setting, adding variables for different commodities may be critical to finding optimal solutions fast. Next we describe some important features of our implementation, regarding mainly the search and inclusion of path variables into the RMP.

Number of columns.

In each pricing step, we generate as many columns as possible to be added to the RMP. However, in order to control the size of the RMP, we define the following two parameters for the algorithm:

• mck: maximum number of columns related to a single commodity that can be added per pricing step.

• mc: maximum number of columns (in total) that can be added per pricing step.

Commodities' pricing. The set of commodities is analysed in an arbitrary order. Also, with the goal of avoiding arbitrary asymmetries on the number of columns generated per commodity, at each pricing step we perform individual searches for each commodity separately, always limiting to mck the number of generated columns for a given commodity. Furthermore, if at some time of one pricing step the global limit mc is reached (hence interrupting this pricing step), we store the current commodity k being analyzed and we use this commodity as a starting point for the next pricing step (shifting the order in a round-robin fashion).

Pricing strategies We implemented 3 different strategies to price each commodity. As each commodity is priced separately, all strategies are applied to the sub network associated with the corresponding commodity.

• dijkstra: we use Dijkstra's t-shortest path algorithm with t = mck on the pricing network to find the best t columns with highest reduced cost.

• first: we perform a straightforward enumeration of all potential columns and keep the first mck columns with a positive reduced cost.

• mix: we perform a straightforward enumeration and keep the first mck columns the reduced cost of which is bigger than a given threshold µ. If less than mck columns are added, we complete up to this number by keeping the best columns among the remaining ones with positive reduced cost. We remark that a variable with positive reduced cost is associated with a path on the pricing network with a cost smaller than 1. With this in mind, we fix the threshold µ = 0.5, which means that the weight of the path should be smaller than 0.5.

MIP solution time.

Our MIP heuristic consists in solving the linear relaxation of the formulation by means of column generation and then solving the MIP formulation restricting the set of path variables to those generated in the column generation phase. The whole algorithm is run within a given time limit and, especially on large instances, even solving the linear relaxation may consume the entire allocated time, causing the heuristic to end without a feasible solution. To avoid this situation, a percentage of the given time is saved as a safe time to run the final MIP. This means that if at some point of the column generation phase, the remaining time reaches the safe time, then the column generation is interrupted and the final MIP is run on this generated set of columns. We observe that when this happens, no upper bound can be deduced from the linear relaxation, as it is not solved to optimality. We further address this issue on Section 4.3.

Lagrangian bounds

We recall that the optimal value of the linear relaxation is an upper bound for the feasible solutions of the SFOT and can be used as a certificate of optimality gap of the obtained heuristic solution. However, if the column generation phase is interrupted before achieving optimality, no upper bound can be deduced. We address this issue proposing a method to derive proper upper bounds from the solutions of the intermediate RMPs. We note that for each k ∈ K, ∈ L, s ∈ S in k and t ∈ S out k , the dual constraint associated with variable w k st is constraint (37), i.e., λ

k + λ (3) + λ (4) s + λ (4) t ≥ 1. (2)
Let ( λ, σ) be a feasible (not necessarily optimal) solution of the dual of the RMP (as modelled in (36)-( 41)) and its associated objective function value be Z. This solution may not be feasible for the dual of the original formulation, as constraints (37) associated with the columns that are left out from the master may be violated by ( λ, σ). For each k ∈ K, let Zk be the maximum reduced cost associated with commodity k, i.e.,

Zk = max

,s,t

1 - λ(2) k + λ(3) + λ(4) s + λ(4) t = 1 -min ,s,t λ(2) k + λ(3) + λ(4) s + λ(4) t . (43)
Then, for each ∈ L, s ∈ S in k and t ∈ S out k , it holds by definition that

Zk ≥ 1 - λ(2) k + λ(3) + λ(4) s + λ(4) t ,
which may be stated as

( λ(2) k + Zk ) + λ(3) + λ(4) s + λ(4) t ≥ 1. ( 44
)
With these definitions, we can obtain a feasible solution ( λ, σ) for the dual of the (LP) master problem, based on ( λ, σ) by setting λ(2) k = max( λ(2) k + Zk , 0), for each k ∈ K (leaving all remaining variables unchanged from ( λ, σ)). We remark that ( λ, σ) is feasible as λ (2) k (i.e., the sole modified variable) appears only in the dual constraint associated to variables w k st , i.e., in constraint (37), and the latter is satisfied by (44). Note that by taking λ(2) k = max( λ(2) k + Zk , 0) we also prevent the non-negativity constraint of this variable being violated (in case of a negative Zk ). As ( λ, σ) is feasible for the dual of the (LP) master problem, then its objective function value

Z = Z + k∈K ( λ(2) k - λ(2) k )P k
is a valid upper bound for the master problem (we recall that P k is the right-hand-side of (2), thus the coefficient of λ

(2)
k in the dual objective function).

Computational experiments

In this section, we present numerical results on several computational experiments. We first compare the three formulations presented in Section 3.2. Secondly, we show the results obtained by the Column Generation (CG) heuristic of Section 4 when using different values for the parameters, in order to find the best configuration. We compare then the CG-heuristic approach with the three formulations of the SFOT on different criteria, such as running time and quality of the obtained solution. We close this section by performing a managerial sensitivity analysis, in order to gain further insight on the strategical issues related to FOT systems.

Our computational framework is a desktop PC running Microsoft Windows 10 on an Intel Xeon CPU at 3.7 GHz and 128 Gb of RAM. We use CPLEX version 20.1 to solve the MIP formulations. The column generation approaches are developed in Java by resorting to CPLEX's Java API. In all approaches, we use CPLEX default configuration, except for the number of parallel threads which is limited to one for a better control on the measuring of running times.

Instance generation

Gathering real data to generate realistic SFOT instances is a challenging issue, involving logistic real data (e.g., public transport lines and stops), geographical data (e.g., depots and destinations locations), realistic demand data (e.g., neighbourhood-dependent demands, etc.) and coherent infrastructure capacities (lines and stops capacities, strategical constraints limits). For this first work on the SFOT, we generate nearlyrealistic instances from synthetic data. We developed an instance generator which generates instances of the SFOT based on the following input parameters: the number n lin of lines to generate, the number n sto of drop-in and drop-out stops on each line, the number n fle of available fleets to which lines will be assigned, and the numbers n ori and n des of origins and destinations to be generated in the instance. For the sake of transparency, we describe next the key aspects of our instance generator.

Each location (i.e., stops, origins and destinations) is generated in a rectangular area composed of a subarea A c of a fictional city and a subarea A s of its surroundings. Our procedure starts by generating the n lin lines, one at a time. Each line starts from an initial stop, which is randomly generated in A s . Then, n sto consecutive drop-in stops are generated for line , towards the central zone of the city area A c (with potential random perturbations in the direction). The distance between each stop and the following stop is given by a predefined step size (with potential random perturbations). Finally, n sto consecutive drop-out stops are generated following the same direction (and step size), towards the surroundings of the city on the opposite side of the initial point. If at any time, a stop is generated close enough to a previously generated stop, then the two stops are merged into a single one (thus, allowing different lines to share stops). For each line , we generate the reverse line by associating the stops in the inverse sequence (drop-in stops as drop-out and vice versa). In terms of the models, the reverse path would be represented as a different line , independent from (although associated with the same fleet). For a better representation of real bus paths, we fix a probability q for each stop of the reverse path to be placed at a different location than the one from (in which case a new stop may be generated). In our implementation we fix q = 30%. The location of the n ori origins are randomly generated in A s and the location of the n des destinations are generated in A c . As an example of the obtained geographical distribution, Figure 6 depicts two instances constructed by our generator; squares represent origins, circles represent destinations, each color line represents a public transport line and triangles represent stops (pointing up are drop-in stops and pointing down are drop-out stops; some stops may be both drop-in and drop-out stops, although for different lines).

As far as the demands of the generated instances, for each origin-destination pair (o, d), a probability p is defined to determine whether a demand k from o to d is generated or not. In our instances, we fix p = 80%. If such demand k is generated, then the set of reachable drop-in stops S in k (resp. drop-out stops S out k ) are those stops within a radius from o (resp. d) of a predefined number of steps (in terms of the mentioned step size). In our implementation we limit this to 5 steps from origin o k for drop-in stops in S in k , and 3 steps from destination d k for drop-out stops S out k . We remark that we do not generate origins or destinations with no reachable stops, however, even when both S in k and S out k are non-empty, demand k might be unreachable if there exists no line ∈ L with both a drop-in stop in S in k and a drop-out stop in S out k , i.e., if for each line ∈ L, S in k = ∅ or S out k = ∅ . We further explore this issue in Section 5.5. All numerical elements generated for the instance (e.g., capacities, demanded packages, etc.) are randomly chosen within predefined ranges according to a uniform distribution. After some preliminary experiments, we set these ranges as follows:

• Demands: P k ∈ [1, 10]. • Capacities: C ∈ [80, 400], C s ∈ [50, 200].
• Constraints:

-maxstops ∈ [0.05 × |S|, 0.2 × |S|]. -maxlines ∈ [0.1 × |L|, 0.3 × |L|]. -maxfleets ∈ [0.3 × |F |, 0.5 × |F |].
For the purpose of this work we generated two sets of instances, with two different objectives. The first set is used to compare the computational performance of the various approaches to solve the SFOT, and it is composed of instances of different sizes: n lin ∈ {20, 40, 60}, n sto ∈ {4, 8}, n ori ∈ {5, 10} and n des ∈ {50, 100, 200, 500}. Every possible combination of these values is considered except the ones having n ori = 10 and n des ≤ 200, which we do not consider as realistic scenarios (recall that n ori refers to the number of CDCs). For each combination of the above parameters, 3 random instances are generated, for a total of 90 instances. We call this set the Group A, which is used in the next three subsections. The second set of instances, used in the managerial analysis, is described later in Section 5.5. All these instances generated for our computational experimentation will be available for public use after publication of this paper.

Empirical comparison of the three MIP models

The first set of experiments is made on the Group A instances and compares the performance of the three MIP models proposed in Section 3.2, with a time limit of 300 seconds. We call path, flow and aggflow the path-based, flow-based, and aggregated flow-based formulations, respectively. Figure 7 shows the average running times and average lower bounds (i.e., best solution found) obtained by each formulation, grouping the instances according to the number of destinations.

Focusing on running times, aggflow is clearly much more time-consuming than the other two formulations. It reaches the time limit for all instances with 200 and 500 destinations (even for 2 instances with 500 destinations which are solved in less than 10 seconds by path and flow). On the other hand, the running times of path and flow are tied; flow seems to be slightly faster on instances having up to 200 destinations, but it is outperformed by path on the group of largest instances.

When analyzing the obtained solutions (right chart on Figure 7), we can see again that aggflow is considerably outperfomed by the other two formulations. In turn, this chart shows that path obtains better solutions than flow as the size of the instances grow. For a better understanding of the differences on the lower bounds obtained by the latter two formulations, we show in Figure 8 a cumulative chart of the dual gaps obtained for all the instances in Group A. The dual gap for an instance is calculated as (ub -lb)/ub, where lb and ub stand for the lower and upper bound, respectively, obtained by the model. In the cumulative chart, for each value t on the horizontal axis, we show the number of instances in which the obtained gap is not bigger than t. For a better visualization, we cut the horizontal axis to include only instances in which the obtained gaps are smaller than 100%. As can be seen, path shows to be more efficient than flow for closing dual gaps.

Based on the results presented above, we conclude that the path formulation is the most promising, especially when the size of the instances grows. In the following subsections we explore the limits of this formulation both as an exact method and as a heuristic one, by using column generation techniques described in Section 4.

Column generation strategies

In Section 4.2 we proposed three strategies to solve the pricing problem when applying column generation to solve the linear relaxation of the path formulation, namely: dijkstra, first and mix. We recall that the pricing problem is run separately for each commodity and the routine is based on the following parameters:

• mck: maximum number of columns associated with a single commodity that can be added per pricing step.

• mc: maximum total number of columns that can be added per pricing step.

Based on the results of exhaustive computational experiments (which are shown in Appendix A), it was possible to determine the best (most stable) configuration for each of the three strategies, which is mc = 500 and mck = 10. In addition, based on the same experiments, we set the safe time parameter (see Section 4.2) to be 10% of the time limit.

Comparing the CG-heuristics with the path formulation

We now assess the performance of the GC-heuristics detailed in Section 4, i.e., solve the LP by means of CG, then solve a reduced MIP with the columns generated during the CG phase. We test the three strategies to solve the pricing problem with the best configuration settings for each of them (as described in Section 5.3).

The experiments are made on instances in Group A, with a time limit of 300 seconds.

In Figure 9 we show the average running times (left) and average values (right) of the solutions obtained by path and the three CG-heuristics, over all instances in Group A. Instances are grouped by number of destinations. It can be seen that the average running times are considerably higher for path for all subgroups of instances. Even if the path formulation obtains better objective values on small instances (i.e., up to 200 destinations), the difference against solutions obtained by the CG-heuristics is not significant (path gets improvements of 3.4%, 2.8% and 4.2% for the three first subgroups of instances, respectively). Moreover, for the subgroup of biggest instances, the quality of the solutions obtained by path decreases and it is considerably outperformed by the CG-heuristics, whose solution values are around 30% better than path.

We further analyze the groups of biggest instances and we show in Figure 10 the same data as in Figure 9 but considering only instances with 200 and 500 destinations (grouped by the number of lines). We can see again a similar behaviour both on the running times and in the quality of the obtained solutions. In terms of solution quality, we note that the largest difference is obtained on the hardest instances (i.e., instances with 60 lines), in which the average solution values found by the MIP heuristic achieve an improvement of around 50% with respect to the average values found by path.

To conclude this section, we analyze the dual gaps represented by the obtained solutions. To this end, for each individual instance, we consider the best (i.e., smallest) upper bound ub * found within the given time limit by any of the three exact MIP formulations (i.e., the optimal solution if found, or the best achieved dual bound when time limit is reached), and we use this value to calculate the dual gap obtained by each of the methods analyzed in this section. We calculate this gap as (ub * -lb)/ub * , where lb is the solution value obtained by the corresponding method. We focus on the set of biggest instances, i.e., those with 500 destinations. The left chart of Figure 11 shows the average dual gap (as defined above) obtained by each method grouping instances by number of available lines. As evidenced by the chart, the CG-heuristics produces much better solutions. Only for the set of instances with 20 lines the obtained gaps are similar to the ones of path. However, when focusing on instances with 40 and 60 lines, we can see that the gap obtained by path is significantly larger than the one obtained by the CG-heuristics. Specifically, for these two groups of instances, the path formulation obtains average gaps of 12.87% and 44.89% while the CGheuristics obtain gaps around 3.70% and 12.03% (average over the three strategies), respectively. With the goal of analyzing this data in a disaggregated fashion, the right chart of Figure 11 arranges this information in a cumulative chart. For each value t on the horizontal axis, we show the number of instances in which the obtained gap is at most t. For a better visualization, we display the horizontal axis in logarithmic scale. As can be seen, the gaps obtained by the CG-heuristics are usually smaller than the ones for path. In particular, we note that path obtains a gap of 100% for 6 of the 36 instances of this group, while the CG-heuristics achieve gaps smaller than 12% for all but one of these instances (and a gap of 100% for the remaining one), thus confirming their superiority.

Finally, we further restrict our attention on the subset of largest instances of this group, namely those with 40 and 60 lines, and we show in Figure 12 CG-heuristics outperform the exact MIP approaches.

Sensitivity analysis and managerial insights

The instances of the SFOT involve many elements and parameters that affect the solutions in different aspects, such as lines and stops capacities, strategical constraints, number of available lines, number of commodities, etc.

In this section we perform a sensitivity analysis with the aim of finding some managerial insights on different elements of the SFOT. To this end, we designed a test set of base instances, from which several scenarios are constructed by relaxing/restricting constraints or increasing/reducing resources availability. Specifically, given a base instance I, we define the modified instance I(f card , f cap , f lin , f dest ) as the instance obtained from I by doing the following modifications:

• strategical limits maxstops, maxlines and maxfleets are scaled by a factor f card ,

• stops and line capacities C s and C are scaled by a factor f cap ,

• the set of lines available for the instance is composed only by the first f lin |L| lines in L (new lines may be generated if f lin > 1), and

• the set of destinations is composed only by the first f dest |D| destinations in D with their corresponding demands (new destinations and demands may be generated if f dest > 1).

For example, given an instance I, we can create instances I(f card , 1, 1, 1), for f card ∈ {1.2, 1.6, 1.8, 2.0}, to explore the improvement in the solution quality when the strategical limits (i.e., maxstops, maxlines and maxfleets) are increased by 20%, 60%, 80% and 100%, respectively.

Our analysis aims to explore the marginal contribution of varying the four key dimensions or parameters mentioned above, namely: increasing lines and stops capacities, relaxing cardinality constraints, expanding the set of available lines and expanding the set of destinations (and its associated demands). Our test set is constructed starting from a set I base of 10 base instances generated by our instance generator setting n lin = 5, n sto = 10, n ori = 5 and n des = 100 (see Section 5.1 for references on these parameters). From each of these instances we obtain 3 additional instances for each of the mentioned dimensions, in which the corresponding dimension is increased and/or relaxed by 20%, 60% and 100%. Specifically, we generate the following sets of instances (which, for completion, include also base instances): We remark that the increments 20%, 60% and 100% are specifically chosen to yield integer values when expanding the discrete sets of available lines and destinations (i.e., I lines and I dest ), to be properly in line with the other two dimensions. Note that the first four set of instances are aimed at studying what happens by modifying one dimension at a time. The last set instead considers the modifications of all elements simultaneously.

• I card = {I(f card , 1, 1,

Effects on the solution value

Figure 13 shows the solution values when each of the considered dimensions is increased by 20%, 60% and 100%. The values shown on the chart represent the average improvement of the solution value with respect to the base instances (the latter are included and labeled as Base in the chart for reference). Each dimension is analyzed individually on the first four groups of bars by resorting to the sets of instances I card , I cap , I lin and I dest , respectivelly. In addition, we show in the fifth group of bars the results obtained when all four dimensions are increased at the same time in the set of instances I all .

As can be seen in Figure 13, the percentage of value improvement is not always in line with the percentage of increase of the parameter w.r.t. the base instances, which is not surprising in sensitivity analysis. The smallest improvements are found when relaxing cardinality constraints (group Cardinality), which seems to have almost no impact on the solution quality, beyond the first increase of 20%, which in turn yields a marginal improvement in the solution value of just 5%. This suggests that cardinality constraints are slack after this point and the limitation is given by other dimensions of the problem. The best improvements are obtained when adding more lines (group Lines), however, even when doubling the number of lines (+100%) the average improvement is just 41%. We remark that adding lines to an instance not only adds delivery capacity, but may also increase the set of reachable demands, thus increasing the potential of the objective function.

Even though individual increase of each dimension do not result in the same increase o the solution value, a remarkable fact is that when all four dimensions are relaxed at the same time (group All at once), significant improvements are obtained. In particular, doubling all dimensions at once (+100%) allows to increase the objective values by a factor 3 compared to those of the base instances.

Effects on the resource saturation levels

For a better comprehension of the limitations imposed by each of the four considered dimensions, we focus now on the level of saturation of different element of the SFOT. Specifically, we analyze the following aspects:

• Covered demand: the percentage of the reachable demand covered by the solution. Recall that for a given set of lines, not necessarily all demand is reachable (see Section 5.1). We distinguish between the covered reachable demand (referred as covered demand in what follows) and covered total demand (which may change, when e.g., additional lines are introduced).

• Lines/stops capacity: the percentage of the full capacity of the lines/stops selected by the solution which is actually used.

• Fleets/lines/stops cardinality: the percentage of the maximum number of fleets/lines/stops imposed by the strategic constraints which has been selected as part of the solution. For example, if maxlines = 4 and only 3 lines are employed in the solution, then the line-cardinality saturation level is 75%.

Figure 14 depicts the average saturation levels of the elements described above, achieved when relaxing each of the four dimensions considered in the former section. We show one chart for each dimension and finally a chart on which all dimensions are relaxed at the same time.

We observe from Figure 14a that when the maximum number of stops, lines and fleets is increased, the capacity of the network cannot be fully exploited. For example, even though we allow to double the number of fleets, the number of actually used ones only marginally increases. Indeed, after doubling the allowed cardinality, we realize that there is only a negligible increase in the covered demand (namely from 79% for the base setting, to 81% for the setting where the cardinality is doubled). This indicates that the cardinality constraints are less constraining than the other elements of the model like capacity constraints, geographical dispersion of the demand and distribution of the lines. Figure 14b depicts the features of solutions obtained when stop-and line-capacities are increased. The effect on the covered demand is much more pronounced in this case: we observe an increase from 79% of the total covered demand (base setting) to 92% (in the setting in which the allowed capacities are doubled). There is also a positive effect on reducing the number of used stops, fleets and lines. The latter reductions are only marginal, which is again due to the geographical features of the input graph, i.e., dispersion of the demand and distribution of the lines.

For decision makers concerned with maximizing the covered total demand, another important factor is the number of available lines and their geographical distribution. Indeed, as Figure 14c indicates, if the number of lines in the city is doubled, even with the same capacity and cardinality constraints, we are able to increase the covered total demand from 60% to 80%. A potentially even larger increase could be achieved, if these lines are ideally distributed in the network. However, such ideal topology is less likely to happen in practice. Concerning the covered reachable demand, we also notice an increase from 79% (base setting) to 89% (when the number of available lines is doubled).

In Figure 14d we analyze the effects caused by increasing the number of potential customers (i.e., the number of destinations). We assume that the available capacities and cardinalities are not changed, and we look at the impact of doubling the number of destinations. Not surprisingly, the covered reachable demand drops from 79% to 50%, if the number of destinations gets doubled. On the other hand, there is almost no effect on the size of the fleet and the number of lines used, which can be explained by the aggregation of demand. Indeed, more demand is getting served by the existing stops and existing lines, which can be seen from the increase in the utilization of the stop-and line-capacities.

To summarize, we observe that major effects on increasing the covered (total or reachable) demand are due to improving the available infrastructure , i.e., adding additional lines and spreading them closer to the demand points. This goes in line with the analysis provided in Section 5.5.1. Additional but much smaller increase of the covered demand can be achieved by increasing the stop and line capacities. Finally, a marginal increase is also possible by allowing to select a larger number of lines, stops or fleets.

Conclusions

In this paper we studied a new problem related to a recent logistic concept in last-mile delivery, in which freight and passenger transportation are combined. We proposed three mathematical formulations for supporting strategic decisions, i.e., deciding which public lines to use and where to locate drop-out and drop-in stops. The three models are computationally assessed on a set of instances of medium size. For the bestperforming of the three models, we then developed column generation-based MIP heuristics, capable of solving instances of realistic size. These heuristics provide high-quality solutions in acceptable times. A detailed analysis of managerial insights indicates that the topology of the public transportation network, in combination with the distribution and density of the demand, plays the most important role in maximizing the demand coverage.

We emphasize that the accuracy of the strategic model could be improved by considering the uncertainty and variability of customer demand using stochastic optimization. Our results indicate that the path-based formulation is a promising candidate for developing a two-stage stochastic model to tackle this problem.

Many open questions related to tactical and operational decisions still remain to be answered. In the strategic approach, we do not model the last leg of transportation. However, once the public transportation infrastructure is made available, one has to decide on the modes of transportation for the last leg and the associated fleet. In addition, daily paths for the package delivery have to be determined. To this end, new mathematical models and approaches need to be developed to solve these challenging combinatorial optimization problems in real time.

A Column generation strategies

In Section 4.2 we proposed 3 strategies to solve the pricing problem when applying column generation to solve the linear relaxation of path, namely, dijkstra, first and mix. In this appendix we show the results of an exhaustive computational experimentation comparing their performance under different configuration values for parameters mck and mc. We tested the 3 presented strategies on all instances in Group A (with a time limit of 300 seconds) with different configurations, setting mc ∈ {10, 100, 200, 500}, and mck ∈ {1, 5, 10}. Figure 15 reports the results obtained by each configuration for each strategy. Specifically, we show the average running time to solve the linear relaxation of path (top chart), the average number of added columns (middle chart) and the average number of pricing steps (bottom chart). These averages are taken among all instances in Group A.

The results suggest that dijkstra seems to be the fastest strategy, also reducing both the number of columns and iterations needed. The first strategy obtains fairly similar running times but it usually demand more columns and iterations to finish. We note that the running times reported for the mix strategy are significantly higher than the other two strategies. We recall that when there are not enough columns with a reduced cost greater than 0.5, the mix strategy needs to check every possible column to select the best columns within the remaining ones. This explains why significantly more time is needed for the column generation of the mix strategy, even though its number of added columns is comparable to the one of first (in which no threshold is needed to use a positive reduce column). As we will see later, even though this property may sound disadvantageous, it brings a diversity of the columns that helps us to find feasible solutions of the higher quality.

As far as the setting of mc and mck is concerned, it is clear that using mc = 10 does not yield the best results, but there seems to be just slight differences for values greater than 10. Regarding mck, all strategies seem to benefit from setting mck = 10, although the differences are not significant. Based on these results, we decide to choose the best (most stable) configuration for each of the three strategies, which seems to be to use mc = 500 and mck = 10.

C O N T A C T Prénom NOMNOMNOM

Fonction fonctionr Tel. +33 (0)1 34 43 00 00 nnomnom@essec.edu

3. 1

 1 Problem statement An instance of the Strategic Freight-on-Public-Transport problem (SFOT) is composed of several sets of elements and parameters. The city is represented by a set of origins O, a set of destinations D and a set • O, D: Sets of origins and destinations • L: Set of lines • F : Set of fleets • L f ⊆ L: Set of lines of fleet f • f( ) : Fleet of line

  represents the flow of packages of commodity k from origin o k to stop s for line (resp. from stop t of line to destination d k ). To model multi-commodity flows on G , we shall also use variables to indicate the amount of flow transported along each arc (s i , s i+1 ) ∈ A and for each commodity k ∈ K.

1

 1

Figure 3 :

 3 Figure 3: Example network illustrating the flow-based formulation variables.

Figure 4 :

 4 Figure 4: Example network illustrating the aggregated flow-based formulation variables.

Figure 5 :

 5 Figure 5: Example of a pricing network associated with the pricing problem.

Figure 6 :

 6 Figure 6: Two examples of instances created with our SFOT instance generator.

Figure 7 :

 7 Figure 7: Average times (left) and lower bounds (right) obtained by the 3 MIP formulations on instances from Group A, grouping the instances by number of destinations.

Figure 8 :

 8 Figure 8: Number of instances with a dual gap smaller than a given value (on the horizontal axis).

Figure 9 :

 9 Figure 9: Average running times (left) and average values (right) of solutions obtained by path and the three CG-heuristics.

Figure 10 :

 10 Figure 10: Average running times (left) and average values (right) of the solutions obtained by path and the three CG-heuristics (instances with 200 and 500 destinations only).

Figure 11 :

 11 Figure 11: Dual gaps obtained by each approach on Group A instances with 500 destinations. Left chart shows the average dual gaps (grouped by number of public lines); right chart depicts the number of instances having a dual gap not greater than a given value.

Figure 12 :

 12 Figure 12: Box plot diagram of the dual gaps obtained by each method on the set of instances with 500 destinations. Left chart use instances with 40 and 60 lines; right chart only those with 60 lines.

  1) : I ∈ I base and f card ∈ {1, 1.2, 1.6, 2.0}} • I cap = {I(1, f cap , 1, 1) : I ∈ I base and f cap ∈ {1, 1.2, 1.6, 2.0}}

Figure 13 :

 13 Figure 13: Average improvement on solution value when each key parameter is increased by 20%, 60%, 100%.

  (a) Relaxation of cardinality restrictions (b) Increase of capacities (c) Increase of number of lines available (d) Increase of number of destinations

Figure 14 :

 14 Figure 14: Changes on the covered demands and resource utilisation when different elements of the problem are modified.

Acknowledgements

This work was funded by CY Initiative of Excellence (grant "Investissements d'Avenir" ANR-16-IDEX-0008"). This support is greatly appreciated.