Robust transit line planning based on demand estimates obtained from mobile phones

Lee, Chungmok; Nair, Rahul Article Robust transit line planning based on demand estimates obtained from mobile phones EURO Journal on Transportation and Logistics (EJTL) Provided in Cooperation with: Association of European Operational Research Societies (EURO), Fribourg Suggested Citation: Lee, Chungmok; Nair, Rahul (2021) : Robust transit line planning based on demand estimates obtained from mobile phones, EURO Journal on Transportation and Logistics (EJTL), ISSN 2192-4384, Elsevier, Amsterdam, Vol. 10, Iss. 1, pp. 1-13, https://doi.org/10.1016/j.ejtl.2021.100034 This Version is available at: https://hdl.handle.net/10419/325148 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ Robust transit line planning based on demand estimates obtained from mobile phones Chungmok Lee a , Rahul Nair b , * a Department of Industrial and Management Engineering, Hankuk University of Foreign Studies, South Korea b IBM Research Europe, Ireland ARTICLE INFO Keywords: Robust network design Transit line planning Mobile phone data Column generation ABSTRACT The line-planning problem seeks to determine the set of fixed routes (or lines) a transit operator should run, along with associated operation frequencies. We propose an optimization algorithm for the transit-line-planning problem based on bi-level programming that exploits the problem’s structure in conjunction with the range estimation of demands. The issue of conservativeness due to the range estimation is mitigated by adopting the robust optimization approach. The model was inspired by a real-world application that leverages big data available from telecommunications operators to estimate city-wide mobility patterns. Demand estimates from such sources are based on large sample sizes (often orders of magnitude larger than those used in survey-based approaches), and capture day-to-day variability in travel demand as ranges. The validity of the proposed algorithm is demonstrated by using real-world data derived from 2.5 billion call data records from Abidjan, C^ ote d’Ivoire. 1. Introduction We present a robust bi-level program that addresses the transit-lineplanning problem with capacity constraints. The line-planning problem seeks to determine a set of fixed routes on which an operator should run along with the frequencies for each route. In real life, the travelers do not just take the shortest path to their destination, as it often involves a long waiting time for public transport such as buses and/or subways. In order to reduce travelers’total travel times, the operator needs to determine, whilst taking resource limitations in to account, which public transport routes to operate and how frequently they should operate, since the more frequently the transport operates, the less waiting time is expected. Line planning is one stage in the hierarchy of the overall transitplanning process that transit operators conduct. Several recent texts provide a context for this process, (Vuchic, 2017;Ceder, 2016) gives an overview. The first stage is the design of a fixed route infrastructure (stops, streets, tracks for rail/tram services). The second stage is the transit operator’s determination of the lines that should run on that infrastructure. A line is defined, in our context, as a fixed route to be operated during a particular time of day. More specifically, a line is an ordered sequence of stops. Each line has an associated frequency, which denotes the number of services the travelers see every hour. The subsequent stages of transit planning involve timetabling and rostering of crews and vehicles. In this paper, we focus on the line-planning problem faced by operators who have a large number of routes to operate. The travel-demand model used in this research is motivated by a realworld application for estimation of mobility demand from big data available from telecommunications operators. Demand estimations from such sources leverage large sample sizes, typically orders of magnitude larger than traditional surveys. In contrast to classical demandestimation methods that seek a point estimate between an origindestination (OD) pair for a specific time period, mobile-phone-based estimates provide a distribution that captures inherent day-to-day uncertainties in the travel-demand process. Accurate estimation of travel demand can be problematic on two counts. First, classical econometric demand-estimation methods based on surveys can be prohibitively expensive for many cities in the growth markets. Second, such surveys might not be timely, since cities are growing so rapidly the changes in travel needs often outpace the frequency at which travel surveys can be conducted. African cities, for example, are growing at an annual rate of 3–5% (Kumar and Barrett, 2008). This rapid urbanization has been accompanied by similarly rapid mobile-phone adoption. The data obtainable from mobile phones can * Corresponding author. E-mail addresses: [email protected]c.kr (C. Lee), [email protected] (R. Nair). Contents lists available at ScienceDirect EURO Journal on Transportation and Logistics journal homepage: www.journals.elsevier.com/euro-journal-on-transportation-and-logistics https://doi.org/10.1016/j.ejtl.2021.100034 Received 9 January 2020; Received in revised form 12 March 2021; Accepted 15 March 2021 2192-4376/©2021 The Author(s). Published by Elsevier B.V. on behalf of Association of European Operational Research Societies (EURO). This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). EURO Journal on Transportation and Logistics 10 (2021) 100034 offer key insights into mobility patterns; moreover, it is infrastructurefree and timely. Such data, representative of large user samples, can provide rich characterizations of travel demand. Data resolution in space and time typically enables activity-sequence reconstruction that allows for a detailed assessment of movement patterns and captures day-to-day variability inherent in the demand process. For example, the data can reveal changes in OD flows on account of special events (Dong et al., 2015) as well as people density estimation by location for security applications (Calabrese et al., 2015). A key big-data insight from telecommunications data is that repeated observations over a period of time allow for a richer characterization of mobility-demand processes. Travel demand between two zones not only can be considered as a point estimate (as a function of time of day), but also as a random variable with a known distribution. This characterization should then be exploited in the context of transit-network design which requires a novel solution approach. For optimization of the line-planning problem with respect to demand uncertainties, there are two dominant approaches, namely stochastic programming and robust optimization. In the stochastic programming, uncertain demand is represented by probability distributions according to certain estimated parameters. Based on these distributions, many demand scenarios are generated and considered in order to achieve solution robustness. Explicit consideration of several such discrete scenarios often results in a very large optimization problem. To address this problem, optimization algorithms such as Benders decomposition can be used, though the number of required scenarios often grows too quickly with the number of origin-destination (OD) pairs. By contrast, the other approach to the line-planning problem, robust optimization does not explicitly take each scenario into account. Rather, instead of discretized scenarios, a set of all (infinitely many) possible realizations of demand data, termed an uncertainty set,isdefined. The goal, simply, is to protect the solution from the worst-case realization of demand over the elements of the uncertainty set. Using uncertainty sets is less computationally demanding than the traditional scenario-based approaches, if the uncertainty sets are defined carefully, which is a very advantageous property, especially in real-world applications (Ben-Tal and Nemirovski, 2007;Bertsimas and Sim, 2004). The worst-case realization may impact not only the feasibility of constraints but also the objective function. When the coefficients of the objective function are uncertain, the so-called objective robustness seeks the best worst-case solution in the context of the max-min framework. Nonetheless, the objective functions with uncertain coefficients defined by the uncertainty set can be trivially reformulated to constraints by adding auxiliary decision variables. Several alternative measures for ensuring solution robustness to input uncertainty are available. We propose a method that models the demand uncertainty as a range. There are two reasons for using range estimates and uncertainty sets to capture robustness. First, the proposed rangebased characterization allows us to exploit the Γ-parameterization approach (Bertsimas and Sim, 2004), which enables control of the conservatism of the solutions. For realistically sized networks, competing approaches such as scenario-based approaches can result in intractably large models. Second, the apportionment of an uncertainty budget across various demand ranges generates a conservative solution that is valid for all realizations from that set. Because it is unlikely that all demands are simultaneously at their upper bound, the solution is robust against any joint-demand realizations that are encountered. Our solution approach is based on bi-level programming, because the direct single-level formulation involves nonlinear constraints. In the bilevel formulation developed, the operator’s line-planning (routes and frequencies) decisions are considered at the upper level, while the passengers’route choices are made at the lower level. For modeling of passenger route choices, we adapt the more realistic optimal strategy model originally proposed by Spiess and Florian (1989). We then extend the model to include capacity constraints, specifically, by taking into account the delays incurred from over-congested routes. This paper’s contribution to the literature is a capacity-constrained robust transit-line planning model wherein demand estimation is available as a range. We first propose a mathematical formulation for the problem by employing the robust optimization and the bi-level programming approach, which yields a so-called flow-based formulation. To address problem sizes arising in practice, the column generation method is adapted in the context of the bi-level programming by decomposing the flow-based formulation into many shortest path problems. The column generation based approach shows significantly better results in terms of computational times than the original flow-based counterpart, which enables us to tackle a real-world problem in Abidjan, C^ ote d’Ivoire. A Monte-Carlo simulation designed to show the value of the robust approach is also presented. The paper proceeds as follows. Section 2reviews the existing literature on the line-planning problem, and Section 3introduces the proposed modeling approach and solution method. Section 4presents the Abidjan usage case, and Section 5discusses the relevant computational experiments. Finally, Section 6draws conclusions. 2. Literature review There is an extensive literature on the transit-planning process and, more specifically, the line-planning problem. The problem has been studied in several contexts and with a variety of considerations, both in the objectives and constraints. Readers are referred to the review papers on the subject (GuihaireJin-Kao, 2008;Guy and Hickman, 2007;Ibarra-Rojas et al., 2015) as well as a recent survey on deterministic methods (ZanjiraniFarahani et al., 2013). We focus herein on recent studies that have considered strategic planning and/or robustness, used similar structures of bi-level models, and exploited opportunistic sensing for transit design. All of the previous methods address one or more of the following main challenges associated with transit planning: (a) computational challenges associated with practical-sized networks leading to special solution methods (Bornd€ orfer et al., 2007), heuristic (Mauttone and Urquhart, 2009) or meta-heuristic approaches (AsadiBagloee and AviCeder, 2011;Cipriani et al., 2012); (b) the large feasible search space for network design (Bornd€ orfer et al., 2007); (c) competing trade-offs between operators and users (Ceder and Wilson, 1986;Matisziw et al., 2006); (d) the balancing of design attributes such as routes, frequencies, travel directness, transfer penalties, and (e) models of user behavior and responses to designs. Bi-level formulations have been considered in several studies. At the upper level, the design decisions are made, and at the lower level, the passengers make service choices based on service offerings (e.g., routes and operation frequencies). (Szeto and Jiang, 2014) formulated a model with upper-level decisions relating to route structure and frequencies. The lower-level problem typically models passenger route choice based on optimal strategies (Spiess and Florian, 1989), and additionally accounts for congestion and in-vehicle costs. (Constantin and Florian, 1995) proposed a bi-level program formulation based on the optimal strategies model and conducted computational experiments for three (small) real-world networks. (Goerigk and Schmidt, 2017) considers line planning when passenger route choice is integrated. This contrasts with classical approaches which impose sequence - design first and passenger assignment second. The problem is formulated as a bi-level optimization problem and solved with constraint generation. This approach is conceptually similar to the one presented here. Some works have looked to integrate all levels of transit planning (Sch€ obel, 2017). To the best of our knowledge, all of the previous approaches have assumed that demand, as based on point estimation, is deterministic. Other work has considered the lower-level assignment problem in the context of variational inequalities and solved them using an agent-based simulation (Ma, 2013). Some studies have argued that solutions to line-planning problems need to account for the uncertainty of key inputs such as travel times C. Lee, R. Nair EURO Journal on Transportation and Logistics 10 (2021) 100034 2 (Yan et al., 2013) and link failures (Marín et al., 2009). Other studies, considering demand data to be uncertain, have employed a scenario-based demand-realization perspective. For example, (Amiripour et al., 2014), studying seasonal demand variations, argued that transit designs that are valid for flows during one season are likely to be sub-optimal during others. Certainly, formulations are robust only to the extent that the strategic design is valid for all scenarios considered. Utilizing a bi-level model incorporating genetic algorithm based meta-heuristics, (Ukkusuri et al., 2007) approached the problem of long-term demand uncertainty by considering a discrete scenario set to minimize the expected travel costs and the standard deviation. Unfortunately, determining the scenario set is not straightforward for planning applications that arise in practice, and moreover, the number of scenarios needed to achieve a representative sample of demand uncertainty can be quite large. From the data mining perspective, there are extensive studies on extraction of disaggregated and/or aggregate mobility patterns obtained from opportunistic sensing. Some of those researchers have sought to address transit-planning questions. For example, (Bastani et al., 2011) used GPS trajectory data from taxis to suggest minibus routes. (Chen et al., 2013) also used the same kind of data to derive OD matrices that they then used to plan nighttime services. (Calabrese et al., 2010), having conducted a detailed comparison of mobile-phone and vehicular data, showed that the trip lengths of vehicles and movements suggested by mobile phones are linearly related. (Alexander et al., 2014) compared flows derived from mobile phones, on the basis of which they suggested scaling based on census data for estimation of true demand processes. (Calabrese et al., 2011) presented methods for extraction of OD flows from mobile data and compared the respective results with gravity models. (Berlingerio et al., 2013) utilized point estimates of demand from mobile data to address the frequency-setting problem therein. In contrast to these previous studies, the presented work aims to leverage demand distributions as opposed to point estimates. These can be estimated from big data sources like large mobile-phone data samples. The subsequently derived robust optimization model, acknowledging the inherent uncertainty associated with travel-demand uncertainty due to day-to-day variability and any potential biases associated with the data source, explicitly considers demand uncertainty as a range. To the best of our knowledge, this is the first practical-solution algorithm for the transit-network design that uses the robust optimization methodology. The proposed algorithm is implicitly able to take into account infinitely many demand scenarios (defined by an uncertainty set), and so it does not suffer from the scalability issue that often appears when the uncertain demands are considered according to a set of discrete demand scenarios. 3. Model Given (a) a transit network of stops and arcs, (b) a set of current and candidate routes to operate, (c) a fleet size constraint in terms of number of vehicles, (d) estimated demand distributions for OD pairs, and (e) a model of traveler behavior, we seek to determine the set of lines to operate along with their frequencies such that the total of system-wide travel times is minimized. 3.1. Notations We consider a transit network GðV;AÞthat consists of a set of nodes V and a set of arcs A. The passenger (traveler) demands set Kis defined by pairs of origin and destination nodes. The sets of origins and destinations are denoted by S⊂Vand T⊂V, respectively. For node i2V, let Aþ i(or A i) be a set of arcs originating from (or entering to) i2V. The head and tail of an arc aare denoted by headðaÞand tailðaÞ, i.e., headðaÞ¼jand tailðaÞ¼iwhere a¼ði;jÞ. Each arc is associated with travel time (or cost) ca. There is a subset of nodes b V⊂Vat which waiting occurs (stop nodes). Any arc emanating from waiting node i2b Vto any nodes not in S[Tis called an access arc (or link). An access arc has an associated exit arc that is parallel but opposite in direction to the access arc. We let Ai⊆Aþ ibe the set of access links of a waiting node i2b V.Rrepresents a set of all routes (lines), and each route’s operating frequency is denoted by fr. Each route r2Rpasses an ordered set of stop nodes, while each stop node is connected by the access and exit arcs to the corresponding waiting node of the route. Any two consequent stop nodes for the same route are connected with a service link. Let Aac and Asv be the sets of access and service links, respectively. Sets of access links and service links for a route r2R are denoted by Ar ac and Ar sv, respectively. Similarly, waiting nodes for a Fig. 1. Sample network with four OD nodes and routes. The waiting occurs at the waiting nodes. The service links incur the travel time and delays if they are over capacity. Each route is described by a (frequency, capacity) tuple. This network gives the sets defined in Section 3.1 as follows: R¼f1;2;3;4g;b V¼ fa;b;c;dg;Aa¼fða;a1Þ;ða;a2Þg;…;A1 ac ¼fða;a1Þ;ðd;d1Þg;…;A1 sv ¼fða1;d1Þg and so on. Fig. 2. An example optimal strategy of demand s→t. The optimal strategy is essentially a union of three distinct routes. C. Lee, R. Nair EURO Journal on Transportation and Logistics 10 (2021) 100034 3 route r2Rare denoted by b Vr. Any access or service link a2Arof route r2Rhas operation frequency fasuch that fa¼fr, i.e., fa¼frfor all a2 Ar. A service link a2Asv also has a capacity ha(the capacity of the vehicle operating on the arc). Let Cdenote the total number of vehicles available. Fig. 1 illustrates an example network with various types of arcs and nodes. Additional notations are defined in the relevant sections of the text. 3.2. User behavior model We assume that travelers employ the so-called optimal strategy in navigating the transportation network (Spiess and Florian, 1989). The idea behind the optimal strategy is that by considering multiple routes, the waiting time can be reduced even if the routes in consideration are not the shortest ones. InthesamplenetworkshowninFig. 1, consider a traveler from sto twhose optimal strategy states that it is best to consider two routes (route 1 and route 2) at the same time to reduce the waiting time at node s. If the traveler took route 2 at node s, (s)he should transfer to either route 3 or route 4 at node cto get to node t. This optimal strategy is illustrated in Fig. 2, which is essentially a union of three paths between sand t.Thefirst path is s→a→a1→d1→d→t;the second is s→a→a2→b2→c2→c→c3→d3→d→t,andthethirdis s→a→a2→b2→c2→c→c4→d4→d→t. It is important to note that theoptimalstrategymightinvolvemultiplepathsthatarenotnecessarilythe shortest in terms of travel time. Multiple paths can reduce the waiting time, because the traveler can take whatever transport arrives first at the waiting node. Besides the waiting times, a traveler also needs to take into account the capacities of routes. The travel time of an over-congested route can be considerably longer, which results in a sub-optimal strategy. Therefore, a traveler, in estimating the possible delay, makes an implicit assumption about how many other travelers will share the same route. As it is impossible for a traveler to precisely predict the congestion levels of the routes, establishing the optimal strategy using a point estimation of the traffic level might yield a strategy that will perform poorly under the different realizations of traffic conditions. This motivates a need for a robustly optimal strategy that is guaranteed to be optimal against fluctuations in traffic demands. To derive the robustly optimal strategy, first we need to formally characterize the demand uncertainty. 3.3. Demand characterization We assume that the demand between any two nodes is estimated by the following range. Definition 1. For any s2Sand t2T,~ dst 2½ bdst dst ;bdst þdst , where bdst dst The underlying demand uncertainty and systematic biases from opportunistic sensing are likely to cause some OD pairs to show greater variance than others. For example, areas that include stadiums or large institutional buildings like hospitals are associated with greater variations in demand, which can be modeled by a wider range. Examples of this are shown in the box plot, Fig. 4, wherein changes in variability for different OD pairs can be seen. It is unrealistic to suppose that all demands are at the upper bounds simultaneously; therefore, the model makes no such assumption. Rather, a parameter Γis introduced to account for the degree of demand variability. The possible demands as constrained by Γeventually determine the degree to which the obtained solution should be robust; this explains why the following definition is often called the budget of uncertainty (Bertsimas and Sim, 2004). Note that the problem becomes a deterministic case when Γ→∞because it implies that all demands are at the upper bounds at the same time. Definition 2. Budget of uncertainty: For a given Γ, we want the solution to be guaranteed valid for any demand realization in a set D:¼8 > > < > > : d2RjSjjTj þ dst ¼bdst þdstγst ; P s2S;t2T γst Γ; 1γst 1;8s2S;t2T 9 > > > = > > > ; (1) where bdst and dst are the nominal and maximum deviation demands, respectively. 3.4. Model development Naturally the problem can be seen as a two-player game: operator vs. travelers. The operator decides which lines to operate for a given fleet of transports, where the number of vehicles operating and round trip time tr on a given route reventually determine the frequency of the route fr.A transport can be a bus, subway, tram, or anything deployable for the public transportation. Then the feasible set of frequencies for a given set of routes Rcan be stated as F:¼8 > > < > > : f2RjRj þ X r2R tr 60frC; frfrfr;8r2R: (2) The first constraint ensures that the total number of vehicles in travel is not greater than the size of the fleet C. We assume that the frequency of a transport should be within fr;frfor any r2R, which is guaranteed by the second constraints. When fr¼fr¼0, it means the route ris not operated at all. Therefore, the frequency setting fdetermines which routes to employ and how frequently the routes in use will operate at the same time. We assume the operator is a public sector body interested in social welfare and not profit maximization. The operator goal, in this case, is to find the optimal frequency setting f2Fthat minimizes the system-wide travel time. That objective function can be stated as Fðv;w;qÞ:¼X s2S 0 B B B B @X a2A caX t2Tbdstvst a |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} journey time þX i2bV ws i |fflffl{zfflffl} waiting time 1 C C C C AþβX a2Asv caqa: |fflfflfflfflffl{zfflfflfflfflffl} delay (3) The total travel time Fðv;w;qÞis the sum of the total journey time, the total waiting time, and the total delay time due to over-congested service links. The parameter βcontrols the trade-off between the travel time and the congestion-incurred travel delay. Traffic volume vst arepresents the fraction of flow from source node sto destination node ton arc a. Other terms, for example those related to operational costs, can also be included. However at this stage of transit planning, these costs are generally not known precisely without full roster information on crews and schedules. Note that the total journey time and total waiting time are considered with the nominal demands. This study focuses on the impact of demand uncertainty on the capacity because, without the capacity constraints, the optimal strategy of the individual traveler remains the same at higher demands. Moreover, the objective function implicitly depends on the frequency setting: the waiting time wis a function of the frequency setting because the frequency and optimal strategy uniquely determine the waiting time at the nodes. Then, the feasible set of flow v, waiting time w, and delay time qfor a given transport frequency fcan be stated as C. Lee, R. Nair EURO Journal on Transportation and Logistics 10 (2021) 100034 4 Yf:¼fðv;w;qÞ2RjAjjSjjTjþjbVjjSjþjAsvj þ X a2A i vst aX a2Aþ i vst a¼8 > < > : 1;if i¼s; 1;if i¼t; 0;otherwise: ;8i2V;s2S;t2T; X t2Tbdstvst aX r2R δr afrws tailðaÞ;8a2Aac;s2S; X s2SX t2Tbdstvst aþmax Ps2S;t2Tγst ¼Γ; 0γst 1;8s2S;t2T 8 > > > > < > > > > :X s2SX t2T dstvst aγst 9 > > > > = > > > > ; haX r2R δr afrþqa;8a2Asvg where δr ais 1 if a2Ar, 0 otherwise. The above represents a traveler’s set of feasible behaviors under a given operating frequency f. The first set of constraints, which is often referred to as flow-balance, ensures that all demand from the origins to destinations is served, i.e. no demand is lost. The second set of constraints relates the waiting time for a demand sat a node ito the frequencies of vehicles to ride. The last set of constraints measures the delay time on arc aaccording to the total sum of flows and the total capacity of the arc. Note that the capacity constraints are protected by the sum of worst-case demands defined by the budget of uncertainty. Due to this protection term, the traveler’s behavior becomes robust to the variation of demands. In other words, the traveler is guaranteed to travel the optimal strategy route with a delay time less than or equal to qa, which yields the robustly optimal strategy. It is well-known that the last set of constraints defining Yfcan be stated as below by using the reformulation trick from Bertsimas and Sim (2004): X s2SX t2Tbdstvst aþΓyaþX s2SX t2T ust ahaX r2R δr afrþqa;8a2Asv; yaþust adstvst a;8a2Asv;s2S;t2T; ya0;8a2Asv; ust a0;8a2Asv;s2S;t2T which is a system of linear inequalities. We denote ~ Yfas Yfwith replacing the last constraints with the reformulated ones presented above. Let PYfand ~ PYfdenote the polyhedra of Yfand ~ Yf, respectively. It is clear that the projection of ~ PYfonto the ðv;w;qÞ-space is equivalent to PYf. Here, delays are assumed to grow linearly with link volumes. This approach to handling delays is close to the “line-specific overload delays” strategy proposed by Hing-Keung Lam et al. (1999), albeit in an equilibrium setting. Several other alternative approaches to handling the delays have been proposed in the literature, including: (a) a convex optimization by changing the link cost to a function of the flows, i.e., ca¼ caðvaÞ(e.g., (Spiess and Florian, 1989)); (b) making the waiting time dependent on link flows (e.g., (De Cea and Fern andez, 1993) modeled route sections and used waiting links for each route section); (c) the “effective frequency”approach by which crowded buses have f→0 (since it is harder for passengers to board them; when transports are empty (low congestion) the effective frequency becomes the same as the nominal frequency), and (d) deploying more advanced schemes, wherein capacity constraints →higher delays →higher boarding/dwell times → lower frequencies. The optimal strategy under the given frequency fminimizing the travelers’travel times can be obtained by solving min ðv;w;qÞ2Yf fFðv;w;qÞg:(4) 3.5. Mathematical formulation We now present two alternative formulations to solve the lineplanning problem. The first one aims to jointly determine frequencies of services and flows, i.e., it makes the operator and traveler decisions jointly. The second formulation is a bi-level program that naturally separates those decisions. 3.5.1. Nonlinear programming formulation The optimal line planning problem can be stated as: ðP0Þmin X s2S X a2A caX t2Tbdstvst aþX i2bV ws i!þβX a2Asv caqa(5) s:t:X a2A i vst aX a2Aþ i vst a¼8 < : 1;if i¼s; 1;if i¼t; 0;otherwise: 8i2V;s2S;t2T;(6) X t2Tbdstvst aX r2R δr afrws tailðaÞ;8a2Aac;s2S;(7) X s2SX t2Tbdstvst aþΓyaþX s2SX t2T ust ahaX r2R δr afrþqa;8a2Asv;(8) yaþust adstvst a;8a2Asv;s2S;t2T;(9) ya0;8a2Asv;(10) ust a0;8a2Asv;s2S;t2T;(11) X r2R tr 60frC;(12) frfrfr;8r2R;(13) vst a0;8a2A;s2S;t2T;(14) ws i0;8i2b V;s2S;(15) qa0;8a2Asv;(16) where δr ais 1 if a2Ar, 0 otherwise. The objective minimizes the total system travel time and is the same as Equation (3). Constraints (6)–(11) are from the user-behavior model Yfwith the reformulation presented in the previous section. Constraints (7) relate the frequency of vehicles f (operator decision) to the waiting time w(traveler-side decision). The nonlinearity of these constraints makes solving the problem difficult. One possible approach is using linearization techniques such as RLT (Sherali and Shetty, 1980;SheraliAmine, 1992) and McCormick linearization (McCormick, 1976). However, those reformulation techniques would involve the adding of (possibly exponentially) many constraints and variables, while still not being able to guarantee the optimality of the solutions, because constraint (7) is not convex. Constraints (12) and (13) ensure the frequencies are feasible as defined in (2). 3.5.2. Bi-level programming formulation The two-player interpretation of the above problem makes the following alternative, bi-level-programming-based formulation possible: ðP1Þmin f2F GðfÞ s:t:GðfÞ¼ min ðv;w;qÞ2Yf Fðv;w;qÞ: C. Lee, R. Nair EURO Journal on Transportation and Logistics 10 (2021) 100034 5 Here, the optimal value function GðfÞis defined by an optimal solution of a lower-level optimization problem from the user perspective. The upper-level denotes the transit-operator-determined frequencies. The lower-level problem asserts the given frequencies by evaluating the total travel time and penalty for over-capacity on the link. Bi-level programming problems are generally NP-hard (Colson et al., 2005;Chinchuluun et al., 2009). Hereafter, we denote the upper-level problem as (P1) U while the lower-level problem as (P1) L . The problem (P1) can be seen as a decomposition of the problem (P0) so that both problems (P1) U and (P1) L are linear. Therefore, for the lower-level problem, the frequency setting f2Fshould be given from the upper-level problem, and the optimal solution of the lower-level problem would provide an evaluation of the goodness of the frequency setting to the upper-level problem. 3.6. Solution approach for Bi-level programming formulation The solution approach is summarized as follows. First, we set an initial frequency setting f0; then, a gradient direction to reduce the total travel time is identified by solving the traveler-side problem (P1) L , and the frequency fis updated using the descent direction and step size. This procedure is repeated until the frequency solution converges. It is noteworthy that the most computationally demanding step is obtaining the descent direction, which involves solving a large-sized user-behavior problem for the robustly optimal strategy. We propose the column generation based approach, which enables us to benefit from demand-wise decomposition and the use of very fast specialized algorithms for the shortest path problem. More detailed explanations follow below. Algorithm 1. Frank-Wolfe algorithm for (P1) We consider Algorithm 1 to solve (P1). Let b σ be a descent direction of the optimal value function GðfÞat a given frequency bf. In order to use Algorithm 1, the b σ should be known for any f2F, which is given by Theorem 1. For a given (or fixed) f, let π ,ϕ, and μ be the dual variables associated with constraints (7), (8), and the remaining constraints of the problem (P1) L , respectively. Theorem 1.For any f 2F, let b w and ðb π ;b ϕÞbe the primal and dual optimal solutions of (P1) L , respectively. Consider a vector. b σ :¼"X s2SX a2Ar ac b ws tailðaÞb π s aX a2Ar sv hab ϕa#r2R : Then,b σ is a descent direction if b σ 6¼ 0. Moreover,the descent direction b σ becomes the negative of the gradient,i.e.rGðfÞ,if the primal an dual optimal solutions are unique,where rGðfÞis the gradient of GðfÞat f 2F. Proof.We let b Yf⊂Yfand b Zfdenote the set of primal and dual optimal solutions,respectively,of the problem minðv;w;qÞ2YfFðv;w;qÞfor the given f 2 F. For a given f 2F,we consider the Lagrangian function Lðf;v;w;q; π ;ϕ; μ Þfor the problem minðv;w;qÞ2YfFðv;w;qÞ.Because,for a given f 2F the problem minðv;w;qÞ2YfFðv;w;qÞis an LP satisfying the constraint qualification,as shown in Still (2018) the directional derivative D σ GðfÞin the direction σ such that Pr2R σ 2 r¼1at the point f is given as D σ GðfÞ¼ min ðv;w;qÞ2bYf max ð π ;ϕ; μ Þ2bZfX r2R σ r ∂ Lðf;v;w;q; π ;ϕ; μ Þ ∂ fr : Note that Lðf;v;w;q; π ;ϕ; μ Þcan be seen as a function of f,which gives the partial derivatives of Lðf;v;w;q; π ;ϕ; μ Þ.For any r 2R,we have, ∂ Lðf;v;w;q; π ;ϕ; μ Þ ∂ fr ¼X s2SX a2Ar ac ws tailðaÞ π s aþX a2Ar sv haϕa0 because π s a0and ϕa0. Now we consider a pair of primal and dual optimal solutions ðbv;b w;bqÞ2 b Yfand ðb π ;b ϕ;b μ Þ2b Zf.We claim that b σ 2RjRj þdefined as below is a descent direction if b σ r6¼ 0: b σ r:¼ ∂ Lðf;bv;b w;bq;b π ;b ϕ;b μ Þ ∂ fr ¼ X s2SX a2Ar ac bws tailðaÞb π s aX a2Ar sv hab ϕa;8r2R: Then,we need to show that Db σ GðfÞ¼ min ðv;w;qÞ2bYf max ð π ;ϕ; μ Þ2bZfX r2Rb σ r ∂ Lðf;v;w;q; π ;ϕ; μ Þ ∂ fr <0; which implies that b σ is a descent direction.Due to (B), it is clear that Db σ GðfÞ0. Assume,for a contradiction,that Db σ GðfÞ¼0. Then,there must be ~ r2R such that b σ ~r>0 and ∂ Lðf;bv;bw;bq;~ π ;~ ϕ;~ μ Þ ∂ f~ r ¼0;for some ð~ π ;~ ϕ;~ μ Þ2b Zf; which implies that ~ π s a¼0for all s 2S and a 2Aac ~ rsuch that b ws tailðaÞ>0. Note that Lðf;bv;b w;bq;~ π ;~ ϕ;~ μ Þhas terms ~ π s a X r2R δa ~rf~rb ws tailðaÞX t2Tbdstbvst a! for s 2S and a 2Aac ~ r.Because ~ π s a¼0for all s 2S and a 2Aac ~ rsuch that b ws tailðaÞ>0, there should be ð~ π ;~ ϕ;~ μ Þ2b Yfsuch that Lðf;~ v;~ w;~ q;~ π ;~ ϕ;~ μ Þ< Lðf;bv;b w;bq;~ π ;~ ϕ;~ μ Þ,which means that ðbv;b w;bqÞ 62 b Yf.This gives a contradiction. Note that (A)becomes D σ GðfÞ¼X r2R σ r ∂ Lðf;bv;b w;bq;b π ;b ϕ;b μ Þ ∂ fr C. Lee, R. Nair EURO Journal on Transportation and Logistics 10 (2021) 100034 6 when the optimal solution ðbv;b w;bq;b π ;b ϕ;b μ Þis unique,i.e., b Yf¼b Zf¼1, which implies that b σ ¼rGðfÞ. Theorem 1 implies that b σ at any f can be obtained easily from the primal and dual optimal solutions of (P1) L .Note that the descent direction may depend on the choice of a pair of primal and dual optimal solutions ðbv;b w;bqÞ2 b Yfand ðb π ;b ϕ;b μ Þ2b Zf.For this study,we take the primal and dual optimal solutions of the problem (P1) L solved by the dual Simplex method.Then,the direction finding problem at f lcan be stated as follows. ðDFÞmin X r2Rb σ rfr(17) s:t:X r2R tr 60frC;(18) frfrfr;8r2R:(19) The above is an LP problem that can be solved easily by LP solvers. Moreover,the problem is also a variant of the fractional Knapsack problem which can be solved in OðjRjlogjRjÞ as shown in the following theorem. Theorem 2.(DF) can be solved in OðjRjlogjRjÞ. Proof.Since (DF)is a fractional Knapsack problem,an optimal solution is obtained by selecting the maximum possible frfor each r 2R in decreasing order of b σ tt=60,where the sorting is done in OðjRjlogjRjÞ. For the step length finding problem,using discretized values between 0and 1for α ,we solved the problem min α 2½0;0:1;0:2;…;1G α f*þð1 α Þfl;(20) whose optimal solution is denoted as α *,where f*is the optimizer of (DF). This implies that we need to solve (P1) L many times.Recall that problem (P1) L is to find an optimal strategy with the capacity constraints under demand uncertainty,which minimizes the total travel time,for a given frequency f.If there are no capacity constraints,the problem (P1) L is separable by the source nodes s2S,which results in jSj-independent problems for the optimal strategies.For these problems,we may use a dynamic programming algorithm developed by Spiess and Florian (1989),which is very fast.Unfortunately,the problem is no longer separable with capacity constraints,so their approach cannot be used. On the other hand,the problem (P1) L can be stated by removing (12)-(13) and fixing f in (P0). According to the linear programming duality,the inner maximization problem in constraints (8) can be reformulated (Bertsimas and Sim, 2004). Then,we have the following formulation,called the “flow-based” formulation. ðGÞminð5Þ s:t:ð6Þ;ð7Þ;ð14Þ;ð15Þ;ð16Þ; X s2SX t2Tbdstvst aþΓyaþX s2SX t2T ust a haX r2R δr afrþqa;8a2Asv;(21) yaþust adstvst a;8a2Asv;s2S;t2T;(22) ya0;8a2Asv;(23) ust a0;8a2Asv;s2S;t2T:(24) Note that (G)is a linear programming problem with a fixed f.However,the size of the problem grows rapidly with the size of the network and especially with the number of demands.Specifically,each demand pair ðs;tÞincreases the size of the problem by jAjþjAsvj.A typical transit network should consider a large number of demands;consequently the size of the problem may be too large to solve as it is.To overcome this,we propose a column generation based approach in the next section. 3.6.1. Column generation for solving (G) For a given source and sink pair st, any feasible solution of (G) can be represented as a union of paths between sand tdue to the flow balance constraints (6). Let Pst be a set of all paths between sand t. Then it is easily seen that the following relation vst a¼X p2Pst σ p azst p; holds, where σ p a¼1 if path p2Pst uses arc a, 0 otherwise. We introduce new extended variables 0 zst p1 representing the fraction of demand between sand tusing path p. By using the above relation, we consider the extended formulation. ðEGÞmin X s2S X a2A caX t2T dst X p2Pst σ p azst pþX i2bV ws i!þβX a2Asv caqa(25) s:t:X p2Pst zst p¼1;8s2S;t2T;(26) X t2T dst X p2Pst σ p azst pX r2R δr afrws tailðaÞ;8a2Aac;s2S;(27) X s2SX t2Tbdst X p2Pst σ p azst pþΓyaþX s2SX t2T ust ahaX r2R δr afrþqa;8a2Asv;(28) yaþust adst X p2Pst σ p azst p; 8a2Asv;s2S;t2T;(29) ya0;8a2Asv;(30) ust a0;8a2Asv;s2S;t2T;(31) zst a0;8a2A;s2S;t2T;(32) ws i0;8i2b V;s2S;(33) qa0;8a2Asv:(34) Since there can be exponentially many paths, we incorporate the column generation based approach. Consider a column generation master problem with a limited set of paths. Let ðbθ;b π ;b ϕ;bλÞbe the dual optimal solution of the column generation master problem, where θ, π ,ϕ, and λ are the dual variables associated with (26), (27), (28), and (29), respectively. The column generation master problem is optimal if the condition bdst0 @X a2AðpÞ caX a2AðpÞ\Aac b π s aX a2AðpÞ\Asv b ϕaþdst bdst X a2AðpÞ\Asvbλst a1 Abθst 0; 8p2Pst;s2S;t2T; holds, where AðpÞ⊆Ais a set of arcs in path p. Then, the column generation sub-problem for demand stis given as follows: C. Lee, R. Nair EURO Journal on Transportation and Logistics 10 (2021) 100034 7 ðSPÞξst ¼min X a2AðpÞ ~ ca(35) s:t:p2Pst;(36) where ~ ca¼ 8 > > > > > > > < > > > > > > > : cab π s a;if a2Aac; cab ϕaþdst bdstbλst a;if a2Asv; ca;otherwise: Note that ~ cais nonnegative for any a2Aprovided that ca0 since b π 0 and b ϕ0. A new path p*is added to the master problem if dst ξst < bθst , where p*is an optimizer of ξst . Since every cost of the graph for the column generation problem is nonnegative, very fast shortest path algorithms such as Dijkstra’s can be used. In practice, it is sufficient to solve the column generation subproblem jSjtimes, because the Dijkstra’s algorithm will provide all shortest paths from a source to all destinations. 3.6.2. Early termination of column generation Algorithm 1 entails solving (EG) many times. Specifically, at each iteration, we solve problem (EG) twelve times by changing f(one for obtaining ~ rGðflÞand eleven for the step-length search, since the steplength problem is solved for discrete values in the unit interval). Therefore, we keep all generated columns so that the subsequent column generation will require only a small number of new columns. Moreover, the column generation can be terminated early if the condition ~ ZEG þX s2S;t2Tdstξst bθstb Z;(37) holds, where ~ ZEG and b Zare the objective value of the current column generation master problem and the best objective value (i.e., incumbent solution), respectively. The validity of the above condition is easily seen from the so-called Lagrangian bound (i.e., the optimal value of a column generation master the current objective value of the column generation master þthe minimum of reduced costs). 4. The case of Abidjan, C^ ote d’Ivoire The mobile-phone data, made available by Orange (2014), consisted of 2.5 billion call-detail records (CDRs) for C^ ote d’Ivoire. A CDR is generated each time a call or a text-message transaction occurs. Around the city of Abidjan, this dataset consisted of anonymized CDRs for half a million users over an arbitrary two-week period within a five month interval. Each CDR contains a time and location when and where a call/text transaction occurred along with an associated ID. The spatial component is encoded as the mobile-phone antenna that handled the transaction: typically, but not necessarily, the closest antenna to the user. In the city of Abidjan, a set of 407 antennas are considered. These towers, which are analogous to Traffic Analysis Zones (TAZ), serve as the reference for production and attraction zones for travel. Each individual in the dataset therefore represents a sequence of time-space events between these zones. CDRs by themselves are mostly infrequent snapshots of disaggregated behavior. (Berlingerio et al., 2013) demonstrated, for the same dataset, how such infrequent snapshots can be transformed into (partial) activity sequences. The main intuition presented therein, is that repeat observations at the disaggregated level allow for activity-sequence reconstruction, even if the constituent daily activity patterns are partial. At the aggregate level, this variability detected in individual behaviors along with their likelihood, translates to demand, as a distribution, between any two points on the network. The same OD matrix derived by Berlingerio et al. (2013) is used here, with the exception that the range values are specified. Readers are referred to their paper for the data-extraction algorithms used. Fig. 4 shows the variability of the top 100 OD pairs during peak periods along with a box plot. As shown in Fig. 3, the data during the peak periods have a higher spread. In order to specify the random demand variable ~ dk, as a range, for each of the 107,761 OD pairs derived from the CDR data, the median of the CDR-based flows is considered as bdkand one standard deviation as dk. In practice, some scale up of demand as a function of mobile operator market share as well as population statistics such as those suggested in the literature can be considered (Alexander et al., 2014). From the supply-side, Sub-Saharan cities across Africa, with the exception of Lagos, have seen organized public transport deteriorate over the past few decades (Gwilliam, 2011). Due to cost inefficiencies and higher fares, public operators have ceded ground to informal operators who have the majority of the transport mode share. In Abidjan, a city with 4.5 million citizens, the public operator SOTRA has a fleet of 439 buses. This is complemented by roughly 5,000 mini-buses and roughly 11,000 shared taxis (woro-woro). Data on the 92 SOTRA routes were considered as the baseline network shown in Fig. 5b. The frequencies of these services were extracted from their public website. Additionally, a set of future potential routes, called the candidate route set, were considered, as shown in Fig. 5c. Roughly half of those routes were generated based on highdensity corridors as observed in the CDR data (Berlingerio et al., 2013); the other half of the candidates were generated by a shortest path algorithm between zones with high volume of passenger flows. Fig. 3. Flow distribution for top 100 OD flows by hour of day. Fig. 4. Box plot of passenger flow distribution for top 100 OD flows (inset shows top 1000 OD flows). C. Lee, R. Nair EURO Journal on Transportation and Logistics 10 (2021) 100034 8