# The Integration of Multimodal Networks: The Generalized Modal Split and Collaborative Optimization of Transportation Hubs

Authors: Yifei Cai; Jun Chen; Da Lei; Jiang Yu
Journal: Journal of Advanced Transportation
 (2022)

---

## Abstract

Integrated transportation is one of the most important methods to encourage the modal shift from car to public transportation (PT). However, as most cities have an existing multimodal network, it is difficult to expand the current networks by building more PT routes. Thus, integrating different modes through the optimization of hubs is a cost-efficient way to promote sustainable mobility. This paper develops a bilevel multimodal network design problem based on the collaborative optimization of urban transportation hubs. The upper-level problem is formulated as a mixed-integer nonlinear program to achieve a modal shift from congested subnetworks to underutilized subnetworks to realize a balanced use of the entire network. The decision variables are classified into location-based (hub locations) and route-based (route layouts and frequency setting) ones. The lower-level problem is a generalized modal split/traffic assignment problem (GMS/TAP), which captures the mode choices of all modes in the path set. The GMS/TAP is formulated as a nonlinear optimization problem (NLP) and is solved using a hybrid method of the successive average (MSA) algorithm. A hybrid genetic search with advanced diversity control (HGSADC) is developed to solve the bilevel model, where the exploration of the search space is expanded using the biased fitness function and diversification mechanism. The solution properties of the hybrid MSA and HGSADC are demonstrated in two modified nine-node networks. The model performance is illustrated in a real-size network in Jianye district, Nanjing. 9.2% decrease of travel time, 25.7% increase of service level, and a significant modal shift from car to PT are obtained.


---

## Body

## 1. Introduction

The rapid development of technologies has led to the emergence of various new transport modes and mobility services. The multimodal transport network, which consisted of these new forms of urban mobility and traditional transport, plays an overriding role in urban development [1, 2]. However, the spatial in-coordination between the demand and multimodal capacity leads to some negative effects, such as traffic congestion, air pollution, and social equality problems, that makes the system unsustainable [3]. Concerning these problems, cities have realized that the key issue to solve the problems is to encourage a modal shift from car to public transportation (PT) [4]. Consequently, integrated transport has received a lot of attention from policymakers and researchers in recent 30 years because of its potential to foster sustainable travel behaviours [5, 6].As most cities have an existing multimodal network, it can be time and cost consuming for governments to create the changes in the current network [7]. From the infrastructure perspective, there are two ways to promote the attractiveness of PT. One is to provide more direct routes to reduce transfers [8]. The other is to connect the routes with one another to create access a wider range of destinations [9]. Extended to the entire network, for cities with mature multimodal networks, it would be cost-effective to integrate all the travel modes through hubs rather than building more PT routes [10]. Urban transportation hubs, such as parking lots, metro stations, bus stops, and bicycle stations, are the places where transfers occur. They play the key role in integrating the multimodal network [11, 12]). Thus, this paper focuses particularly on integrating the physical multimodal transport network via the optimization of transportation hubs to achieve a sustainable modal shift.In the multimodal network, private car and conventional PT services are often considered as the major travel modes, complemented by shared, on-demand services and other travel modes [13]. They all have their own optimal ranges and advantages. The integration aims to provide the users with all the available travel options to develop a customizable door-to-door combined trip that could increase the attractiveness of PT [14]. Recently, mobility as a service (MaaS) is considered as an evolutionary continuation in terms of integrated transport, and it is often described as a personalised user-oriented platform that unifies the creation, purchase, and delivery of door-to-door passenger transport solutions integrating traditional and innovative modes of transport [2, 15, 16]. Similarly, the success of MaaS depends on a shift in behaviour away from reliance upon private car ownership [17].However, the effects of the combined mobility services are often questioned by the recent research works and applications [18], that MaaS should be complementary rather than in substitution to a private car. Even if in principle it may reduce the car usage of MaaS users, it may imply bigger use of the private car by non-MaaS users because of the crowding of transport services by MaaS users [19]. This reminds us that the sustainable mobility not only relies on a modal shift from car to PT but also on a balance between the demand and multimodal capacity. Since the effect of the integration did not happen immediately but gradually over the years, comparing with the short-term modal shift stimulated by new policies or reforms, the balanced use between all travel modes can keep it remain motivated [20, 21]. Thus, the aim of the integration via hubs in this paper is extended to achieve a reasonable modal shift to promote a balanced use of all travel modes.From a user’s perspective, they have many mode choices to finish their trips, including single modes and combined modes. A combined trip is usually less competitive than a single-mode trip, considering the transfer penalty. A user will only choose a combined travel path when the congestion cost of a single-mode path is higher than the transfer cost of the combined trip [22]. Thus, integrating a multimodal network by optimizing the hubs is a process that provides users with competitive combined travel choices by reducing transfer penalties. The mixed choices among all modes of all users can realize a total modal shift from congested subnetwork to underutilized subnetwork while cutting down the usage of the private car.In summary, to achieve a modal shift from congested subnetworks to underutilized subnetworks and promote a balanced use of the entire network, this study conducted the integration of the multimodal network by collaboratively optimizing the hubs of different subnetworks. This problem is referred to as the integrated multimodal network design problem (IMNDP) in this paper, focusing on the two key attributes of IMNDP, the combined modes and optimization of hubs. A new bilevel mathematical formulation is proposed based on the network design problem in this paper. To handle with the numerous modes options the network can provide, a generalized modal/split traffic assignment problem (GMS/TAP) is proposed in the lower-level problem. Meanwhile, for integrating the multimodal network through a collaborative optimization of transportation hubs, the methodologies for the assessment and optimal design of a hub integration scheme are provided in the upper-level problem.The remainder of this paper is organized as follows. The literature review and research gaps are introduced in Section2. Section 3 introduces the formulation of the lower-level GMS/TAP and upper-level problem along with their solution methods. Section 4 verifies the solution properties of the lower-level and upper-level problems along with the application in a real-size network. The conclusions are finally presented in Section 5.


## 2. Literature Review

Most of the previous researches on the integrated transportation focus on the policy side. UK Government addressed integration as “integration within and between different types of transport so that each contributes its full potential and people can move easily between them” [23]. Hull [24] defined eight levels for an integrated transport policy: (8) the integration of policy measures, (7) the integration of policy sectors, (6) institutional and administrative integration, (5) the integration of environmental issues, (4) integration with social objectives, (3) integration with market needs, (2) modal integration, and (1) physical and operational integration. The different elements involved, especially the policy and institutional integration from approximately level 4 to 8, are very difficult to achieve, that makes integrated transport to be described as an ineffective practice [25].Among the 8 integration levels, the integration levels 1-2 are achievable [26]. Chowdhury et al. [7] stated that transfer is the key component of an integrated system, and proposed five main attributes of an integrated system: (1) network integration, (2) fare integration, (3) information integration, (4) the physical integration of stations, and (5) coordinated schedules. Chowdhury et al. [27] showed that physical integration is more important than information integration for current users. Ibrahim [11] proposed that the network integration is a fully multimodal integrated system that optimizes the available resources and therefore reduces the wasteful duplication of services while still maintaining an adequate coverage area. Studies have also shown that network integration can increase the usage of PT services. Matas [28] investigated the significant increase in PT use (>40%) in Madrid. The study discussed that integrated fare system and network integration had the most impact on ridership. Hidalgo and King [29] discussed the impact of the gradual implementation of an integrated PT system, which requires a combination of political leadership, technical capacity, and funding.Following the prior efforts of the integrated transportation, with the development of mobile Internet and smartphone apps, MaaS is considered as an emerging opportunity to integrate the conventional PT, new on-demand, and shared mobility services. Sochor et al. [30] put forward three major attributes of the descriptions:(1)
Offering a service with customer/user/traveller/consumer transport needs as the main focus(2)
Offering (multimodal) mobility rather than transport(3)
Offering integration of transport services, information, payment, and ticketingThe way to offer multimodal mobility services via a mobile app, not the transport modes itself, makes the combined trip more acceptable for the users. The combination of transport modes with different strengths and weaknesses is believed to be more efficient than one single mode [18]. However, even if MaaS has received a lot of attention, MaaS applications are still not a fully fledged reality. Many problems prevent its implementation and make the gap between theory and practice still high. Polydoropoulou et al. [2] elaborated on possible MaaS barriers and enablers related to infrastructures (integrated payments and open application programming interfaces), hard institutions (regulations, e.g., related to fares and security concerns), soft institutions (acceptability and trust between operators), and capabilities (need for transport investments). Some authors also highlight that bundling of mobility services is typically seen as the next step after the full integration of operations, information, and payments [6] since the infrastructure is the foundation of MaaS, providing the physical means to get from the origin to destination. The network integration via the optimization of hubs in this paper can be viewed as a hardware support for the MaaS applications.Quantitative studies to facilitate the integration scheme on modal shift are relatively sparse. Practical field tests are regularly used in the later phases of the development process, but in the earlier phases, one is generally dependent on less powerful methods, such as checklists and expert judgements [31]. Solecka [32] proposed a common variant ranking method, Electre III, to assist policy makers in deciding the most suitable version of an integrated PT system for their city. Kash and Hidalgo [33] developed a framework to include consultation with the users and incorporate their needs into a regional PT plan. Data were collected by undertaking semistructured interviews with transportation experts. Considering the rigorous amounts of data and analysis which consumes time and resources, El-Adaway et al. [34] used the social network analysis method as a complementary tool for improved transportation planning. Other studies focused on the simulation method, Santi et al. [35] proposed a method of shareability networks, which formulates the spatiotemporal sharing problem using a graph-theoretic framework. Simulations on a data set of New York City taxi trips indicate that the cumulative trip length can be reduced by 40% through ride sharing. Shen et al. [36] designed an integrated autonomous vehicles (AV)-PT system by repurposing low-demand bus routes using the shared AV service while retaining the high-demand bus routes. Through the agent-based simulations, the authors calculate the best number of AV vehicles to replace bus routes.Since the plan of the network infrastructure affects the long-term passenger distributions. The simulation methods concentrate on the dynamic and short-term evolution of the user’s behaviour, which is good at predicting the usage of on-demand or real-time traffics. However, the backbone of the multimodal network is metro, scheduled bus and car, integrating the hubs mainly serves these subnetworks. Thus, the multimodal network design problem (MNDP) which uses static traffic assignment problem is more suitable for capturing the final equilibrium flow among the entire network.The MNDP is often formulated as a bilevel mathematical program with equilibrium constraints [37–39], at the upper-level, where the optimization variables focus on [40] (1) route design [41], (2) frequency setting [42, 43], (3) timetabling [44], (4) vehicle scheduling [45], and (5) crew scheduling. MNDP is also classified according to the aimed design period into strategic, tactical, and operational stages. Route design constitutes the strategic stage, and the frequency setting is undertaken at the tactical stage, and the remaining activities are handled at the operational stage. The lower-level problem is usually a traffic assignment problem (TAP). The indexes generated by the TAP will be fed back to the upper-level problem to assess the optimization schemes.As a widely used tool to predict the flow distribution of a network, Sheffi [46] described the TAP of capturing the mode choices of users as a modal split/traffic assignment problem (MS/TAP), in which the flow of each mode was equivalent to postulating a travel demand function, where the number of trips was proportional to a negative exponential function of the travel cost and linked to route choice models based on the Wardropian user equilibrium [47]. The single-mode model was first successfully formulated and solved by Evans [48]. Then, Leblanc [49]; Safwat and Magnanti [50]; Boyce and Zhang [51]; and many others studied various [52–54] aspects of the problem, including its formulations, algorithms, cost functions, and calibration methods. Among these studies, mode choice was often made between two or three single modes (car and transit) and combined with the route choice to describe the users’ spatial distribution in the network [55, 56].As the single-mode models were being perfected, the research attention gradually moved to combined modes [57–59]. Many studies about specific combined modes, including P&R (park and ride) [60], remote P&R [61], and transit and e-hailing [62], were proposed. The mode-choice structure then changes from a single level to multiple levels and is usually formulated with a nested logit (NL) or cross-nested logit (CNL) model [63]).The upper-level part of the MNDP applies a process similar to that from a single network [64–66] to multiple networks. Only a few researchers have simultaneously tackled over two subnetwork facilities or services. Chen and Nie [62] developed two different relative spatial position designs for an integrated e-hailing/fixed-route transit system. Song et al. [67] proposed an integrated planning framework to locate park and ride (P&R) facilities and simultaneously optimize their capacities and transit service frequencies. These studies had the goal of improving the attractiveness of one combined trip by simultaneously optimizing the services of the two modes involved.As discussed above, the IMNDP of this paper focuses on integrating the physical multimodal network via the optimization of transportation hubs to achieve a balanced flow distribution within the entire network. This brings three problems in adopting it into the MNDP model structure:(1)
IMNDP involves the overall optimization of hubs in different subnetworks, where the change in one subnetwork’s hubs will influence the flow distribution among other subnetworks, so the mode choice in the MS/TAP should consider all the modes that the multimodal network can provide. In the traditional MS/TAP, the mode choice only covers the aimed design modes. However, with the development of multimodal networks, more basic travel modes are introduced in while the increase in travel distance also brings in an increase of transfers, which leads to the combined modes the multimodal network can provide are almost incalculable.(2)
Unlike the MNDP, which focuses on increasing one subnetwork’s performance, the IMNDP concentrates on the interaction between different modes. The objectives of the upper-level problem should reflect the modal shift and balanced use of the entire network.(3)
Considering the collaborative optimization of different subnetwork hubs, a large search space requires highly efficient algorithms to solve the bilevel problem.For the first problem, if we use the mode structure of the traditional MS/TAP, the large number of modes will make the mathematical formulation extremely complicated, which can be found in Liu et al. [61]. In addition, for a path-based algorithm, the listed modes require a corresponding number of paths, which further increases the calculation burden although there may be countless possible combined modes in the entire multimodal network. From the user’s perspective, they will only choose more competitive modes from a certain number of feasible alternatives. There is no need to list all modes for all OD pairs. These travel modes can be found indirectly from the path set. The modal split can be built among the modes in the path set. Using this method, the number of variables can be constrained by the number of paths, which can significantly reduce the calculation burden. Following this idea, a generalized modal split/traffic assignment problem (GMS/TAP) is proposed in this paper.The second problem is to reflect the balanced use of the entire network while cutting down the usage of the private car. Because the balanced use can be considered as no links and hubs becoming over congested, in this study, an exponential service-level function was formulated with the sum of the flow/capacity ratio of all the links and hubs. Meanwhile, in case the balanced use is achieved by increasing the car usage or other unexpected modal shift, a weighted modal shift value is formulated to guide the expected modal shift.For the third problem, the IMNDP is a mixed-integer nonlinear bilevel NP-hard problem (MIBLP) [68]. There are usually two approaches to solve an MIBLP: exact search methods (exact optimization) and nonexact methods (heuristic), including genetic algorithms (GAs), evolutionary methods, ant colony optimization, tabu search, and simulated annealing [69]. However, in a real-size network, the IMNDP experiences a combinatorial explosion in the size of the solution space as the problem size increases. The evaluation of potential candidate solutions requires the execution of a GMS/TAP, which requires a higher-order polynomial time [70]. To date, heuristic approaches have proved to be the only practical methods for solving realistic and large-sized problem instances [41].Among the heuristic algorithms, the GA has been proven to be an efficient algorithm for solving an MIBLP because of its simplicity and minimal problem restrictions [71]. Many different types of GAs have been successfully adopted for MNDP [72–75]. Since IMNDP dealt with over three sets of decision variables, combining these different variables leads to challenges regarding infeasible solution management and neighbourhood evaluation. Thus, a hybrid genetic algorithm with adaptive diversity management (HGSADC) is proposed based on the work of Vidal et al. [76]. HGSADC is a hybrid metaheuristic that combines the exploration capabilities of GAs with efficient local search-based improvement procedures and diversity management mechanisms [77, 78].The local search and diversity management in HGSADC is problem specific and deals with vehicle routing problems with time windows. To adopt HGSADC for the IMNDP, a local search with a crossover, mutation, and repair operator; diversity management with a Hamming distance-based objective; and a dynamic penalty were introduced based on the previous research [73, 79].To summarize, this study proposed a new IMNDP concept to achieve a modal shift and balanced use of the entire network by the collaborative optimization of the hubs. The main aim of this study was to provide new methodologies for the assessment and optimal design of a hub integration scheme. The major contributions of this study are as follows:(1)
This paper extended the MS/TAP to describe the mode choice for all modes in the path set, and a unified optimization problem was formulated by combining the mode choice and route choice together.(2)
A method for the optimal integration scheme was proposed, in which the modal shift and balanced use were chosen as the objectives, while the decision variables allow a collaborative optimization of hub locations, route layouts, and frequency setting for hubs with different modes. An HGSADC was proposed to address the optimal solution of the IMNDP.The abbreviations used in this paper are shown in Table1.Table 1
List of abbreviations.
Full nameAbbreviationPublic transportPTMobility as a serviceMaaSAutonomous vehiclesAVPark and rideP&RMultimodal network design problemMNDPIntegrated multimodal network design problemIMNDPTraffic assignment problemTAPModal split/traffic assignment problemMS/TAPGeneralized modal split/traffic assignment problemGMS/TAPMixed-integer nonlinear bilevel NP-hard problemMIBLPGenetic algorithmsGAHybrid genetic algorithm with adaptive diversity managementHGSADCMethod of successive averagesMSAMultinomial logitMLNested logitNLCross-nested logitCNLUser equilibriumUEVariational inequalityVINonlinear optimization problemNLPBureau of public roadBPR


## 3. Materials and Methods

In this section, the bilevel IMNDP is formulated for both the lower-level and upper-level problem. The network representation is built in Section3.1 at first. Then the formulation of lower-level GMS/TAP is divided into two parts: the mode-choice structure and the user equilibrium (UE) principle. The method of building the mode-choice structure is shown in Section 3.2. After combining the mode-choice structure with the user equilibrium principle, a unified GMS/TAP is formulated in Section 3.3 along with its solution method in Section 3.4. The upper-level model is introduced in Section 3.5, and the HGSADC and its operators are provided in Section 3.6.


### 3.1. Network Representation

Five basic modes were considered in this study, i.e., bicycle, car hailing, car, bus, and metro. The car, scheduled bus, and metro are the backbone of the urban transportation, complemented by bicycle and the online car hailing services such as Uber and Didi. The mode set is denoted byM=bic,ch,car,bus,mtr. Denote the subnetworks as GmNm,Amsub,Hm, where each subnetwork corresponds to a basic travel mode, m, where m∈M. Nm, Amsub, and Hm are the node, subnetwork link, and hub set, respectively. These subnetworks are combined with three link subsets. The set of auxiliary links, Aaux, connects the subnetwork nodes and their hubs within an acceptable walking distance. The sets of embarking and alighting links (Aemb and Aali, respectively) connect the OD nodes with the subnetwork nodes and represent the embarking and alighting activities. The set of transfer links, Atra, connects different hubs and represents the transfer activities from one mode to another.The cost of the subnetwork link is the in-vehicle time, which is a link-flow-dependent function. The alighting/embarking link and auxiliary link comprise the walking time, waiting time, and time consumed in the hub, which is related to the length of the link, route frequency, and hub flow. The impedance of the transfer link adds a pure penalty based on the above three times, which represents the user’s natural resistance to transfer. The impedance functions of the four links are continuously differentiable and monotonically increasing functions.


### 3.2. Mode-Choice Structure

To build a GMS/TAP, a mode-choice structure is first required to represent the travel modes. A modified NL model is used in this paper, see Figure1. Suppose the maximum transfer number is K. In level −1, all modes are classified according to their number of transfers and divided into K + 1 columns (from 0 to K transfers). Then, for each column k, a sub-NL model is built by modifying the NL structure for modes with k transfers. Then, from level 0 to level k, each basic mode of the k transfers combined mode is listed level by level. In the mode-choice structure, nests in column k, level k are called mode nests (nests in blue boxes), which represent complete combined modes, while the nests in level 0 to k − 1 are called intermediate nests (nests in yellow boxes), which are the first n modes of one complete combined mode.Figure 1
Mode-choice structure based on the nested logit model.To build a NL choice structure from the path set. A given path set of OD pairw is denoted as Pw, where Pw can be divided into K + 1 subsets Pkw according to the number of transfers, k=0∼K. For the sth path pk,sw in Pkw, s=0∼Skw. The mode of pk,sw can be denoted with a mode chain mck,swm0,m1,…,mn,…,mk, where mn is the nth basic mode, n=0∼k. For convenience, a complete mode chain, mck,swm0,m1,…,mn,…,mk is denoted as mck,swk, and an intermediate mode chain with the first n modes, mck,swm0,m1,...,mn is denoted as mck,swn.For the NL mode-choice structure, denote theqth nest in column k and level n as ek,qwm0,m1,…,mn, which can be simplified to ek,qwn. In this study, if the mode sequence in mck,swn is the same as the mode sequence in ek,qwn, then mck,swn=ek,qwn. If the first n modes in ek,qwn+1 are the same as ek,qwn, then ek,qwn+1∈ek,qwn, meaning ek,qwn+1 is the lower-level nest of ek,qwn.Following this corresponding relationship from mode chain to nest, for a given path set, the NL choice-structure can be built by matching the mode chains and the nests one by one.


### 3.3. Mathematical Formulation of the GMS/TAP

#### 3.3.1. The Three Conditions

After building the mode-choice structure, the NL rules should be formulated among the flow of each nest. Meanwhile, to combine the mode choice and route choice together, this paper assumes that the path flow follows the user equilibrium principle. Thus, the following 3 conditions should be satisfied at the same time for a unified formulation.C1:Choice of Transfer Number. After the proportions are satisfied, no user will change the transfer number. Suppose the demand of OD pair w, gw is fixed. In the mode-choice structure, the flow of nests in column k, level −1 is denoted as gk,qw−1. Then, C1 is formulated as ∀w∈W, k=0∼K, and gk,qw−1, and gw should satisfy the following constraints:(1)gk,qw−1=exp−α−1+β−1Uk,qw−1∑k=0Kexp−α−1+β−1Uk,qw−1gw,whereUk,qw−1 is the equilibrium utility of the nest in column k, level −1, α−1, and β−1 are the logit model parameters associated with level −1.C2: Choice of Mode. After the proportions are satisfied, no user will change the mode sequence. Denote the flow of each intermediate nest as gk,qwn and its lower-level nest flow as gk,qwn+1. C2 is represented as ∀w∈W, k = 0∼K and n = 0∼k − 1, and the following constraints should be obeyed:(2)gk,qwn+1=exp−αkn+1+βkn+1Uk,qwn+1∑q=1Qkwexp−αkn+1+βkn+1Uk,qw¯n+1jk,qwn+1gk,qwn,(3)jk,qwn+1=1,ifek,qwn+1∈ek,qwn,0,otherwise,whereUk,qwn+1 is the equilibrium utility of nest ek,qwn+1, Uk,qw¯n+1 is the equilibrium utility set of all nests in level n + 1, and jk,qwn+1 is a binary variable for determining whether ek,qwn+1∈ek,qwn.C3: Choice of Path. Paths in each mode nest should follow the UE principle. Denote the cost of path pk,sw as ck,sw. C3 is represented as ∀w∈W, k=0∼K, s=1∼Skw, and mck,swk=ek,qwk:(4)ck,sw−Uk,qwk=0,iffk,sw∗>0,≥0,iffk,sw∗=0.Equation (4) states that only when path cost ck,sw equals its corresponding mode nest’s equilibrium utility Uk,qwk, then can path pk,sw carry flow; otherwise, the path flow is 0. It should be noted that the condition mck,swk=ek,qwk is to find the corresponding mode nest ek,qwk of the complete mode chain mck,swk, and only paths in the same mode nest should obey the UE principle.However, the variable of C1 and C2 is the nest flow, while for C3, it is the path flow. To keep the variables consistent in the three conditions so that a unified formulation can be built to combine the mode choice and route choice together, a simple transformation is built.Denoteg¯k,sw−1 as the flow of path pk,sw’s corresponding nest in level −1 and g¯k,swn as the flow of mck,swn’s corresponding nest in level n, which means g¯k,sw−1=gk,qw−1 and g¯k,swn=gk,qwn. Then, C1 can be reformulated as ∀w∈W, k=0∼K, s=1∼Skw, and mck,sw−1=ek,qw−1:(5)gk,sw¯−1=exp−α−1+β−1Ukw∑k=0Kexp−α−1+β−1Ukwgw.While C2 can be reformulated as∀w∈W, k = 0∼K, s=1∼Skw, n = 0∼k – 1, and mck,swn=ek,qwn:(6)g¯k,swn+1=exp−αkn+1+βkn+1Uk,qwn+1∑q=1Qkwexp−αkn+1+βkn+1Uk,qw¯n+1jk,qwn+1g¯k,swn.Equations (5) and (6) only replaces gk,qw−1 and gk,qwn in equations (1) and (2) with g¯k,sw−1 and g¯k,swn, and the other variables are the same.


#### 3.3.2. Unified Optimization Problem

Denote the equilibrium flow vector of all OD pairs asF∗. The following variational inequality (VI) formulation was built to unify the above three conditions. The proof is presented in Appendix A.Proposition 1.
The optimal solutionF∗ of the proposed VI problem gives the NL mode choice solution and UE route choice solution together only if there exists the set of values λw∗, for∀w∈W,k=0∼K,n=0∼k,s=1∼Skw, mck,sw−1=ek,qw−1, and mck,swn=ek,qwn;(7)ck,swF∗+Λk,swF∗−λw∗F∗T⋅Fw−Fw∗≥0,where(8)Λk,swF∗=lng¯k,sw−1+α−1β−1,whenk=0,lng¯k,sw−1+α−1β−1−∑n=0k−lng¯k,swn+lng¯k,swn−1−αknβkn,whenk=1,2,…,K,λw∗F∗=Uw+lngwβ−1.
The solution space of the VI problem is(9)∑k=0K∑s=1Sfk,sw∗=gw,(10)F∗≥0.
Equation (9) ensures that the total path flow of OD w equals demand gw, and equation (10) indicates that the flow is nonnegative.
The VI problem can be reformulated as a single-level NLP according to Aghassi et al. [80]. Since the solution space (equations (9) and (10)) of the VI problem can be reformulated as a nonempty polyhedron, D=F∈ℝn|B⋅F=G,F≥0, where F∈ℝn indicates that flow vector F is an n-dimensional real number vector, G=g1,…,gw,…,gWT, B is the incident matrix of equation (9), denoting ck,swF∗+Λk,swF∗−λw∗F∗ as HF∗. The VI problem in equation (7) can be transformed into the following NLP.Proposition 2.
Suppose thatD is a nonempty polyhedron, D=F∈ℝn|B⋅F=G,F≥0. Then, F∗ solves the VI problem in equation (7) if the following optimization problem has an optimal value of 0 and ∃ξ∗∈ℝW, such that F∗,ξ∗ is an optimal solution:(11)minHFTF−GTξ,s.t.B⋅F=G,F≥0,BTξ≤HF.


### 3.4. Hybrid MSA Algorithm

The NLP can be solved using the currently available software such as CPLEX, Gurobi, and MATLAB fmincon solver. Meanwhile, the standard formulation provides an excellent coverage criterion, where the minimum objective function should be zero. However, the number of variables in the two vectors,F and ξ, which are related to the number of paths and number of OD pairs, respectively, can still be very large in real-size network. A highly efficient algorithm is necessary to apply the GMS/TAP to the IMNDP. Thus, a method of successive averages (MSA)-based method is proposed, which is commonly used to solve traffic assignment problems [81].


### 3.5. Formulation of the Upper-Level Problem

The notations used in the upper-level problem are shown in Table2.Table 2
Notations of the upper-level problem.
SetCNSet of candidate nodesASet of linksHSet of hubsMSet of modesParametersmsl,mThe weight when flow of modem with direct OD distance l changes one unitLwThe direct distance between OD pairwcaphm, capamThe capacity of the hub or link with modemψxA binary function which equals 1 ifx = 1, and 0, otherwisedcmsub, ccmsubThe demolition cost or construction cost of modem hubBPThe purchase cost of one busϑxA binary function which equals 1 ifx > 0, and 0, otherwiseocThe operational cost for one bus traveling one kilometerRmaxThe maximum number of bus routesLminstop, LmaxstopThe minimum and maximum bus stop distanceLmin, LmaxThe minimum and maximum bus route lengthfremin, fremaxThe minimum and maximum bus frequencyFSmin, FSmaxThe minimum and maximum of bus fleet sizeHNmin, HNmaxThe minimum and maximum number of hubs with modemB1∼B4The weight of four objectivesVariablesmfl,m, mfl,moriAmong OD pairs with direct distance betweenl − 1 and l km, the total flow of mode m in the optimization scheme or in the original networkslhm, slamThe service level value of the hub or link with modemXojr1 if router starts at node j, and 0 otherwiseXior1 if router ends at node i, and 0 otherwiseXijr1 if router (r=1toRmax) passes through node j immediately after node i, and 0 otherwiseXihmter, Xihmint1 if terminal or interchange is allocated to candidate nodei, and 0, otherwiseXihm1 if modem hub is allocated in node i, and 0 otherwiselijrThe length between nodei and j of route rlrThe length of routerfrerThe frequency of routerfsrThe fleet size of routerTTThe total travel time of all usersTCThe total construction cost of the schemeMSThe modal shift value of the schemeBUThe balanced use value of the schemeIMNDP aims to realize a modal shift and obtain a balanced use among the entire multimodal network under budget limits. This results in three objectives: construction cost, modal shift, and balanced use. Because the balanced use can be considered as no links and hubs becoming over congested, an exponential service-level function was formulated with the sum of the flow/capacity ratio of all the links and hubs. Meanwhile, in case the balanced use is achieved by increasing the car usage or other unexpected modal shift, a weighted modal shift value is formulated. As a widely used objective, the total time of all users was also adopted to protect the users’ interests. Then, the objective function of the problem is(12)minZ=B1TT+B2TC+B3MS+B4BU.Since it is beyond the work of this paper to cover all hubs of different modes, for simplification, only the five modes in Section3.1 are considered. The optimization of hubs is classified into two types. One is the location-based optimization, which are parking lot and bicycle station. The optimization only needs to determine their locations. The other is the route-based optimization, which is the bus stop. The change of their hub locations will also influence the route layouts and frequencies. The integration scheme is to optimize these variables simultaneously to build a well-connected multimodal network. The objective functions are calculated with equations (13)–(16) using these results:(13)TT=∑w=1W∑k=0K∑s=1Skwfk,swck,sw,(14)TC=∑i∈CN∑m∈MdcmψXihmori−Xihm+∑i∈CN∑m∈MccmψXihm−Xihmori+BP∑r=1Rmaxϑfsr−fsrori+oc∑r=1Rmaxlrfsr,(15)MS=∑m=1M∑l=1Lmfl,m−mfl,morimsl,m,(16)MS=∑hm∈Hexpfhmcaphm+θ+∑am∈Aexpfamcapam+θ.Equation (13) is the total travel time. Equation (14) is the total construction cost, which consists of the demolition and construction costs of the hubs, the bus purchase cost, and the bus operational cost. Equation (15) is the modal shift value. To analyse the modal shift of OD pairs with different lengths, mfl,m is denoted as the total flow of mode m among all OD pairs with distance l. The predefined weight parameter msl,m is used to value the 1 unit of flow change compared to the original network if msl,m is set as a minus value, which means the modal shift of this mode in this distance is encouraged to increase, and vice versa. Equation (16) is the service-level value, which is an exponential service-level function with the sum of the flow/capacity ratio of all the links and hubs. To ensure that highly loaded hubs or links will have higher service level value penalty, a parameter θ is added in equation (16).For location-based optimization, the constraints include those as follows:(17)∑i∈CN,m∈car,bicXihm=1,(18)HNmin≤∑i∈Hm,m∈car,bicXihm≤HNmax.Equation (17) ensures that for each candidate node, only one hub with mode m can be placed. Equation (18) is the hub number constraint for hubs with mode m.For route-based optimization, where the decision variables are bus layouts and frequencies, referring to [8, 73], the constraints include those as follows:(19)∑j∈HbusterXojr=1forr=1toRmax,(20)∑i∈HbusterXior=1forr=1toRmax,(21)∑i∈Hbusint,i≠jXijr−∑i∈Hbusint,i≠jXjir=0j∈Hbusint,forr=1toRmax,(22)∑i∈Hbusint,i≠jXijr≤1j∈Hbusint,forr=1toRmax,(23)∑j∈Hbusint,i≠jXijr≤1forj∈Hbusint,r=1toRmax,(24)Xjjr=0forj∈Hbusint,r=1toRmax,(25)Lminstop≤Xijrlijr≤Lmaxstopfori∈Hbusint,i≠j,j∈Hbusint,forr=1toRmax,(26)Lmin≤lr≤Lmaxforr=1toRmax,(27)fremin≤frer≤fremaxforr=1toRmax,(28)FSmin≤∑r=1Rmaxfsr≤FSmax,(29)frer=fsrlr2vbus.Among the constraints, equations (19) and (20) ensure that all bus routes start and end at a bus terminal. Equation (21) ensures that each bus interchange has one preceding node and one subsequent node. Equations (22)–(24) ensure that each bus node can be visited by one particular bus route only once at most (bus node repeat constraint). Equations (25)–(28) are the bus stop distance, bus route length, bus frequency, and bus fleet size constraints, respectively. The bus frequency of one route is calculated using the assigned number of buses in equation (29).


### 3.6. Solution Method of the Upper-Level Problem

The upper-level problem is a MIBLP with 4 sets of variables. To solve the problem, the HGSADC framework is introduced as follows.The local search and diversity management in HGSADC is problem specific. Thus, the local search with a crossover, mutation, and repair operator; diversity management with a Hamming distance-based objective; and a dynamic penalty were introduced one by one in the following sections.


#### 3.6.1. Initialization Operator

For route-based optimization, the initialization operator determines the initialization locations of the bus stops and their sequences in each bus route. The initialization scheme should satisfy the stop distance constraint, route length constraint, and bus node repeat constraint.The bus route frequency is initialized by randomly assigning buses to each route according to equation (29). If the frequency violates frequency constraint and the number of buses violates fleet size constraint, the process will repeat until the constraints are satisfied. It should be noted that lines 4 ∼ 12 is called the bus sequence generation function in the rest of paper. This function can be used to generate a bus sequence with defined starting and ending nodes.For location-based optimization, the parking lot and bicycle station initialization should determine the number and position of each hub. The initialization operator randomly adds or deletes several hubs and ensures that the total number of hubs satisfies the hub number constraint.


#### 3.6.2. Crossover Operator

The crossover of the bus route consists of two levels: combinations of different routes and different sequences of bus stops in one route. To increase the exploration capability and ensure the satisfaction of the constraints for the generated routes after crossover, referring to Ngamchai and Lovell [79], two crossover methods for route-based optimization, i.e., line crossover and stop crossover, are used. Line crossover exchanges two different routes in two solutions, while the stop crossover operator is triggered when two routes in one solution have segments that can be spliced.For location-based crossover, two hubs in two solutions are randomly selected; if they are different and not repeated in the other solution, then exchange.


#### 3.6.3. Mutation Operator

The mutation process generates different genetic materials to facilitate exploration in the search space. According to Szeto and Wu [73], four types of mutation operators for route-based optimization are proposed: insert, delete, swap, and transfer.The insert mutation is used to insert an interchange node in a randomly selected route. The deletion mutation is similar. The swap mutation is to exchange the nodes of two routes if the stop distance constraint can be satisfied after the swap, while transfer mutation is to transfer one node to the other route, and similarly, the stop distance constraint should be obeyed after the transfer.The location-based mutation only contains the insert and delete mutations. The insert mutation inserts a random hub that is not in the solution, whereas the delete mutation deletes one random hub in the solution.


#### 3.6.4. Repair Operator

After the crossover and mutation processes, the generated children may violate the constraints, and a repair operator is needed to fix these infeasible solutions. The constraints for route-based optimization are route length, stop distance, node repeat, and frequency, while the constraint for location-based optimization is the number of hubs. With four constraints, there is a high probability that the bus route will violate other constraints after one constraint is repaired. A new four-stage repair strategy was proposed in this study. The solution set will undergo four repair operators in the sequence of the route length operator, stop distance operator, node repeat operator, and route frequency operator.(1)
Route Length Repair Operator. See Figure 2, when route r1 is too short, the algorithm selects two nodes, n1a and n1b, with a sufficient distance and inserts an intermediate node n1c between them. Then, two stop sequences from n1a to n1c and from n1c to n1b are generated with the bus stop sequence generation function in Algorithm1, and the original stop sequence between n1a and n1b is replaced. The situation when the route is too long is similar, and two nodes with a sufficient distance are chosen and replaced with the generated stop sequence using the bus stop sequence generation function.(2)
Stop Distance Repair Operator. A stop sequence that violates the stop distance constraint may contain only two nodes or consist of several consecutive nodes in which all the adjacent nodes violate the stop distance constraint. When the stop distance is too long, and the stop sequence has only two nodes, see Figure 3, the algorithm finds an intermediate node to insert between the two nodes. When the stop sequence has over two nodes, a new stop sequence will be generated to replace it. The situation when stop distance is too short is similar.(3)
Route Repeat Repair Operator. For repeated nodes, one node is randomly selected to be unchanging. Then, for each repeated node, n1a, a node is selected from CN that satisfies the stop distance constraint with both n1a − 1 and n1a + 1 to replace n1a.(4)
Route Frequency Repair Operator. The frequency repair operator adds or deletes one bus until the bus frequency constraint is satisfied.Algorithm 1: Bus route initialization.
(1)
for route r from 1 to Rmax(2)
Whether router needs to construct a new terminal with probability Pters and Ptere(3)
Determine the starting and ending node: if new terminal should be built, select a node fromCN as the starting node nrs or ending node nre of route r; if no, then select a node from TN as the starting node or ending node(4)
while the tth node of roue r, nrt is not the nre, and the route satisfies the route length constraint, bus node repeat constraint, and time < Tmaxdo(5)
Find the node setNrstopt with distance to nrt satisfying the stop distance constraint(6)
ifNrstopt is empty or all nodes have been visited in Nrstopt(7)
nrt=0, t = t-1(8)
elseifNrstopt is not empty(9)
Randomly select a nodenrstopt from Nrstopt that has not been visited(10)
ifnre−nrtempt<nre−nrtnrt+1=nrtemptelseifnre−nrtempt≥nre−nrtMarknrstopt has been visited and go back to line 6end(11)
if route r is feasible thenInsert router into initialization solution(12)
elseif route r is infeasible thenGo back to line 2.end


#### 3.6.5. Diversification Mechanism and Infeasible Solution Management

To ensure the diversity of the solutions, the Hamming distance is introduced to represent the differences between one solution and the best solution, which is also called the diversity contribution functionΔS. It is defined as the average distance from S to its closest neighbours in the subpopulation. ΔS consists of three parts, which are the bus route, parking lot, and bicycle station Hamming distance. According to Szeto and Wu [73], the bus Hamming distance is the number of different consecutive node pairs for each pair of corresponding routes compared to the best solution. For the location-based one, it is the number of different hubs compared to the best solution, and ΔS is the sum of the above three distances.Figure 2
Situation 1 in the length repair operator.Figure 3
Situation 1 of the stop distance repair operator.A dynamic penalty is then introduced to handle the infeasible solutions. Suppose the number of feasible and infeasible solutions areNinf and Nfea, respectively. A basic penalty value, ϖ, is adjusted for every μ iterations, and when the number of Ninf increases by 1, ϖ decreases by λ%, and vice versa. (see Algorithm 2)Algorithm 2: HGSADC.
(1)
Initialization(2)
while number of iterations without improvements < ItNIdo(3)
while maximum solution size is not reached doselect parent solutionsS1 and S2, create offspring O using route-based and location-based crossover operator. Insert O into solution set (crossover)end while(4)
select solutionS1 and mutate with route-based and location-based mutation operator (mutation)(5)
for each solutionS1, if infeasible then(6)
RepairS1 with probability Prep (repair)if repairable theninsert into feasible subpopulation.elseinsert into infeasible subpopulationend if(7)
elseifS1 feasible theninsert into feasible subpopulationend if(8)
calculate objective values with GMS/TAP(9)
calculate the biased fitness valueSelect survivors(selection)(10)
if best solution not improved for Itdiv iterations, thenDiversify population(11)
Adjust penalty parameters for infeasibilityend whileReturn best feasible solutionThen, the survivors of each iteration are selected based on the rank of objective functionfitS and the rank of Hamming distance hdrS, which is called the biased fitness in this paper. For infeasible solutions, the objective function should multiply dynamic penalty values, ϖ. Referring to Vidal et al. [76], the biased fitness, BFS, is formulated:(30)BFS=fitS+1−NeliteNmaxhdrS,where Nelite is a parameter ensuring elitism properties during survivor selection.After building the whole bilevel model, the extracts of the algorithms can be found in Appendix B.


## 3.1. Network Representation

Five basic modes were considered in this study, i.e., bicycle, car hailing, car, bus, and metro. The car, scheduled bus, and metro are the backbone of the urban transportation, complemented by bicycle and the online car hailing services such as Uber and Didi. The mode set is denoted byM=bic,ch,car,bus,mtr. Denote the subnetworks as GmNm,Amsub,Hm, where each subnetwork corresponds to a basic travel mode, m, where m∈M. Nm, Amsub, and Hm are the node, subnetwork link, and hub set, respectively. These subnetworks are combined with three link subsets. The set of auxiliary links, Aaux, connects the subnetwork nodes and their hubs within an acceptable walking distance. The sets of embarking and alighting links (Aemb and Aali, respectively) connect the OD nodes with the subnetwork nodes and represent the embarking and alighting activities. The set of transfer links, Atra, connects different hubs and represents the transfer activities from one mode to another.The cost of the subnetwork link is the in-vehicle time, which is a link-flow-dependent function. The alighting/embarking link and auxiliary link comprise the walking time, waiting time, and time consumed in the hub, which is related to the length of the link, route frequency, and hub flow. The impedance of the transfer link adds a pure penalty based on the above three times, which represents the user’s natural resistance to transfer. The impedance functions of the four links are continuously differentiable and monotonically increasing functions.


## 3.2. Mode-Choice Structure

To build a GMS/TAP, a mode-choice structure is first required to represent the travel modes. A modified NL model is used in this paper, see Figure1. Suppose the maximum transfer number is K. In level −1, all modes are classified according to their number of transfers and divided into K + 1 columns (from 0 to K transfers). Then, for each column k, a sub-NL model is built by modifying the NL structure for modes with k transfers. Then, from level 0 to level k, each basic mode of the k transfers combined mode is listed level by level. In the mode-choice structure, nests in column k, level k are called mode nests (nests in blue boxes), which represent complete combined modes, while the nests in level 0 to k − 1 are called intermediate nests (nests in yellow boxes), which are the first n modes of one complete combined mode.Figure 1
Mode-choice structure based on the nested logit model.To build a NL choice structure from the path set. A given path set of OD pairw is denoted as Pw, where Pw can be divided into K + 1 subsets Pkw according to the number of transfers, k=0∼K. For the sth path pk,sw in Pkw, s=0∼Skw. The mode of pk,sw can be denoted with a mode chain mck,swm0,m1,…,mn,…,mk, where mn is the nth basic mode, n=0∼k. For convenience, a complete mode chain, mck,swm0,m1,…,mn,…,mk is denoted as mck,swk, and an intermediate mode chain with the first n modes, mck,swm0,m1,...,mn is denoted as mck,swn.For the NL mode-choice structure, denote theqth nest in column k and level n as ek,qwm0,m1,…,mn, which can be simplified to ek,qwn. In this study, if the mode sequence in mck,swn is the same as the mode sequence in ek,qwn, then mck,swn=ek,qwn. If the first n modes in ek,qwn+1 are the same as ek,qwn, then ek,qwn+1∈ek,qwn, meaning ek,qwn+1 is the lower-level nest of ek,qwn.Following this corresponding relationship from mode chain to nest, for a given path set, the NL choice-structure can be built by matching the mode chains and the nests one by one.


## 3.3. Mathematical Formulation of the GMS/TAP

### 3.3.1. The Three Conditions

After building the mode-choice structure, the NL rules should be formulated among the flow of each nest. Meanwhile, to combine the mode choice and route choice together, this paper assumes that the path flow follows the user equilibrium principle. Thus, the following 3 conditions should be satisfied at the same time for a unified formulation.C1:Choice of Transfer Number. After the proportions are satisfied, no user will change the transfer number. Suppose the demand of OD pair w, gw is fixed. In the mode-choice structure, the flow of nests in column k, level −1 is denoted as gk,qw−1. Then, C1 is formulated as ∀w∈W, k=0∼K, and gk,qw−1, and gw should satisfy the following constraints:(1)gk,qw−1=exp−α−1+β−1Uk,qw−1∑k=0Kexp−α−1+β−1Uk,qw−1gw,whereUk,qw−1 is the equilibrium utility of the nest in column k, level −1, α−1, and β−1 are the logit model parameters associated with level −1.C2: Choice of Mode. After the proportions are satisfied, no user will change the mode sequence. Denote the flow of each intermediate nest as gk,qwn and its lower-level nest flow as gk,qwn+1. C2 is represented as ∀w∈W, k = 0∼K and n = 0∼k − 1, and the following constraints should be obeyed:(2)gk,qwn+1=exp−αkn+1+βkn+1Uk,qwn+1∑q=1Qkwexp−αkn+1+βkn+1Uk,qw¯n+1jk,qwn+1gk,qwn,(3)jk,qwn+1=1,ifek,qwn+1∈ek,qwn,0,otherwise,whereUk,qwn+1 is the equilibrium utility of nest ek,qwn+1, Uk,qw¯n+1 is the equilibrium utility set of all nests in level n + 1, and jk,qwn+1 is a binary variable for determining whether ek,qwn+1∈ek,qwn.C3: Choice of Path. Paths in each mode nest should follow the UE principle. Denote the cost of path pk,sw as ck,sw. C3 is represented as ∀w∈W, k=0∼K, s=1∼Skw, and mck,swk=ek,qwk:(4)ck,sw−Uk,qwk=0,iffk,sw∗>0,≥0,iffk,sw∗=0.Equation (4) states that only when path cost ck,sw equals its corresponding mode nest’s equilibrium utility Uk,qwk, then can path pk,sw carry flow; otherwise, the path flow is 0. It should be noted that the condition mck,swk=ek,qwk is to find the corresponding mode nest ek,qwk of the complete mode chain mck,swk, and only paths in the same mode nest should obey the UE principle.However, the variable of C1 and C2 is the nest flow, while for C3, it is the path flow. To keep the variables consistent in the three conditions so that a unified formulation can be built to combine the mode choice and route choice together, a simple transformation is built.Denoteg¯k,sw−1 as the flow of path pk,sw’s corresponding nest in level −1 and g¯k,swn as the flow of mck,swn’s corresponding nest in level n, which means g¯k,sw−1=gk,qw−1 and g¯k,swn=gk,qwn. Then, C1 can be reformulated as ∀w∈W, k=0∼K, s=1∼Skw, and mck,sw−1=ek,qw−1:(5)gk,sw¯−1=exp−α−1+β−1Ukw∑k=0Kexp−α−1+β−1Ukwgw.While C2 can be reformulated as∀w∈W, k = 0∼K, s=1∼Skw, n = 0∼k – 1, and mck,swn=ek,qwn:(6)g¯k,swn+1=exp−αkn+1+βkn+1Uk,qwn+1∑q=1Qkwexp−αkn+1+βkn+1Uk,qw¯n+1jk,qwn+1g¯k,swn.Equations (5) and (6) only replaces gk,qw−1 and gk,qwn in equations (1) and (2) with g¯k,sw−1 and g¯k,swn, and the other variables are the same.


### 3.3.2. Unified Optimization Problem

Denote the equilibrium flow vector of all OD pairs asF∗. The following variational inequality (VI) formulation was built to unify the above three conditions. The proof is presented in Appendix A.Proposition 1.
The optimal solutionF∗ of the proposed VI problem gives the NL mode choice solution and UE route choice solution together only if there exists the set of values λw∗, for∀w∈W,k=0∼K,n=0∼k,s=1∼Skw, mck,sw−1=ek,qw−1, and mck,swn=ek,qwn;(7)ck,swF∗+Λk,swF∗−λw∗F∗T⋅Fw−Fw∗≥0,where(8)Λk,swF∗=lng¯k,sw−1+α−1β−1,whenk=0,lng¯k,sw−1+α−1β−1−∑n=0k−lng¯k,swn+lng¯k,swn−1−αknβkn,whenk=1,2,…,K,λw∗F∗=Uw+lngwβ−1.
The solution space of the VI problem is(9)∑k=0K∑s=1Sfk,sw∗=gw,(10)F∗≥0.
Equation (9) ensures that the total path flow of OD w equals demand gw, and equation (10) indicates that the flow is nonnegative.
The VI problem can be reformulated as a single-level NLP according to Aghassi et al. [80]. Since the solution space (equations (9) and (10)) of the VI problem can be reformulated as a nonempty polyhedron, D=F∈ℝn|B⋅F=G,F≥0, where F∈ℝn indicates that flow vector F is an n-dimensional real number vector, G=g1,…,gw,…,gWT, B is the incident matrix of equation (9), denoting ck,swF∗+Λk,swF∗−λw∗F∗ as HF∗. The VI problem in equation (7) can be transformed into the following NLP.Proposition 2.
Suppose thatD is a nonempty polyhedron, D=F∈ℝn|B⋅F=G,F≥0. Then, F∗ solves the VI problem in equation (7) if the following optimization problem has an optimal value of 0 and ∃ξ∗∈ℝW, such that F∗,ξ∗ is an optimal solution:(11)minHFTF−GTξ,s.t.B⋅F=G,F≥0,BTξ≤HF.


## 3.3.1. The Three Conditions

After building the mode-choice structure, the NL rules should be formulated among the flow of each nest. Meanwhile, to combine the mode choice and route choice together, this paper assumes that the path flow follows the user equilibrium principle. Thus, the following 3 conditions should be satisfied at the same time for a unified formulation.C1:Choice of Transfer Number. After the proportions are satisfied, no user will change the transfer number. Suppose the demand of OD pair w, gw is fixed. In the mode-choice structure, the flow of nests in column k, level −1 is denoted as gk,qw−1. Then, C1 is formulated as ∀w∈W, k=0∼K, and gk,qw−1, and gw should satisfy the following constraints:(1)gk,qw−1=exp−α−1+β−1Uk,qw−1∑k=0Kexp−α−1+β−1Uk,qw−1gw,whereUk,qw−1 is the equilibrium utility of the nest in column k, level −1, α−1, and β−1 are the logit model parameters associated with level −1.C2: Choice of Mode. After the proportions are satisfied, no user will change the mode sequence. Denote the flow of each intermediate nest as gk,qwn and its lower-level nest flow as gk,qwn+1. C2 is represented as ∀w∈W, k = 0∼K and n = 0∼k − 1, and the following constraints should be obeyed:(2)gk,qwn+1=exp−αkn+1+βkn+1Uk,qwn+1∑q=1Qkwexp−αkn+1+βkn+1Uk,qw¯n+1jk,qwn+1gk,qwn,(3)jk,qwn+1=1,ifek,qwn+1∈ek,qwn,0,otherwise,whereUk,qwn+1 is the equilibrium utility of nest ek,qwn+1, Uk,qw¯n+1 is the equilibrium utility set of all nests in level n + 1, and jk,qwn+1 is a binary variable for determining whether ek,qwn+1∈ek,qwn.C3: Choice of Path. Paths in each mode nest should follow the UE principle. Denote the cost of path pk,sw as ck,sw. C3 is represented as ∀w∈W, k=0∼K, s=1∼Skw, and mck,swk=ek,qwk:(4)ck,sw−Uk,qwk=0,iffk,sw∗>0,≥0,iffk,sw∗=0.Equation (4) states that only when path cost ck,sw equals its corresponding mode nest’s equilibrium utility Uk,qwk, then can path pk,sw carry flow; otherwise, the path flow is 0. It should be noted that the condition mck,swk=ek,qwk is to find the corresponding mode nest ek,qwk of the complete mode chain mck,swk, and only paths in the same mode nest should obey the UE principle.However, the variable of C1 and C2 is the nest flow, while for C3, it is the path flow. To keep the variables consistent in the three conditions so that a unified formulation can be built to combine the mode choice and route choice together, a simple transformation is built.Denoteg¯k,sw−1 as the flow of path pk,sw’s corresponding nest in level −1 and g¯k,swn as the flow of mck,swn’s corresponding nest in level n, which means g¯k,sw−1=gk,qw−1 and g¯k,swn=gk,qwn. Then, C1 can be reformulated as ∀w∈W, k=0∼K, s=1∼Skw, and mck,sw−1=ek,qw−1:(5)gk,sw¯−1=exp−α−1+β−1Ukw∑k=0Kexp−α−1+β−1Ukwgw.While C2 can be reformulated as∀w∈W, k = 0∼K, s=1∼Skw, n = 0∼k – 1, and mck,swn=ek,qwn:(6)g¯k,swn+1=exp−αkn+1+βkn+1Uk,qwn+1∑q=1Qkwexp−αkn+1+βkn+1Uk,qw¯n+1jk,qwn+1g¯k,swn.Equations (5) and (6) only replaces gk,qw−1 and gk,qwn in equations (1) and (2) with g¯k,sw−1 and g¯k,swn, and the other variables are the same.


## 3.3.2. Unified Optimization Problem

Denote the equilibrium flow vector of all OD pairs asF∗. The following variational inequality (VI) formulation was built to unify the above three conditions. The proof is presented in Appendix A.Proposition 1.
The optimal solutionF∗ of the proposed VI problem gives the NL mode choice solution and UE route choice solution together only if there exists the set of values λw∗, for∀w∈W,k=0∼K,n=0∼k,s=1∼Skw, mck,sw−1=ek,qw−1, and mck,swn=ek,qwn;(7)ck,swF∗+Λk,swF∗−λw∗F∗T⋅Fw−Fw∗≥0,where(8)Λk,swF∗=lng¯k,sw−1+α−1β−1,whenk=0,lng¯k,sw−1+α−1β−1−∑n=0k−lng¯k,swn+lng¯k,swn−1−αknβkn,whenk=1,2,…,K,λw∗F∗=Uw+lngwβ−1.
The solution space of the VI problem is(9)∑k=0K∑s=1Sfk,sw∗=gw,(10)F∗≥0.
Equation (9) ensures that the total path flow of OD w equals demand gw, and equation (10) indicates that the flow is nonnegative.
The VI problem can be reformulated as a single-level NLP according to Aghassi et al. [80]. Since the solution space (equations (9) and (10)) of the VI problem can be reformulated as a nonempty polyhedron, D=F∈ℝn|B⋅F=G,F≥0, where F∈ℝn indicates that flow vector F is an n-dimensional real number vector, G=g1,…,gw,…,gWT, B is the incident matrix of equation (9), denoting ck,swF∗+Λk,swF∗−λw∗F∗ as HF∗. The VI problem in equation (7) can be transformed into the following NLP.Proposition 2.
Suppose thatD is a nonempty polyhedron, D=F∈ℝn|B⋅F=G,F≥0. Then, F∗ solves the VI problem in equation (7) if the following optimization problem has an optimal value of 0 and ∃ξ∗∈ℝW, such that F∗,ξ∗ is an optimal solution:(11)minHFTF−GTξ,s.t.B⋅F=G,F≥0,BTξ≤HF.


## 3.4. Hybrid MSA Algorithm

The NLP can be solved using the currently available software such as CPLEX, Gurobi, and MATLAB fmincon solver. Meanwhile, the standard formulation provides an excellent coverage criterion, where the minimum objective function should be zero. However, the number of variables in the two vectors,F and ξ, which are related to the number of paths and number of OD pairs, respectively, can still be very large in real-size network. A highly efficient algorithm is necessary to apply the GMS/TAP to the IMNDP. Thus, a method of successive averages (MSA)-based method is proposed, which is commonly used to solve traffic assignment problems [81].


## 3.5. Formulation of the Upper-Level Problem

The notations used in the upper-level problem are shown in Table2.Table 2
Notations of the upper-level problem.
SetCNSet of candidate nodesASet of linksHSet of hubsMSet of modesParametersmsl,mThe weight when flow of modem with direct OD distance l changes one unitLwThe direct distance between OD pairwcaphm, capamThe capacity of the hub or link with modemψxA binary function which equals 1 ifx = 1, and 0, otherwisedcmsub, ccmsubThe demolition cost or construction cost of modem hubBPThe purchase cost of one busϑxA binary function which equals 1 ifx > 0, and 0, otherwiseocThe operational cost for one bus traveling one kilometerRmaxThe maximum number of bus routesLminstop, LmaxstopThe minimum and maximum bus stop distanceLmin, LmaxThe minimum and maximum bus route lengthfremin, fremaxThe minimum and maximum bus frequencyFSmin, FSmaxThe minimum and maximum of bus fleet sizeHNmin, HNmaxThe minimum and maximum number of hubs with modemB1∼B4The weight of four objectivesVariablesmfl,m, mfl,moriAmong OD pairs with direct distance betweenl − 1 and l km, the total flow of mode m in the optimization scheme or in the original networkslhm, slamThe service level value of the hub or link with modemXojr1 if router starts at node j, and 0 otherwiseXior1 if router ends at node i, and 0 otherwiseXijr1 if router (r=1toRmax) passes through node j immediately after node i, and 0 otherwiseXihmter, Xihmint1 if terminal or interchange is allocated to candidate nodei, and 0, otherwiseXihm1 if modem hub is allocated in node i, and 0 otherwiselijrThe length between nodei and j of route rlrThe length of routerfrerThe frequency of routerfsrThe fleet size of routerTTThe total travel time of all usersTCThe total construction cost of the schemeMSThe modal shift value of the schemeBUThe balanced use value of the schemeIMNDP aims to realize a modal shift and obtain a balanced use among the entire multimodal network under budget limits. This results in three objectives: construction cost, modal shift, and balanced use. Because the balanced use can be considered as no links and hubs becoming over congested, an exponential service-level function was formulated with the sum of the flow/capacity ratio of all the links and hubs. Meanwhile, in case the balanced use is achieved by increasing the car usage or other unexpected modal shift, a weighted modal shift value is formulated. As a widely used objective, the total time of all users was also adopted to protect the users’ interests. Then, the objective function of the problem is(12)minZ=B1TT+B2TC+B3MS+B4BU.Since it is beyond the work of this paper to cover all hubs of different modes, for simplification, only the five modes in Section3.1 are considered. The optimization of hubs is classified into two types. One is the location-based optimization, which are parking lot and bicycle station. The optimization only needs to determine their locations. The other is the route-based optimization, which is the bus stop. The change of their hub locations will also influence the route layouts and frequencies. The integration scheme is to optimize these variables simultaneously to build a well-connected multimodal network. The objective functions are calculated with equations (13)–(16) using these results:(13)TT=∑w=1W∑k=0K∑s=1Skwfk,swck,sw,(14)TC=∑i∈CN∑m∈MdcmψXihmori−Xihm+∑i∈CN∑m∈MccmψXihm−Xihmori+BP∑r=1Rmaxϑfsr−fsrori+oc∑r=1Rmaxlrfsr,(15)MS=∑m=1M∑l=1Lmfl,m−mfl,morimsl,m,(16)MS=∑hm∈Hexpfhmcaphm+θ+∑am∈Aexpfamcapam+θ.Equation (13) is the total travel time. Equation (14) is the total construction cost, which consists of the demolition and construction costs of the hubs, the bus purchase cost, and the bus operational cost. Equation (15) is the modal shift value. To analyse the modal shift of OD pairs with different lengths, mfl,m is denoted as the total flow of mode m among all OD pairs with distance l. The predefined weight parameter msl,m is used to value the 1 unit of flow change compared to the original network if msl,m is set as a minus value, which means the modal shift of this mode in this distance is encouraged to increase, and vice versa. Equation (16) is the service-level value, which is an exponential service-level function with the sum of the flow/capacity ratio of all the links and hubs. To ensure that highly loaded hubs or links will have higher service level value penalty, a parameter θ is added in equation (16).For location-based optimization, the constraints include those as follows:(17)∑i∈CN,m∈car,bicXihm=1,(18)HNmin≤∑i∈Hm,m∈car,bicXihm≤HNmax.Equation (17) ensures that for each candidate node, only one hub with mode m can be placed. Equation (18) is the hub number constraint for hubs with mode m.For route-based optimization, where the decision variables are bus layouts and frequencies, referring to [8, 73], the constraints include those as follows:(19)∑j∈HbusterXojr=1forr=1toRmax,(20)∑i∈HbusterXior=1forr=1toRmax,(21)∑i∈Hbusint,i≠jXijr−∑i∈Hbusint,i≠jXjir=0j∈Hbusint,forr=1toRmax,(22)∑i∈Hbusint,i≠jXijr≤1j∈Hbusint,forr=1toRmax,(23)∑j∈Hbusint,i≠jXijr≤1forj∈Hbusint,r=1toRmax,(24)Xjjr=0forj∈Hbusint,r=1toRmax,(25)Lminstop≤Xijrlijr≤Lmaxstopfori∈Hbusint,i≠j,j∈Hbusint,forr=1toRmax,(26)Lmin≤lr≤Lmaxforr=1toRmax,(27)fremin≤frer≤fremaxforr=1toRmax,(28)FSmin≤∑r=1Rmaxfsr≤FSmax,(29)frer=fsrlr2vbus.Among the constraints, equations (19) and (20) ensure that all bus routes start and end at a bus terminal. Equation (21) ensures that each bus interchange has one preceding node and one subsequent node. Equations (22)–(24) ensure that each bus node can be visited by one particular bus route only once at most (bus node repeat constraint). Equations (25)–(28) are the bus stop distance, bus route length, bus frequency, and bus fleet size constraints, respectively. The bus frequency of one route is calculated using the assigned number of buses in equation (29).


## 3.6. Solution Method of the Upper-Level Problem

The upper-level problem is a MIBLP with 4 sets of variables. To solve the problem, the HGSADC framework is introduced as follows.The local search and diversity management in HGSADC is problem specific. Thus, the local search with a crossover, mutation, and repair operator; diversity management with a Hamming distance-based objective; and a dynamic penalty were introduced one by one in the following sections.


### 3.6.1. Initialization Operator

For route-based optimization, the initialization operator determines the initialization locations of the bus stops and their sequences in each bus route. The initialization scheme should satisfy the stop distance constraint, route length constraint, and bus node repeat constraint.The bus route frequency is initialized by randomly assigning buses to each route according to equation (29). If the frequency violates frequency constraint and the number of buses violates fleet size constraint, the process will repeat until the constraints are satisfied. It should be noted that lines 4 ∼ 12 is called the bus sequence generation function in the rest of paper. This function can be used to generate a bus sequence with defined starting and ending nodes.For location-based optimization, the parking lot and bicycle station initialization should determine the number and position of each hub. The initialization operator randomly adds or deletes several hubs and ensures that the total number of hubs satisfies the hub number constraint.


### 3.6.2. Crossover Operator

The crossover of the bus route consists of two levels: combinations of different routes and different sequences of bus stops in one route. To increase the exploration capability and ensure the satisfaction of the constraints for the generated routes after crossover, referring to Ngamchai and Lovell [79], two crossover methods for route-based optimization, i.e., line crossover and stop crossover, are used. Line crossover exchanges two different routes in two solutions, while the stop crossover operator is triggered when two routes in one solution have segments that can be spliced.For location-based crossover, two hubs in two solutions are randomly selected; if they are different and not repeated in the other solution, then exchange.


### 3.6.3. Mutation Operator

The mutation process generates different genetic materials to facilitate exploration in the search space. According to Szeto and Wu [73], four types of mutation operators for route-based optimization are proposed: insert, delete, swap, and transfer.The insert mutation is used to insert an interchange node in a randomly selected route. The deletion mutation is similar. The swap mutation is to exchange the nodes of two routes if the stop distance constraint can be satisfied after the swap, while transfer mutation is to transfer one node to the other route, and similarly, the stop distance constraint should be obeyed after the transfer.The location-based mutation only contains the insert and delete mutations. The insert mutation inserts a random hub that is not in the solution, whereas the delete mutation deletes one random hub in the solution.


### 3.6.4. Repair Operator

After the crossover and mutation processes, the generated children may violate the constraints, and a repair operator is needed to fix these infeasible solutions. The constraints for route-based optimization are route length, stop distance, node repeat, and frequency, while the constraint for location-based optimization is the number of hubs. With four constraints, there is a high probability that the bus route will violate other constraints after one constraint is repaired. A new four-stage repair strategy was proposed in this study. The solution set will undergo four repair operators in the sequence of the route length operator, stop distance operator, node repeat operator, and route frequency operator.(1)
Route Length Repair Operator. See Figure 2, when route r1 is too short, the algorithm selects two nodes, n1a and n1b, with a sufficient distance and inserts an intermediate node n1c between them. Then, two stop sequences from n1a to n1c and from n1c to n1b are generated with the bus stop sequence generation function in Algorithm1, and the original stop sequence between n1a and n1b is replaced. The situation when the route is too long is similar, and two nodes with a sufficient distance are chosen and replaced with the generated stop sequence using the bus stop sequence generation function.(2)
Stop Distance Repair Operator. A stop sequence that violates the stop distance constraint may contain only two nodes or consist of several consecutive nodes in which all the adjacent nodes violate the stop distance constraint. When the stop distance is too long, and the stop sequence has only two nodes, see Figure 3, the algorithm finds an intermediate node to insert between the two nodes. When the stop sequence has over two nodes, a new stop sequence will be generated to replace it. The situation when stop distance is too short is similar.(3)
Route Repeat Repair Operator. For repeated nodes, one node is randomly selected to be unchanging. Then, for each repeated node, n1a, a node is selected from CN that satisfies the stop distance constraint with both n1a − 1 and n1a + 1 to replace n1a.(4)
Route Frequency Repair Operator. The frequency repair operator adds or deletes one bus until the bus frequency constraint is satisfied.Algorithm 1: Bus route initialization.
(1)
for route r from 1 to Rmax(2)
Whether router needs to construct a new terminal with probability Pters and Ptere(3)
Determine the starting and ending node: if new terminal should be built, select a node fromCN as the starting node nrs or ending node nre of route r; if no, then select a node from TN as the starting node or ending node(4)
while the tth node of roue r, nrt is not the nre, and the route satisfies the route length constraint, bus node repeat constraint, and time < Tmaxdo(5)
Find the node setNrstopt with distance to nrt satisfying the stop distance constraint(6)
ifNrstopt is empty or all nodes have been visited in Nrstopt(7)
nrt=0, t = t-1(8)
elseifNrstopt is not empty(9)
Randomly select a nodenrstopt from Nrstopt that has not been visited(10)
ifnre−nrtempt<nre−nrtnrt+1=nrtemptelseifnre−nrtempt≥nre−nrtMarknrstopt has been visited and go back to line 6end(11)
if route r is feasible thenInsert router into initialization solution(12)
elseif route r is infeasible thenGo back to line 2.end


### 3.6.5. Diversification Mechanism and Infeasible Solution Management

To ensure the diversity of the solutions, the Hamming distance is introduced to represent the differences between one solution and the best solution, which is also called the diversity contribution functionΔS. It is defined as the average distance from S to its closest neighbours in the subpopulation. ΔS consists of three parts, which are the bus route, parking lot, and bicycle station Hamming distance. According to Szeto and Wu [73], the bus Hamming distance is the number of different consecutive node pairs for each pair of corresponding routes compared to the best solution. For the location-based one, it is the number of different hubs compared to the best solution, and ΔS is the sum of the above three distances.Figure 2
Situation 1 in the length repair operator.Figure 3
Situation 1 of the stop distance repair operator.A dynamic penalty is then introduced to handle the infeasible solutions. Suppose the number of feasible and infeasible solutions areNinf and Nfea, respectively. A basic penalty value, ϖ, is adjusted for every μ iterations, and when the number of Ninf increases by 1, ϖ decreases by λ%, and vice versa. (see Algorithm 2)Algorithm 2: HGSADC.
(1)
Initialization(2)
while number of iterations without improvements < ItNIdo(3)
while maximum solution size is not reached doselect parent solutionsS1 and S2, create offspring O using route-based and location-based crossover operator. Insert O into solution set (crossover)end while(4)
select solutionS1 and mutate with route-based and location-based mutation operator (mutation)(5)
for each solutionS1, if infeasible then(6)
RepairS1 with probability Prep (repair)if repairable theninsert into feasible subpopulation.elseinsert into infeasible subpopulationend if(7)
elseifS1 feasible theninsert into feasible subpopulationend if(8)
calculate objective values with GMS/TAP(9)
calculate the biased fitness valueSelect survivors(selection)(10)
if best solution not improved for Itdiv iterations, thenDiversify population(11)
Adjust penalty parameters for infeasibilityend whileReturn best feasible solutionThen, the survivors of each iteration are selected based on the rank of objective functionfitS and the rank of Hamming distance hdrS, which is called the biased fitness in this paper. For infeasible solutions, the objective function should multiply dynamic penalty values, ϖ. Referring to Vidal et al. [76], the biased fitness, BFS, is formulated:(30)BFS=fitS+1−NeliteNmaxhdrS,where Nelite is a parameter ensuring elitism properties during survivor selection.After building the whole bilevel model, the extracts of the algorithms can be found in Appendix B.


## 3.6.1. Initialization Operator

For route-based optimization, the initialization operator determines the initialization locations of the bus stops and their sequences in each bus route. The initialization scheme should satisfy the stop distance constraint, route length constraint, and bus node repeat constraint.The bus route frequency is initialized by randomly assigning buses to each route according to equation (29). If the frequency violates frequency constraint and the number of buses violates fleet size constraint, the process will repeat until the constraints are satisfied. It should be noted that lines 4 ∼ 12 is called the bus sequence generation function in the rest of paper. This function can be used to generate a bus sequence with defined starting and ending nodes.For location-based optimization, the parking lot and bicycle station initialization should determine the number and position of each hub. The initialization operator randomly adds or deletes several hubs and ensures that the total number of hubs satisfies the hub number constraint.


## 3.6.2. Crossover Operator

The crossover of the bus route consists of two levels: combinations of different routes and different sequences of bus stops in one route. To increase the exploration capability and ensure the satisfaction of the constraints for the generated routes after crossover, referring to Ngamchai and Lovell [79], two crossover methods for route-based optimization, i.e., line crossover and stop crossover, are used. Line crossover exchanges two different routes in two solutions, while the stop crossover operator is triggered when two routes in one solution have segments that can be spliced.For location-based crossover, two hubs in two solutions are randomly selected; if they are different and not repeated in the other solution, then exchange.


## 3.6.3. Mutation Operator

The mutation process generates different genetic materials to facilitate exploration in the search space. According to Szeto and Wu [73], four types of mutation operators for route-based optimization are proposed: insert, delete, swap, and transfer.The insert mutation is used to insert an interchange node in a randomly selected route. The deletion mutation is similar. The swap mutation is to exchange the nodes of two routes if the stop distance constraint can be satisfied after the swap, while transfer mutation is to transfer one node to the other route, and similarly, the stop distance constraint should be obeyed after the transfer.The location-based mutation only contains the insert and delete mutations. The insert mutation inserts a random hub that is not in the solution, whereas the delete mutation deletes one random hub in the solution.


## 3.6.4. Repair Operator

After the crossover and mutation processes, the generated children may violate the constraints, and a repair operator is needed to fix these infeasible solutions. The constraints for route-based optimization are route length, stop distance, node repeat, and frequency, while the constraint for location-based optimization is the number of hubs. With four constraints, there is a high probability that the bus route will violate other constraints after one constraint is repaired. A new four-stage repair strategy was proposed in this study. The solution set will undergo four repair operators in the sequence of the route length operator, stop distance operator, node repeat operator, and route frequency operator.(1)
Route Length Repair Operator. See Figure 2, when route r1 is too short, the algorithm selects two nodes, n1a and n1b, with a sufficient distance and inserts an intermediate node n1c between them. Then, two stop sequences from n1a to n1c and from n1c to n1b are generated with the bus stop sequence generation function in Algorithm1, and the original stop sequence between n1a and n1b is replaced. The situation when the route is too long is similar, and two nodes with a sufficient distance are chosen and replaced with the generated stop sequence using the bus stop sequence generation function.(2)
Stop Distance Repair Operator. A stop sequence that violates the stop distance constraint may contain only two nodes or consist of several consecutive nodes in which all the adjacent nodes violate the stop distance constraint. When the stop distance is too long, and the stop sequence has only two nodes, see Figure 3, the algorithm finds an intermediate node to insert between the two nodes. When the stop sequence has over two nodes, a new stop sequence will be generated to replace it. The situation when stop distance is too short is similar.(3)
Route Repeat Repair Operator. For repeated nodes, one node is randomly selected to be unchanging. Then, for each repeated node, n1a, a node is selected from CN that satisfies the stop distance constraint with both n1a − 1 and n1a + 1 to replace n1a.(4)
Route Frequency Repair Operator. The frequency repair operator adds or deletes one bus until the bus frequency constraint is satisfied.Algorithm 1: Bus route initialization.
(1)
for route r from 1 to Rmax(2)
Whether router needs to construct a new terminal with probability Pters and Ptere(3)
Determine the starting and ending node: if new terminal should be built, select a node fromCN as the starting node nrs or ending node nre of route r; if no, then select a node from TN as the starting node or ending node(4)
while the tth node of roue r, nrt is not the nre, and the route satisfies the route length constraint, bus node repeat constraint, and time < Tmaxdo(5)
Find the node setNrstopt with distance to nrt satisfying the stop distance constraint(6)
ifNrstopt is empty or all nodes have been visited in Nrstopt(7)
nrt=0, t = t-1(8)
elseifNrstopt is not empty(9)
Randomly select a nodenrstopt from Nrstopt that has not been visited(10)
ifnre−nrtempt<nre−nrtnrt+1=nrtemptelseifnre−nrtempt≥nre−nrtMarknrstopt has been visited and go back to line 6end(11)
if route r is feasible thenInsert router into initialization solution(12)
elseif route r is infeasible thenGo back to line 2.end


## 3.6.5. Diversification Mechanism and Infeasible Solution Management

To ensure the diversity of the solutions, the Hamming distance is introduced to represent the differences between one solution and the best solution, which is also called the diversity contribution functionΔS. It is defined as the average distance from S to its closest neighbours in the subpopulation. ΔS consists of three parts, which are the bus route, parking lot, and bicycle station Hamming distance. According to Szeto and Wu [73], the bus Hamming distance is the number of different consecutive node pairs for each pair of corresponding routes compared to the best solution. For the location-based one, it is the number of different hubs compared to the best solution, and ΔS is the sum of the above three distances.Figure 2
Situation 1 in the length repair operator.Figure 3
Situation 1 of the stop distance repair operator.A dynamic penalty is then introduced to handle the infeasible solutions. Suppose the number of feasible and infeasible solutions areNinf and Nfea, respectively. A basic penalty value, ϖ, is adjusted for every μ iterations, and when the number of Ninf increases by 1, ϖ decreases by λ%, and vice versa. (see Algorithm 2)Algorithm 2: HGSADC.
(1)
Initialization(2)
while number of iterations without improvements < ItNIdo(3)
while maximum solution size is not reached doselect parent solutionsS1 and S2, create offspring O using route-based and location-based crossover operator. Insert O into solution set (crossover)end while(4)
select solutionS1 and mutate with route-based and location-based mutation operator (mutation)(5)
for each solutionS1, if infeasible then(6)
RepairS1 with probability Prep (repair)if repairable theninsert into feasible subpopulation.elseinsert into infeasible subpopulationend if(7)
elseifS1 feasible theninsert into feasible subpopulationend if(8)
calculate objective values with GMS/TAP(9)
calculate the biased fitness valueSelect survivors(selection)(10)
if best solution not improved for Itdiv iterations, thenDiversify population(11)
Adjust penalty parameters for infeasibilityend whileReturn best feasible solutionThen, the survivors of each iteration are selected based on the rank of objective functionfitS and the rank of Hamming distance hdrS, which is called the biased fitness in this paper. For infeasible solutions, the objective function should multiply dynamic penalty values, ϖ. Referring to Vidal et al. [76], the biased fitness, BFS, is formulated:(30)BFS=fitS+1−NeliteNmaxhdrS,where Nelite is a parameter ensuring elitism properties during survivor selection.After building the whole bilevel model, the extracts of the algorithms can be found in Appendix B.


## 4. Results and Discussion

This section discusses three experiments that were performed to validate the proposed models and algorithms. In experiment 1, the exact results of the standard NLP of the GMS/TAP were verified, along with the nonexact results of the hybrid MSA algorithm. In experiment 2, the results of the upper-level problems were validated by testing the diversification mechanism and biased fitness function. The application of the IMNDP in a real-size network was then tested in experiment 3.


### 4.1. Experiment 1

The first experiment was conducted to show the solution properties of the GMS/TAP. The nine-node network used in Wu et al. [82] was adopted, as shown in Figure 4. The link impedance functions of the road links and bus links were BPR (Bureau of public road) functions. The link free-flow travel time and capacity were noted in the two-tuple next to each link.Figure 4
Nine-node network.To adopt the original network into GMS/TAP, a modified network with the same properties was proposed. As shown in Figure5, numbers 1–8 are the road nodes, 9–17 are the bus nodes, and 18–21 are the OD nodes. The transfer links and embarking/alighting links were added at a cost of zero. The road and bus subnetwork links were maintained with the same parameters. Meanwhile, although the car + bus combined mode is difficult to observe in real life, they were allowable in this experiment, and car users could directly transfer from a car node to a bus stop, and vice versa.Figure 5
Modified nine-node network.For simplification, thecar and bus modes are denoted as c and b. Although it was unnecessary for a path set to cover all possible modes, in this experiment, six modes were selected to show the solution properties of the GMS/TAP. The corresponding modes included those as follows: (1) 0 transfers, where the travel mode is either c or b; (2) 1 transfer, where the travel modes are either b + c or c + b, and (3) 2 transfers, where the travel modes comprise b + c + b and c + b + c. Three paths for each mode with the lowest path costs were selected. If there were fewer than three paths, only the available paths were selected.The model was calculated when maximum transfer numberK = 0, 1, and 2 using the MATLAB fmincon solver. Meanwhile, when K = 2, the results were also calculated using the hybrid MSA algorithm. The logit parameters at each level were α−1=0.5, β−1=1, α0=1, β0=0.9, α1=0.5, β1=0.85, α2=0.5, and β2=0.8. The path set for is presented in Table 3. A bold number in the path column represents the transfer hub.Table 3
Path set of experiment 1.
OD: 18–20NumberModePath1c18, 1, 5, 7, 3, 20218, 1, 6, 5, 7, 3, 20318, 1, 5, 7, 8, 3, 204b18, 10, 12, 14, 16, 20518, 10, 12, 9, 14, 16, 206b + c18, 10,1, 5, 7, 3, 20718, 10, 12,5, 7, 3, 20818, 10, 12, 14,7, 3, 209c + b18, 1, 5, 7,14, 16, 201018, 1, 5, 7, 3,16, 201118, 1, 5,12, 14, 16, 2012b + c + b18, 10,1, 5, 7, 14, 16, 201318, 10,1, 5, 7, 3, 16, 201418, 10, 12,5, 7, 14, 16, 2015c + b + c18, 1,10, 12, 5, 7, 3, 201618, 1, 5,12, 14, 7, 3, 201718, 1,10, 12, 14, 7, 3, 20OD: 18–2118c18, 1, 5, 7, 4, 211918, 1, 5, 7, 8, 4, 212018, 1, 6, 5, 7, 4, 2121b18, 10, 12, 9, 15, 17, 2122b + c18, 10,1, 5, 7, 4, 212318, 10, 12,5, 7, 4, 212418, 10, 12, 14,7, 4, 2125c + b18, 1, 5, 7, 4,17, 212618, 1, 5, 7, 8,15, 17, 212718, 1, 5, 7, 8, 4,17, 2128b + c + b18, 10,1, 5, 7, 4, 17, 212918, 10, 12,5, 7, 4, 17, 213018, 10, 12, 14,7, 4, 17, 2131c + b + c18, 1,10, 12, 5, 7, 4, 213218, 1, 5,12, 14, 7, 4, 213318, 1,10, 12, 14, 7, 4, 21OD: 19-2034c19, 2, 5, 7, 3, 203519, 2, 5, 7, 8, 3, 203619, 2, 6, 5, 7, 3, 2037b19, 11, 13, 9, 14, 16, 2038b + c19, 11,2, 5, 7, 3, 203919, 11,2, 5, 7, 8, 3, 204019, 11, 13,6, 5, 7, 3, 2041c + b19, 2, 5, 7,14, 16, 204219, 2, 5, 7, 3,16, 204319, 2, 5,12, 14, 16, 2044b + c + b19, 11,2, 5, 7, 14, 16, 204519, 11,2, 5, 7, 3, 16, 204619, 11,2, 5, 12, 14, 16, 2047c + b + c19, 2, 5,12, 14, 7, 3, 204819, 2, 5, 7,14, 16, 3, 204919, 2, 5,12, 14, 16, 3, 20OD: 19–2150c19, 2, 5, 7, 4, 215119, 2, 5, 7, 8, 4, 215219, 2, 6, 5, 7, 4, 2153b19, 11, 13, 15, 17, 215419, 11, 13, 9, 15, 17, 2155b + c19, 11,2, 5, 7, 4, 215619, 11,2, 5, 7, 8, 4, 215719, 11, 13,6, 5, 7, 4, 2158c + b19, 2, 5, 7, 4,17, 215919, 2, 5, 7, 8,15, 17, 216019, 2, 5, 7, 8, 4,17, 2161b + c + b19, 11,2, 5, 7, 4, 17, 216219, 11,2, 5, 7, 8, 15, 17, 216319, 11,2, 5, 7, 8, 4, 17, 2164c + b + c19, 2, 5,12, 14, 7, 4, 216519, 2, 5,12, 14, 7, 8, 4, 216619, 2, 5, 7,8, 15, 17, 4, 21The optimal flow,F∗, and optimal objective, ξ∗, makes the NLP have the optimal objective, and equilibrium path cost when K = 2 is based on the fmincon solver, and the hybrid MSA are listed in Tables 4 and 5.Table 4
Results of the fmincon solver.
OD: 18–20OD: 18–21OD: 19-20OD: 19–21No.FlowCostNo.FlowCostNo.FlowCostNo.FlowCost1059.618064.7342.652.3500.157.421.358.419064.735059.3513.057.43066.7208.163.536061.352066.44081.7212.564.8378.750.95315.755.5527.455.0222.364.8382.252.354060.760.159.623086.839059.3550.757.47081.724085.240085.0561.957.48080.1250.164.741053.957090.19061.226071.7422.252.3580.457.410059.6272.264.7430.152.359064.311059.6282.864.744053.8602.157.412061.229086.8452.652.3611.257.4130.259.630085.1460.152.362064.314083.331086.84711.550.6631.957.415081.73212.063.148053.9643.355.7161.058.033085.249052.3659.755.717080.166064.3ξ3.750911.0166−4.47420.6910Objective function8.6484e − 6CPU time6.925 sTable 5
Results of the hybrid MSA.
OD: 18–20OD: 18–21OD: 19-20OD: 19–21No.FlowCostNo.FlowCostNo.FlowCostNo.FlowCost1059.418064.6343.052.2500.857.321.358.419064.635059.3513.057.33066.5208.463.536061.352066.44081.8212.664.8378.851.05315.855.7527.354.8221.964.6381.952.254060.960.159.523086.739059.3550.457.37081.624085.340085.1561.657.38080.1250.464.641053.857090.29061.126071.6421.852.2580.457.3100.159.5271.564.643052.459064.411059.6283.264.644053.8601.657.312061.129086.7453.152.2610.857.3130.359.430085.346052.462064.414083.231086.74711.350.8632.857.315081.63212.063.148053.8642.655.9161.058.033085.349052.46510.355.917080.166064.4ξ3.745410.7520−4.74020.4061Objective function3.5948CPU time0.145 sIn Table4, the optimal objectives based on the fmincon solver is 8.6484e − 6, which is close to the coverage criterion zero of Proposition 2. Meanwhile, for all the paths with the same mode, such as paths 18, 19, and 20 with mode c, the paths carrying flow all have the lowest path cost, so the UE principle is satisfied.The other two conditions in section 3.3.1 are also verified. Taking OD 18–20 as an example, the minimum utility values of 0, 1, and 2 transfers are 55.0, 59.6, and 58.0, respectively, and the equilibrium flows are 28.7, 0.1, and 1.2, which obeys condition 1. Among the two transfer trips, the minimum utility values are 59.6 and 58 for modesb + c + b and c + b + c, respectively, and the total flow of 1.2 is divided into 0.2 and 1.0, which also obeys condition 2.Table5 shows that the objective function is 3.5948, and the UE principle is also obeyed. Compared with Table 4, the flow patterns are similar, while the CPU times are 0.145 s and 6.925 s. These results show that the hybrid MSA could obtain an approximate result with a much faster speed.


### 4.2. Experiment 2

#### 4.2.1. Network and Parameters

Experiment 2 tested the performance of the HGSADC, as shown in Figure6. The three subnetworks from left to right are the road, metro, and bus subnetworks. The link and hub free-flow travel times and capacities are denoted by blue and green texts.Figure 6
Network of experiment 2.The following modifications were made based on the network of experiment 1:(1)
A new parking lot was added in node 26 and connected to road node 5 with bidirectional auxiliary links 37 and 51. Users must transfer from a parking lot to a metro station to execute a P&R combined trip.(2)
Two parallel metro lines with a frequency of 0.083 h/veh were added to provide more mode choices. Only one bus route 21–25 was kept at 0.068 h/veh to leave space for optimization.(3)
In experiment 1, the costs of transfer links, embarking/alighting links, and auxiliary links were default to be 0. In this experiment, to verify the results after the optimization of the hubs, the costs of transfer links, embarking/alighting links, and auxiliary links were added with the sum of the walking time of the links, walking time inside the hubs (BPR function with green text as the parameters), waiting time (the route frequency), and pure penalty (a fixed time of 0.083 h).(4)
A Euclidean space coordinate system was built, where thex and y axes represented the position of the node (see Figure 7). Because the four OD pairs were connected to several direct car, bus, or metro services to avoid the situation where no users take combined trips, one unit of length in the network was set at 10 km, and the car, metro, and bus speeds were set at 30, 25, and 20 km/h, respectively.Experiment 2 aimed to optimize the parking lots and bus routes to integrate the network. The maximum number of bus routes was set to two, while the number of parking lots was set to 1–6. The original number of buses was five. The weights of the four objectives such as the total time, construction cost, modal shift value, and service level value, were set at 1, 1, 2, and 3, respectively.Figure 7
Coordinates of experiment 2 subnetworks.
(a)(b)(c)(d)


#### 4.2.2. Results of Experiment 2

Four sets of GAs tested the model performance: GA1 was the normal GA, GA2 was a GA with a diversification mechanism alone, GA3 was a GA with a biased fitness function alone, and GA4 was the HGSADC. Each GA was ended when 500 generations were produced and repeated for 20 runs. The average results of the original network and 4 GAs are listed in Table6.Table 6
Average results of 4 GAs.
GA1GA2GA3HGSADCWeighted objective21.105817.364117.364117.3641Travel time292.3744289.5535289.5535289.5535Construction cost33016.7230084.7230084.7230084.72Modal shift value−19.1918−28.6002−28.6002−28.6002Service level value451.1113405.3683405.3683405.3683Convergence iterations59.7111.392.447.2CPU time (s)97.1122.5103.4193.4In all 20 runs, GA2, GA3, and HGSADC all converged to the same results with a weighted objective value of 17.3641. However, GA1 was trapped in a local optimal objective with an average objective value of 21.1058, showing that both the diversification mechanism and biased fitness function contributed to improving the exploration capacity. Meanwhile, HGSADC had far fewer convergence iterations than the other three GAs. Because the diversification mechanism involved a reinitialization process to obtain more genetic materials, it also made the average CPU time longer than the other three GAs. Considering it is the total time of 500 iterations, the fast convergence speed can be useful in large-scale problems.


#### 4.2.3. Analyzation of Schemes

The performances of the original, GA1, and HGSADC schemes were then compared. See Table7.Table 7
Results of three schemes.
SchemeParking lotBus routesNumber of busesBus frequencyTotal timeConstruction costModal shiftService levelOriginal521-22-23-24-2550.0678296.2606——478.1152GA15, 626-27-28–2920.2494288.222925187.08−17.5560443.090621-22-23-24-2530.113HGSADC526-27-28-29-3030.1130282.180030084.72−28.6002405.368321-22-23-24-2520.1694(1) Original Network. The flow distribution of the original network is shown in Figure 8. From left to right, the multimodal network comprises a green road subnetwork, brown metro subnetwork, and a purple bus subnetwork. The link flow is noted beside the link in black text, while the hub flow is shown in blue text. Meanwhile, transfer flow is noted beside the transfer link. The pink text shows the mtr + car and mtr + bus transfer flow, while the orange text shows the car + mtr and bus + mtr transfer flow.Figure 8
Flow distribution of the original network.The node flow is the sum of all the links entering the node. The two parallel OD pairs, 9–11 and 10–12, can be served by the car and two parallel metro lines directly. The speed of the car is faster than that of the metro, while the capacity of the road link is 10, which is not sufficient to meet the 30 demands of the parallel OD pairs. This portion of the users will choose the single-car mode first and then the metro. Meanwhile, because a parking lot is built at node 26, and road link 1-5-7-3 is congested (the flow/capacity ratio is higher than 1), while metro line 13-14-15-16 is underutilized (lower than 0.75), some users will transfer from parking lot 26 to metro station 14 to bypass the congested road link, 5-7-3.For the two diagonal OD pairs, OD pair 9–12 can be served by the bus route directly. Among the 30 demands, 10.2 users choose bus link 21-22, and the other 8.9 flow of this route is themtr + bus flow. The bus route relieves the congestion of the road links in this direction, and road links 1–6, 5-6, 7-8, and 7-4 all have a good service level. However, the other diagonal OD pair, 10-11, can only be served by the road subnetwork directly, which makes the service level of related links 2–5 and 8-3 very low.(2) GA1 Scheme. The GA1 scheme adds a new parking lot in node 31, as well as bus route 26-27-2829 with frequency 0.12 h/veh. The original parking lot 30 and bus route 21-22-23-24-25 are kept, but the frequency of the original bus route increases to 0.112 h/veh. The flow distribution is illustrated in Figure 9.Figure 9
Flow distribution of GA1 scheme.The newly added bus route, 26–29, only attracts 3.7 users because it is parallel to metro route 13–16. However, the congestion of road links 1-5-7-3 is still slightly reduced. Meanwhile, the new parking lot at 31 attracts 14.6 users to use the P&R mode, which reduces the congestion of road link 6-8-4.(3) HGSADC Scheme. The HGSADC scheme adds a symmetrical bus route, 26-27-28-29-30, with a frequency of 0.113 h/veh. The parking lot of the original network is kept. See Figure 10.Figure 10
Flow distribution of HGSADC scheme.The added bus route provides a direct bus service for OD pair 10-11, which attracts eight single bus mode users to bus stop 26. This bus route creates two new mode choices for this OD pair,mtr + bus and bus + mtr, and except for 8 direct users, 8.1 users transfer from metro station 18 to bus stop 27 and 2.2 users transfer from bus stop 29 to metro station 13, which significantly relieves the congestion of road links 2–5 and 8-3. With more users choosing PT, the flows of other single-car mode paths, such as 2-5-7-3 and 2-6-8-3, also decrease by a certain amount, which leads to an overall reduction in congestion of road subnetwork. Because the weight of the service level value is the highest, the lowest service level value of 405.4 makes the weighted objective of the HGSADC lower than that of the GA1 scheme even though the total time is similar, and the construction cost is higher.


### 4.3. Experiment 3

#### 4.3.1. Original Network

The model was then tested in a real-size network of Jianye District, Nanjing. See Figure11 (Open street map, 2022). The research area covers 80.87 square kilometers, with a population of 489,800. There are two metro lines running through this district and their intersection node, and Yuantong station lies in the central business district (CBD) of the research area. The yellow isolated node in the upper left corner is Linjiang Station. Since Nanjing is divided into two parts by the Yangtze River, Linjiang Station is an important node which connects the south and north banks of the city, and a typical P&R parking lot is built there.Figure 11
Research area.The network comprises 5 subnetworks, i.e., bicycle, car, car hailing, metro, and bus subnetwork. In Figure12, from (a) to (e), the blue, green, and red subnetworks are bicycle, car, and car hailing subnetworks, and they share the same physical network with 125 nodes and 428 links; meanwhile, car and car hailing subnetwork share the same flow. The blue bicycle subnetwork contains 55 bicycle stations while the green car subnetwork has 45 parking lots. The brown metro subnetwork has 2 metro lines and 18 metro stations. The purple bus subnetwork consists of 11 bus lines and 95 bus stops. The district is divided into 22 traffic zones with 462 OD pairs, see Figure 12(f).Figure 12
Original subnetworks. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork. (f) OD pair.
(a)(b)(c)(d)(e)(f)


#### 4.3.2. Results of Experiment 3

The decision variables were bicycle stations, parking lots, bus layouts, and their frequencies. The iteration number was set as 3000, and the maximum transfer numberK was 2. The number of bus routes was kept at 11. The number of bicycle stations, parking lots, and bus stops was set from 50 to 200, 45 to 100, and 80 to 200, respectively. The bus stop distance constraint was from 0.3 km to 2 km and 0.0833 h/veh to 0.5 h/veh for bus route frequency constraint. The bus route length constraint was from 5 km to 25 km. 40 paths for each OD pair were chosen to ensure the stability of hybrid MSA in obtaining the equilibrium flow. The pure penalty is 5 mins.(1) Travel Time and Construction Cost. The algorithm converges in 2335 iterations with 172417.54 s. GMS/TAP function accounts for 77% of the total time, in which path generation, path utility, and hybrid MSA takes 39.5%, 4.8%, and 32.7%. The HGSADC and other functions take 23% of the total time. The best solution contains 51 parking lots, 74 bicycle stations, and 11 bus routes with 117 bus stops. See Figure 13. With a construction cost of 11566314.32 RMB, the travel time decreases from 749042.8 h to 680345.5 h, which improves 9.2%. Since this study concentrates on the modal shift and balanced use of the whole network, comparing with Szeto et al. [69], in which a 22.7% reduction of travel time in transit subnetwork is obtained, the improvement of travel time is not obvious in this experiment.Figure 13
Subnetworks of the best solution. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Bus subnetwork.
(a)(b)(c)(2) Service Level Value. The service level of the original network is presented in Figure 14, and in the original network, the road and metro subnetworks are overloaded, especially in the CBD area with longitude from 118.7 to 118.72, latitude from 31.98 to 31.99. Meanwhile, the flow between north bank and south bank also concentrates in the car and metro mode, and the service level in the corresponding links and hubs are all lower than level E. With a congested metro and car subnetwork, the bicycle and bus subnetworks in Figures 14(a) and 14(e) are underutilized; to fully use the potential capacity of these subnetworks, a modal shift is necessary.Figure 14
Service level of original network. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork.
(a)(b)(c)(d)(e)In the optimal scheme (see Figure15), the service level value decreases from 40253.82 to 29903.51, which improves 25.71%. Although the main traffic corridor from west to east is still congested in road subnetwork, the number of over congested links and hubs (service level higher than 1) decreases with 21 and 3 for road and metro subnetworks, respectively. These flows transfer to bicycle and bus subnetworks, and an increase in usage with 36 and 44 links and hubs is obtained.Figure 15
Service level of best solution. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork.
(a)(b)(c)(d)(e)(3) Modal Shift Value. The modal shift value is −131305.2. After deleting the modes with total flow lower than 2000, the mode flow of single, one transfer, and two transfers modes are shown in Figures 16(a)–16(c), respectively. In the x axis, the mode bic, car, ch, mtr, and bus are denoted as b, c, h, m, and B, respectively. The combined trip, such as car + mtr is denoted as cm. While on the y axis, the mode flow is classified per kilometers according to the direct OD distance. The z axis is the amount of mode flow.Figure 16
Mode flow of original network. (a) Single-mode flow. (b) One transfer mode flow. (c) Two transfers mode flow.
(a)(b)(c)In the original network, 77.33% of the flow is single-mode within which 37.17% comes from the private car mode. Since there are only two metro lines while a complete bus network is built, the metro subnetwork is mostly used in combined trips. The bicycle is normally used alone in short distance trips and in long-distance trips to complement bus and metro subnetworks, and the car hailing subnetwork plays a similar role as the private car but with a lower flow.Considering the original road subnetwork is over congested, the single-modec is encouraged to be reduced with positive weight value. Meanwhile, one transfer modes in short distance and 2 transfers modes in short and middle distance needs to be cut down as well. The modal shift weight matrix is divided into seven columns and three rows, see Table 8.Table 8
The modal shift weight matrix.
bcmhBOne transfer modeTwo transfers modeShort distance (0∼4 km)−33−1−1−212Middle distance (5∼8 km)−22−1−1−3−11Long distance (9∼11 km)−11−1−1−3−2−2The flow of each mode in the optimal scheme is shown in Figure17. The 0 transfer modes and two transfers modes users decrease with 6885 and 12797, respectively, and these users transfer to one transfer mode. Among single-mode trips, a decrease of 27768 car users and an increase of 4853, 14860, and 20853 bic, mtr, and bus users are observed. These four single-modes produce a weighted modal shift value −163681. In one transfer trips, the modal shift does not change as wishes, the increase in 19682 one transfer modes users concentrate in short distance OD pairs, especially the bm and mb modes. But in other perspective, this means the new bicycle station scheme improves the connection between bicycle and metro subnetworks. For two transfers modes, the low flow does not have an obvious influence on the modal shift value.Figure 17
Mode flow of the solution. (a) Single-mode flow. (b) One transfer mode flow. (c) Two transfers mode flow.
(a)(b)(c)


#### 4.3.3. Sensitive Tests of the Parameters

(1) Effect of the Transfer Penalty. This section tested the effect of transfer penalty on the flow distribution of the optimized network. Since there are several types of transfer costs which are related to the connected modes, it is difficult to test these parameters one by one. Considering every transfer involves in a pure penalty with 5 mins, pure penalty represents the natural resistance of the users to the transfer. Thus, the experiment tested the network with a pure penalty from 0 to 15 min, with an interval of 0.25 min. The total mode miles are chosen as the index, which is used to describe the overall usage of the corresponding subnetwork. See Figure 18.Figure 18
Effect of pure penalty on the flow distribution. (a) Effect of pure penalty on the subnetworks. (b) Effect of pure penalty on the number of transfer.
(a)(b)In Figure18(a), when the transfer penalty increases, the miles on car and bicycle subnetwork increases. For the metro and bus subnetworks, which rely more on transfer, the increase of transfer penalty leads to a decrease of their travel miles. In Figure 18(b), we can see that the number of transfers is highly related to the transfer time. Intuitively, the decrease of the number of transfers is good for the transportation system. However, in this test, this is achieved by increasing the transfer time rather than improving the connection of the network. The long-time of waiting pushes the users to the congested links of the single modes, while they cannot use the underutilized links of other subnetworks through combined trip to bypass the congested links. The results suggest that the number of transfers alone may not be a good index to value the integration of the multimodal network.(2) Effect of the Congestion. We further investigated the congestion effect under different demand levels. See Figure 19. The OD demand varies from 5% to 300% of the original demand with an increase of 5%. As the OD demand increases, the users choose road and metro subnetworks at first. When the demand reaches 45% of the original demand, the car subnetwork becomes relatively more congested than other three subnetworks, i.e., bus, bicycle, and car hailing, and their share of travel miles increases while car and metro decreases. As shown in Figure 15, under the original demand, the car and metro subnetworks are already over congested, the decrease speed for these two modes is increasing, and more users choose other three modes.Figure 19
Effect of congestion on the flow distribution.This result suggests that there may be one or two modes that the users will choose at first because their competitions are much better than other modes. Only when these preferred subnetworks are congested, then users will choose the remaining modes. The natural differences of the travel modes are hard to change only via an improvement of the infrastructures. A systematic demand control and operational management are also required to balance the flow distribution.


## 4.1. Experiment 1

The first experiment was conducted to show the solution properties of the GMS/TAP. The nine-node network used in Wu et al. [82] was adopted, as shown in Figure 4. The link impedance functions of the road links and bus links were BPR (Bureau of public road) functions. The link free-flow travel time and capacity were noted in the two-tuple next to each link.Figure 4
Nine-node network.To adopt the original network into GMS/TAP, a modified network with the same properties was proposed. As shown in Figure5, numbers 1–8 are the road nodes, 9–17 are the bus nodes, and 18–21 are the OD nodes. The transfer links and embarking/alighting links were added at a cost of zero. The road and bus subnetwork links were maintained with the same parameters. Meanwhile, although the car + bus combined mode is difficult to observe in real life, they were allowable in this experiment, and car users could directly transfer from a car node to a bus stop, and vice versa.Figure 5
Modified nine-node network.For simplification, thecar and bus modes are denoted as c and b. Although it was unnecessary for a path set to cover all possible modes, in this experiment, six modes were selected to show the solution properties of the GMS/TAP. The corresponding modes included those as follows: (1) 0 transfers, where the travel mode is either c or b; (2) 1 transfer, where the travel modes are either b + c or c + b, and (3) 2 transfers, where the travel modes comprise b + c + b and c + b + c. Three paths for each mode with the lowest path costs were selected. If there were fewer than three paths, only the available paths were selected.The model was calculated when maximum transfer numberK = 0, 1, and 2 using the MATLAB fmincon solver. Meanwhile, when K = 2, the results were also calculated using the hybrid MSA algorithm. The logit parameters at each level were α−1=0.5, β−1=1, α0=1, β0=0.9, α1=0.5, β1=0.85, α2=0.5, and β2=0.8. The path set for is presented in Table 3. A bold number in the path column represents the transfer hub.Table 3
Path set of experiment 1.
OD: 18–20NumberModePath1c18, 1, 5, 7, 3, 20218, 1, 6, 5, 7, 3, 20318, 1, 5, 7, 8, 3, 204b18, 10, 12, 14, 16, 20518, 10, 12, 9, 14, 16, 206b + c18, 10,1, 5, 7, 3, 20718, 10, 12,5, 7, 3, 20818, 10, 12, 14,7, 3, 209c + b18, 1, 5, 7,14, 16, 201018, 1, 5, 7, 3,16, 201118, 1, 5,12, 14, 16, 2012b + c + b18, 10,1, 5, 7, 14, 16, 201318, 10,1, 5, 7, 3, 16, 201418, 10, 12,5, 7, 14, 16, 2015c + b + c18, 1,10, 12, 5, 7, 3, 201618, 1, 5,12, 14, 7, 3, 201718, 1,10, 12, 14, 7, 3, 20OD: 18–2118c18, 1, 5, 7, 4, 211918, 1, 5, 7, 8, 4, 212018, 1, 6, 5, 7, 4, 2121b18, 10, 12, 9, 15, 17, 2122b + c18, 10,1, 5, 7, 4, 212318, 10, 12,5, 7, 4, 212418, 10, 12, 14,7, 4, 2125c + b18, 1, 5, 7, 4,17, 212618, 1, 5, 7, 8,15, 17, 212718, 1, 5, 7, 8, 4,17, 2128b + c + b18, 10,1, 5, 7, 4, 17, 212918, 10, 12,5, 7, 4, 17, 213018, 10, 12, 14,7, 4, 17, 2131c + b + c18, 1,10, 12, 5, 7, 4, 213218, 1, 5,12, 14, 7, 4, 213318, 1,10, 12, 14, 7, 4, 21OD: 19-2034c19, 2, 5, 7, 3, 203519, 2, 5, 7, 8, 3, 203619, 2, 6, 5, 7, 3, 2037b19, 11, 13, 9, 14, 16, 2038b + c19, 11,2, 5, 7, 3, 203919, 11,2, 5, 7, 8, 3, 204019, 11, 13,6, 5, 7, 3, 2041c + b19, 2, 5, 7,14, 16, 204219, 2, 5, 7, 3,16, 204319, 2, 5,12, 14, 16, 2044b + c + b19, 11,2, 5, 7, 14, 16, 204519, 11,2, 5, 7, 3, 16, 204619, 11,2, 5, 12, 14, 16, 2047c + b + c19, 2, 5,12, 14, 7, 3, 204819, 2, 5, 7,14, 16, 3, 204919, 2, 5,12, 14, 16, 3, 20OD: 19–2150c19, 2, 5, 7, 4, 215119, 2, 5, 7, 8, 4, 215219, 2, 6, 5, 7, 4, 2153b19, 11, 13, 15, 17, 215419, 11, 13, 9, 15, 17, 2155b + c19, 11,2, 5, 7, 4, 215619, 11,2, 5, 7, 8, 4, 215719, 11, 13,6, 5, 7, 4, 2158c + b19, 2, 5, 7, 4,17, 215919, 2, 5, 7, 8,15, 17, 216019, 2, 5, 7, 8, 4,17, 2161b + c + b19, 11,2, 5, 7, 4, 17, 216219, 11,2, 5, 7, 8, 15, 17, 216319, 11,2, 5, 7, 8, 4, 17, 2164c + b + c19, 2, 5,12, 14, 7, 4, 216519, 2, 5,12, 14, 7, 8, 4, 216619, 2, 5, 7,8, 15, 17, 4, 21The optimal flow,F∗, and optimal objective, ξ∗, makes the NLP have the optimal objective, and equilibrium path cost when K = 2 is based on the fmincon solver, and the hybrid MSA are listed in Tables 4 and 5.Table 4
Results of the fmincon solver.
OD: 18–20OD: 18–21OD: 19-20OD: 19–21No.FlowCostNo.FlowCostNo.FlowCostNo.FlowCost1059.618064.7342.652.3500.157.421.358.419064.735059.3513.057.43066.7208.163.536061.352066.44081.7212.564.8378.750.95315.755.5527.455.0222.364.8382.252.354060.760.159.623086.839059.3550.757.47081.724085.240085.0561.957.48080.1250.164.741053.957090.19061.226071.7422.252.3580.457.410059.6272.264.7430.152.359064.311059.6282.864.744053.8602.157.412061.229086.8452.652.3611.257.4130.259.630085.1460.152.362064.314083.331086.84711.550.6631.957.415081.73212.063.148053.9643.355.7161.058.033085.249052.3659.755.717080.166064.3ξ3.750911.0166−4.47420.6910Objective function8.6484e − 6CPU time6.925 sTable 5
Results of the hybrid MSA.
OD: 18–20OD: 18–21OD: 19-20OD: 19–21No.FlowCostNo.FlowCostNo.FlowCostNo.FlowCost1059.418064.6343.052.2500.857.321.358.419064.635059.3513.057.33066.5208.463.536061.352066.44081.8212.664.8378.851.05315.855.7527.354.8221.964.6381.952.254060.960.159.523086.739059.3550.457.37081.624085.340085.1561.657.38080.1250.464.641053.857090.29061.126071.6421.852.2580.457.3100.159.5271.564.643052.459064.411059.6283.264.644053.8601.657.312061.129086.7453.152.2610.857.3130.359.430085.346052.462064.414083.231086.74711.350.8632.857.315081.63212.063.148053.8642.655.9161.058.033085.349052.46510.355.917080.166064.4ξ3.745410.7520−4.74020.4061Objective function3.5948CPU time0.145 sIn Table4, the optimal objectives based on the fmincon solver is 8.6484e − 6, which is close to the coverage criterion zero of Proposition 2. Meanwhile, for all the paths with the same mode, such as paths 18, 19, and 20 with mode c, the paths carrying flow all have the lowest path cost, so the UE principle is satisfied.The other two conditions in section 3.3.1 are also verified. Taking OD 18–20 as an example, the minimum utility values of 0, 1, and 2 transfers are 55.0, 59.6, and 58.0, respectively, and the equilibrium flows are 28.7, 0.1, and 1.2, which obeys condition 1. Among the two transfer trips, the minimum utility values are 59.6 and 58 for modesb + c + b and c + b + c, respectively, and the total flow of 1.2 is divided into 0.2 and 1.0, which also obeys condition 2.Table5 shows that the objective function is 3.5948, and the UE principle is also obeyed. Compared with Table 4, the flow patterns are similar, while the CPU times are 0.145 s and 6.925 s. These results show that the hybrid MSA could obtain an approximate result with a much faster speed.


## 4.2. Experiment 2

### 4.2.1. Network and Parameters

Experiment 2 tested the performance of the HGSADC, as shown in Figure6. The three subnetworks from left to right are the road, metro, and bus subnetworks. The link and hub free-flow travel times and capacities are denoted by blue and green texts.Figure 6
Network of experiment 2.The following modifications were made based on the network of experiment 1:(1)
A new parking lot was added in node 26 and connected to road node 5 with bidirectional auxiliary links 37 and 51. Users must transfer from a parking lot to a metro station to execute a P&R combined trip.(2)
Two parallel metro lines with a frequency of 0.083 h/veh were added to provide more mode choices. Only one bus route 21–25 was kept at 0.068 h/veh to leave space for optimization.(3)
In experiment 1, the costs of transfer links, embarking/alighting links, and auxiliary links were default to be 0. In this experiment, to verify the results after the optimization of the hubs, the costs of transfer links, embarking/alighting links, and auxiliary links were added with the sum of the walking time of the links, walking time inside the hubs (BPR function with green text as the parameters), waiting time (the route frequency), and pure penalty (a fixed time of 0.083 h).(4)
A Euclidean space coordinate system was built, where thex and y axes represented the position of the node (see Figure 7). Because the four OD pairs were connected to several direct car, bus, or metro services to avoid the situation where no users take combined trips, one unit of length in the network was set at 10 km, and the car, metro, and bus speeds were set at 30, 25, and 20 km/h, respectively.Experiment 2 aimed to optimize the parking lots and bus routes to integrate the network. The maximum number of bus routes was set to two, while the number of parking lots was set to 1–6. The original number of buses was five. The weights of the four objectives such as the total time, construction cost, modal shift value, and service level value, were set at 1, 1, 2, and 3, respectively.Figure 7
Coordinates of experiment 2 subnetworks.
(a)(b)(c)(d)


### 4.2.2. Results of Experiment 2

Four sets of GAs tested the model performance: GA1 was the normal GA, GA2 was a GA with a diversification mechanism alone, GA3 was a GA with a biased fitness function alone, and GA4 was the HGSADC. Each GA was ended when 500 generations were produced and repeated for 20 runs. The average results of the original network and 4 GAs are listed in Table6.Table 6
Average results of 4 GAs.
GA1GA2GA3HGSADCWeighted objective21.105817.364117.364117.3641Travel time292.3744289.5535289.5535289.5535Construction cost33016.7230084.7230084.7230084.72Modal shift value−19.1918−28.6002−28.6002−28.6002Service level value451.1113405.3683405.3683405.3683Convergence iterations59.7111.392.447.2CPU time (s)97.1122.5103.4193.4In all 20 runs, GA2, GA3, and HGSADC all converged to the same results with a weighted objective value of 17.3641. However, GA1 was trapped in a local optimal objective with an average objective value of 21.1058, showing that both the diversification mechanism and biased fitness function contributed to improving the exploration capacity. Meanwhile, HGSADC had far fewer convergence iterations than the other three GAs. Because the diversification mechanism involved a reinitialization process to obtain more genetic materials, it also made the average CPU time longer than the other three GAs. Considering it is the total time of 500 iterations, the fast convergence speed can be useful in large-scale problems.


### 4.2.3. Analyzation of Schemes

The performances of the original, GA1, and HGSADC schemes were then compared. See Table7.Table 7
Results of three schemes.
SchemeParking lotBus routesNumber of busesBus frequencyTotal timeConstruction costModal shiftService levelOriginal521-22-23-24-2550.0678296.2606——478.1152GA15, 626-27-28–2920.2494288.222925187.08−17.5560443.090621-22-23-24-2530.113HGSADC526-27-28-29-3030.1130282.180030084.72−28.6002405.368321-22-23-24-2520.1694(1) Original Network. The flow distribution of the original network is shown in Figure 8. From left to right, the multimodal network comprises a green road subnetwork, brown metro subnetwork, and a purple bus subnetwork. The link flow is noted beside the link in black text, while the hub flow is shown in blue text. Meanwhile, transfer flow is noted beside the transfer link. The pink text shows the mtr + car and mtr + bus transfer flow, while the orange text shows the car + mtr and bus + mtr transfer flow.Figure 8
Flow distribution of the original network.The node flow is the sum of all the links entering the node. The two parallel OD pairs, 9–11 and 10–12, can be served by the car and two parallel metro lines directly. The speed of the car is faster than that of the metro, while the capacity of the road link is 10, which is not sufficient to meet the 30 demands of the parallel OD pairs. This portion of the users will choose the single-car mode first and then the metro. Meanwhile, because a parking lot is built at node 26, and road link 1-5-7-3 is congested (the flow/capacity ratio is higher than 1), while metro line 13-14-15-16 is underutilized (lower than 0.75), some users will transfer from parking lot 26 to metro station 14 to bypass the congested road link, 5-7-3.For the two diagonal OD pairs, OD pair 9–12 can be served by the bus route directly. Among the 30 demands, 10.2 users choose bus link 21-22, and the other 8.9 flow of this route is themtr + bus flow. The bus route relieves the congestion of the road links in this direction, and road links 1–6, 5-6, 7-8, and 7-4 all have a good service level. However, the other diagonal OD pair, 10-11, can only be served by the road subnetwork directly, which makes the service level of related links 2–5 and 8-3 very low.(2) GA1 Scheme. The GA1 scheme adds a new parking lot in node 31, as well as bus route 26-27-2829 with frequency 0.12 h/veh. The original parking lot 30 and bus route 21-22-23-24-25 are kept, but the frequency of the original bus route increases to 0.112 h/veh. The flow distribution is illustrated in Figure 9.Figure 9
Flow distribution of GA1 scheme.The newly added bus route, 26–29, only attracts 3.7 users because it is parallel to metro route 13–16. However, the congestion of road links 1-5-7-3 is still slightly reduced. Meanwhile, the new parking lot at 31 attracts 14.6 users to use the P&R mode, which reduces the congestion of road link 6-8-4.(3) HGSADC Scheme. The HGSADC scheme adds a symmetrical bus route, 26-27-28-29-30, with a frequency of 0.113 h/veh. The parking lot of the original network is kept. See Figure 10.Figure 10
Flow distribution of HGSADC scheme.The added bus route provides a direct bus service for OD pair 10-11, which attracts eight single bus mode users to bus stop 26. This bus route creates two new mode choices for this OD pair,mtr + bus and bus + mtr, and except for 8 direct users, 8.1 users transfer from metro station 18 to bus stop 27 and 2.2 users transfer from bus stop 29 to metro station 13, which significantly relieves the congestion of road links 2–5 and 8-3. With more users choosing PT, the flows of other single-car mode paths, such as 2-5-7-3 and 2-6-8-3, also decrease by a certain amount, which leads to an overall reduction in congestion of road subnetwork. Because the weight of the service level value is the highest, the lowest service level value of 405.4 makes the weighted objective of the HGSADC lower than that of the GA1 scheme even though the total time is similar, and the construction cost is higher.


## 4.2.1. Network and Parameters

Experiment 2 tested the performance of the HGSADC, as shown in Figure6. The three subnetworks from left to right are the road, metro, and bus subnetworks. The link and hub free-flow travel times and capacities are denoted by blue and green texts.Figure 6
Network of experiment 2.The following modifications were made based on the network of experiment 1:(1)
A new parking lot was added in node 26 and connected to road node 5 with bidirectional auxiliary links 37 and 51. Users must transfer from a parking lot to a metro station to execute a P&R combined trip.(2)
Two parallel metro lines with a frequency of 0.083 h/veh were added to provide more mode choices. Only one bus route 21–25 was kept at 0.068 h/veh to leave space for optimization.(3)
In experiment 1, the costs of transfer links, embarking/alighting links, and auxiliary links were default to be 0. In this experiment, to verify the results after the optimization of the hubs, the costs of transfer links, embarking/alighting links, and auxiliary links were added with the sum of the walking time of the links, walking time inside the hubs (BPR function with green text as the parameters), waiting time (the route frequency), and pure penalty (a fixed time of 0.083 h).(4)
A Euclidean space coordinate system was built, where thex and y axes represented the position of the node (see Figure 7). Because the four OD pairs were connected to several direct car, bus, or metro services to avoid the situation where no users take combined trips, one unit of length in the network was set at 10 km, and the car, metro, and bus speeds were set at 30, 25, and 20 km/h, respectively.Experiment 2 aimed to optimize the parking lots and bus routes to integrate the network. The maximum number of bus routes was set to two, while the number of parking lots was set to 1–6. The original number of buses was five. The weights of the four objectives such as the total time, construction cost, modal shift value, and service level value, were set at 1, 1, 2, and 3, respectively.Figure 7
Coordinates of experiment 2 subnetworks.
(a)(b)(c)(d)


## 4.2.2. Results of Experiment 2

Four sets of GAs tested the model performance: GA1 was the normal GA, GA2 was a GA with a diversification mechanism alone, GA3 was a GA with a biased fitness function alone, and GA4 was the HGSADC. Each GA was ended when 500 generations were produced and repeated for 20 runs. The average results of the original network and 4 GAs are listed in Table6.Table 6
Average results of 4 GAs.
GA1GA2GA3HGSADCWeighted objective21.105817.364117.364117.3641Travel time292.3744289.5535289.5535289.5535Construction cost33016.7230084.7230084.7230084.72Modal shift value−19.1918−28.6002−28.6002−28.6002Service level value451.1113405.3683405.3683405.3683Convergence iterations59.7111.392.447.2CPU time (s)97.1122.5103.4193.4In all 20 runs, GA2, GA3, and HGSADC all converged to the same results with a weighted objective value of 17.3641. However, GA1 was trapped in a local optimal objective with an average objective value of 21.1058, showing that both the diversification mechanism and biased fitness function contributed to improving the exploration capacity. Meanwhile, HGSADC had far fewer convergence iterations than the other three GAs. Because the diversification mechanism involved a reinitialization process to obtain more genetic materials, it also made the average CPU time longer than the other three GAs. Considering it is the total time of 500 iterations, the fast convergence speed can be useful in large-scale problems.


## 4.2.3. Analyzation of Schemes

The performances of the original, GA1, and HGSADC schemes were then compared. See Table7.Table 7
Results of three schemes.
SchemeParking lotBus routesNumber of busesBus frequencyTotal timeConstruction costModal shiftService levelOriginal521-22-23-24-2550.0678296.2606——478.1152GA15, 626-27-28–2920.2494288.222925187.08−17.5560443.090621-22-23-24-2530.113HGSADC526-27-28-29-3030.1130282.180030084.72−28.6002405.368321-22-23-24-2520.1694(1) Original Network. The flow distribution of the original network is shown in Figure 8. From left to right, the multimodal network comprises a green road subnetwork, brown metro subnetwork, and a purple bus subnetwork. The link flow is noted beside the link in black text, while the hub flow is shown in blue text. Meanwhile, transfer flow is noted beside the transfer link. The pink text shows the mtr + car and mtr + bus transfer flow, while the orange text shows the car + mtr and bus + mtr transfer flow.Figure 8
Flow distribution of the original network.The node flow is the sum of all the links entering the node. The two parallel OD pairs, 9–11 and 10–12, can be served by the car and two parallel metro lines directly. The speed of the car is faster than that of the metro, while the capacity of the road link is 10, which is not sufficient to meet the 30 demands of the parallel OD pairs. This portion of the users will choose the single-car mode first and then the metro. Meanwhile, because a parking lot is built at node 26, and road link 1-5-7-3 is congested (the flow/capacity ratio is higher than 1), while metro line 13-14-15-16 is underutilized (lower than 0.75), some users will transfer from parking lot 26 to metro station 14 to bypass the congested road link, 5-7-3.For the two diagonal OD pairs, OD pair 9–12 can be served by the bus route directly. Among the 30 demands, 10.2 users choose bus link 21-22, and the other 8.9 flow of this route is themtr + bus flow. The bus route relieves the congestion of the road links in this direction, and road links 1–6, 5-6, 7-8, and 7-4 all have a good service level. However, the other diagonal OD pair, 10-11, can only be served by the road subnetwork directly, which makes the service level of related links 2–5 and 8-3 very low.(2) GA1 Scheme. The GA1 scheme adds a new parking lot in node 31, as well as bus route 26-27-2829 with frequency 0.12 h/veh. The original parking lot 30 and bus route 21-22-23-24-25 are kept, but the frequency of the original bus route increases to 0.112 h/veh. The flow distribution is illustrated in Figure 9.Figure 9
Flow distribution of GA1 scheme.The newly added bus route, 26–29, only attracts 3.7 users because it is parallel to metro route 13–16. However, the congestion of road links 1-5-7-3 is still slightly reduced. Meanwhile, the new parking lot at 31 attracts 14.6 users to use the P&R mode, which reduces the congestion of road link 6-8-4.(3) HGSADC Scheme. The HGSADC scheme adds a symmetrical bus route, 26-27-28-29-30, with a frequency of 0.113 h/veh. The parking lot of the original network is kept. See Figure 10.Figure 10
Flow distribution of HGSADC scheme.The added bus route provides a direct bus service for OD pair 10-11, which attracts eight single bus mode users to bus stop 26. This bus route creates two new mode choices for this OD pair,mtr + bus and bus + mtr, and except for 8 direct users, 8.1 users transfer from metro station 18 to bus stop 27 and 2.2 users transfer from bus stop 29 to metro station 13, which significantly relieves the congestion of road links 2–5 and 8-3. With more users choosing PT, the flows of other single-car mode paths, such as 2-5-7-3 and 2-6-8-3, also decrease by a certain amount, which leads to an overall reduction in congestion of road subnetwork. Because the weight of the service level value is the highest, the lowest service level value of 405.4 makes the weighted objective of the HGSADC lower than that of the GA1 scheme even though the total time is similar, and the construction cost is higher.


## 4.3. Experiment 3

### 4.3.1. Original Network

The model was then tested in a real-size network of Jianye District, Nanjing. See Figure11 (Open street map, 2022). The research area covers 80.87 square kilometers, with a population of 489,800. There are two metro lines running through this district and their intersection node, and Yuantong station lies in the central business district (CBD) of the research area. The yellow isolated node in the upper left corner is Linjiang Station. Since Nanjing is divided into two parts by the Yangtze River, Linjiang Station is an important node which connects the south and north banks of the city, and a typical P&R parking lot is built there.Figure 11
Research area.The network comprises 5 subnetworks, i.e., bicycle, car, car hailing, metro, and bus subnetwork. In Figure12, from (a) to (e), the blue, green, and red subnetworks are bicycle, car, and car hailing subnetworks, and they share the same physical network with 125 nodes and 428 links; meanwhile, car and car hailing subnetwork share the same flow. The blue bicycle subnetwork contains 55 bicycle stations while the green car subnetwork has 45 parking lots. The brown metro subnetwork has 2 metro lines and 18 metro stations. The purple bus subnetwork consists of 11 bus lines and 95 bus stops. The district is divided into 22 traffic zones with 462 OD pairs, see Figure 12(f).Figure 12
Original subnetworks. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork. (f) OD pair.
(a)(b)(c)(d)(e)(f)


### 4.3.2. Results of Experiment 3

The decision variables were bicycle stations, parking lots, bus layouts, and their frequencies. The iteration number was set as 3000, and the maximum transfer numberK was 2. The number of bus routes was kept at 11. The number of bicycle stations, parking lots, and bus stops was set from 50 to 200, 45 to 100, and 80 to 200, respectively. The bus stop distance constraint was from 0.3 km to 2 km and 0.0833 h/veh to 0.5 h/veh for bus route frequency constraint. The bus route length constraint was from 5 km to 25 km. 40 paths for each OD pair were chosen to ensure the stability of hybrid MSA in obtaining the equilibrium flow. The pure penalty is 5 mins.(1) Travel Time and Construction Cost. The algorithm converges in 2335 iterations with 172417.54 s. GMS/TAP function accounts for 77% of the total time, in which path generation, path utility, and hybrid MSA takes 39.5%, 4.8%, and 32.7%. The HGSADC and other functions take 23% of the total time. The best solution contains 51 parking lots, 74 bicycle stations, and 11 bus routes with 117 bus stops. See Figure 13. With a construction cost of 11566314.32 RMB, the travel time decreases from 749042.8 h to 680345.5 h, which improves 9.2%. Since this study concentrates on the modal shift and balanced use of the whole network, comparing with Szeto et al. [69], in which a 22.7% reduction of travel time in transit subnetwork is obtained, the improvement of travel time is not obvious in this experiment.Figure 13
Subnetworks of the best solution. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Bus subnetwork.
(a)(b)(c)(2) Service Level Value. The service level of the original network is presented in Figure 14, and in the original network, the road and metro subnetworks are overloaded, especially in the CBD area with longitude from 118.7 to 118.72, latitude from 31.98 to 31.99. Meanwhile, the flow between north bank and south bank also concentrates in the car and metro mode, and the service level in the corresponding links and hubs are all lower than level E. With a congested metro and car subnetwork, the bicycle and bus subnetworks in Figures 14(a) and 14(e) are underutilized; to fully use the potential capacity of these subnetworks, a modal shift is necessary.Figure 14
Service level of original network. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork.
(a)(b)(c)(d)(e)In the optimal scheme (see Figure15), the service level value decreases from 40253.82 to 29903.51, which improves 25.71%. Although the main traffic corridor from west to east is still congested in road subnetwork, the number of over congested links and hubs (service level higher than 1) decreases with 21 and 3 for road and metro subnetworks, respectively. These flows transfer to bicycle and bus subnetworks, and an increase in usage with 36 and 44 links and hubs is obtained.Figure 15
Service level of best solution. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork.
(a)(b)(c)(d)(e)(3) Modal Shift Value. The modal shift value is −131305.2. After deleting the modes with total flow lower than 2000, the mode flow of single, one transfer, and two transfers modes are shown in Figures 16(a)–16(c), respectively. In the x axis, the mode bic, car, ch, mtr, and bus are denoted as b, c, h, m, and B, respectively. The combined trip, such as car + mtr is denoted as cm. While on the y axis, the mode flow is classified per kilometers according to the direct OD distance. The z axis is the amount of mode flow.Figure 16
Mode flow of original network. (a) Single-mode flow. (b) One transfer mode flow. (c) Two transfers mode flow.
(a)(b)(c)In the original network, 77.33% of the flow is single-mode within which 37.17% comes from the private car mode. Since there are only two metro lines while a complete bus network is built, the metro subnetwork is mostly used in combined trips. The bicycle is normally used alone in short distance trips and in long-distance trips to complement bus and metro subnetworks, and the car hailing subnetwork plays a similar role as the private car but with a lower flow.Considering the original road subnetwork is over congested, the single-modec is encouraged to be reduced with positive weight value. Meanwhile, one transfer modes in short distance and 2 transfers modes in short and middle distance needs to be cut down as well. The modal shift weight matrix is divided into seven columns and three rows, see Table 8.Table 8
The modal shift weight matrix.
bcmhBOne transfer modeTwo transfers modeShort distance (0∼4 km)−33−1−1−212Middle distance (5∼8 km)−22−1−1−3−11Long distance (9∼11 km)−11−1−1−3−2−2The flow of each mode in the optimal scheme is shown in Figure17. The 0 transfer modes and two transfers modes users decrease with 6885 and 12797, respectively, and these users transfer to one transfer mode. Among single-mode trips, a decrease of 27768 car users and an increase of 4853, 14860, and 20853 bic, mtr, and bus users are observed. These four single-modes produce a weighted modal shift value −163681. In one transfer trips, the modal shift does not change as wishes, the increase in 19682 one transfer modes users concentrate in short distance OD pairs, especially the bm and mb modes. But in other perspective, this means the new bicycle station scheme improves the connection between bicycle and metro subnetworks. For two transfers modes, the low flow does not have an obvious influence on the modal shift value.Figure 17
Mode flow of the solution. (a) Single-mode flow. (b) One transfer mode flow. (c) Two transfers mode flow.
(a)(b)(c)


### 4.3.3. Sensitive Tests of the Parameters

(1) Effect of the Transfer Penalty. This section tested the effect of transfer penalty on the flow distribution of the optimized network. Since there are several types of transfer costs which are related to the connected modes, it is difficult to test these parameters one by one. Considering every transfer involves in a pure penalty with 5 mins, pure penalty represents the natural resistance of the users to the transfer. Thus, the experiment tested the network with a pure penalty from 0 to 15 min, with an interval of 0.25 min. The total mode miles are chosen as the index, which is used to describe the overall usage of the corresponding subnetwork. See Figure 18.Figure 18
Effect of pure penalty on the flow distribution. (a) Effect of pure penalty on the subnetworks. (b) Effect of pure penalty on the number of transfer.
(a)(b)In Figure18(a), when the transfer penalty increases, the miles on car and bicycle subnetwork increases. For the metro and bus subnetworks, which rely more on transfer, the increase of transfer penalty leads to a decrease of their travel miles. In Figure 18(b), we can see that the number of transfers is highly related to the transfer time. Intuitively, the decrease of the number of transfers is good for the transportation system. However, in this test, this is achieved by increasing the transfer time rather than improving the connection of the network. The long-time of waiting pushes the users to the congested links of the single modes, while they cannot use the underutilized links of other subnetworks through combined trip to bypass the congested links. The results suggest that the number of transfers alone may not be a good index to value the integration of the multimodal network.(2) Effect of the Congestion. We further investigated the congestion effect under different demand levels. See Figure 19. The OD demand varies from 5% to 300% of the original demand with an increase of 5%. As the OD demand increases, the users choose road and metro subnetworks at first. When the demand reaches 45% of the original demand, the car subnetwork becomes relatively more congested than other three subnetworks, i.e., bus, bicycle, and car hailing, and their share of travel miles increases while car and metro decreases. As shown in Figure 15, under the original demand, the car and metro subnetworks are already over congested, the decrease speed for these two modes is increasing, and more users choose other three modes.Figure 19
Effect of congestion on the flow distribution.This result suggests that there may be one or two modes that the users will choose at first because their competitions are much better than other modes. Only when these preferred subnetworks are congested, then users will choose the remaining modes. The natural differences of the travel modes are hard to change only via an improvement of the infrastructures. A systematic demand control and operational management are also required to balance the flow distribution.


## 4.3.1. Original Network

The model was then tested in a real-size network of Jianye District, Nanjing. See Figure11 (Open street map, 2022). The research area covers 80.87 square kilometers, with a population of 489,800. There are two metro lines running through this district and their intersection node, and Yuantong station lies in the central business district (CBD) of the research area. The yellow isolated node in the upper left corner is Linjiang Station. Since Nanjing is divided into two parts by the Yangtze River, Linjiang Station is an important node which connects the south and north banks of the city, and a typical P&R parking lot is built there.Figure 11
Research area.The network comprises 5 subnetworks, i.e., bicycle, car, car hailing, metro, and bus subnetwork. In Figure12, from (a) to (e), the blue, green, and red subnetworks are bicycle, car, and car hailing subnetworks, and they share the same physical network with 125 nodes and 428 links; meanwhile, car and car hailing subnetwork share the same flow. The blue bicycle subnetwork contains 55 bicycle stations while the green car subnetwork has 45 parking lots. The brown metro subnetwork has 2 metro lines and 18 metro stations. The purple bus subnetwork consists of 11 bus lines and 95 bus stops. The district is divided into 22 traffic zones with 462 OD pairs, see Figure 12(f).Figure 12
Original subnetworks. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork. (f) OD pair.
(a)(b)(c)(d)(e)(f)


## 4.3.2. Results of Experiment 3

The decision variables were bicycle stations, parking lots, bus layouts, and their frequencies. The iteration number was set as 3000, and the maximum transfer numberK was 2. The number of bus routes was kept at 11. The number of bicycle stations, parking lots, and bus stops was set from 50 to 200, 45 to 100, and 80 to 200, respectively. The bus stop distance constraint was from 0.3 km to 2 km and 0.0833 h/veh to 0.5 h/veh for bus route frequency constraint. The bus route length constraint was from 5 km to 25 km. 40 paths for each OD pair were chosen to ensure the stability of hybrid MSA in obtaining the equilibrium flow. The pure penalty is 5 mins.(1) Travel Time and Construction Cost. The algorithm converges in 2335 iterations with 172417.54 s. GMS/TAP function accounts for 77% of the total time, in which path generation, path utility, and hybrid MSA takes 39.5%, 4.8%, and 32.7%. The HGSADC and other functions take 23% of the total time. The best solution contains 51 parking lots, 74 bicycle stations, and 11 bus routes with 117 bus stops. See Figure 13. With a construction cost of 11566314.32 RMB, the travel time decreases from 749042.8 h to 680345.5 h, which improves 9.2%. Since this study concentrates on the modal shift and balanced use of the whole network, comparing with Szeto et al. [69], in which a 22.7% reduction of travel time in transit subnetwork is obtained, the improvement of travel time is not obvious in this experiment.Figure 13
Subnetworks of the best solution. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Bus subnetwork.
(a)(b)(c)(2) Service Level Value. The service level of the original network is presented in Figure 14, and in the original network, the road and metro subnetworks are overloaded, especially in the CBD area with longitude from 118.7 to 118.72, latitude from 31.98 to 31.99. Meanwhile, the flow between north bank and south bank also concentrates in the car and metro mode, and the service level in the corresponding links and hubs are all lower than level E. With a congested metro and car subnetwork, the bicycle and bus subnetworks in Figures 14(a) and 14(e) are underutilized; to fully use the potential capacity of these subnetworks, a modal shift is necessary.Figure 14
Service level of original network. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork.
(a)(b)(c)(d)(e)In the optimal scheme (see Figure15), the service level value decreases from 40253.82 to 29903.51, which improves 25.71%. Although the main traffic corridor from west to east is still congested in road subnetwork, the number of over congested links and hubs (service level higher than 1) decreases with 21 and 3 for road and metro subnetworks, respectively. These flows transfer to bicycle and bus subnetworks, and an increase in usage with 36 and 44 links and hubs is obtained.Figure 15
Service level of best solution. (a) Bicycle subnetwork. (b) Car subnetwork. (c) Car hailing subnetwork. (d) Metro subnetwork. (e) Bus subnetwork.
(a)(b)(c)(d)(e)(3) Modal Shift Value. The modal shift value is −131305.2. After deleting the modes with total flow lower than 2000, the mode flow of single, one transfer, and two transfers modes are shown in Figures 16(a)–16(c), respectively. In the x axis, the mode bic, car, ch, mtr, and bus are denoted as b, c, h, m, and B, respectively. The combined trip, such as car + mtr is denoted as cm. While on the y axis, the mode flow is classified per kilometers according to the direct OD distance. The z axis is the amount of mode flow.Figure 16
Mode flow of original network. (a) Single-mode flow. (b) One transfer mode flow. (c) Two transfers mode flow.
(a)(b)(c)In the original network, 77.33% of the flow is single-mode within which 37.17% comes from the private car mode. Since there are only two metro lines while a complete bus network is built, the metro subnetwork is mostly used in combined trips. The bicycle is normally used alone in short distance trips and in long-distance trips to complement bus and metro subnetworks, and the car hailing subnetwork plays a similar role as the private car but with a lower flow.Considering the original road subnetwork is over congested, the single-modec is encouraged to be reduced with positive weight value. Meanwhile, one transfer modes in short distance and 2 transfers modes in short and middle distance needs to be cut down as well. The modal shift weight matrix is divided into seven columns and three rows, see Table 8.Table 8
The modal shift weight matrix.
bcmhBOne transfer modeTwo transfers modeShort distance (0∼4 km)−33−1−1−212Middle distance (5∼8 km)−22−1−1−3−11Long distance (9∼11 km)−11−1−1−3−2−2The flow of each mode in the optimal scheme is shown in Figure17. The 0 transfer modes and two transfers modes users decrease with 6885 and 12797, respectively, and these users transfer to one transfer mode. Among single-mode trips, a decrease of 27768 car users and an increase of 4853, 14860, and 20853 bic, mtr, and bus users are observed. These four single-modes produce a weighted modal shift value −163681. In one transfer trips, the modal shift does not change as wishes, the increase in 19682 one transfer modes users concentrate in short distance OD pairs, especially the bm and mb modes. But in other perspective, this means the new bicycle station scheme improves the connection between bicycle and metro subnetworks. For two transfers modes, the low flow does not have an obvious influence on the modal shift value.Figure 17
Mode flow of the solution. (a) Single-mode flow. (b) One transfer mode flow. (c) Two transfers mode flow.
(a)(b)(c)


## 4.3.3. Sensitive Tests of the Parameters

(1) Effect of the Transfer Penalty. This section tested the effect of transfer penalty on the flow distribution of the optimized network. Since there are several types of transfer costs which are related to the connected modes, it is difficult to test these parameters one by one. Considering every transfer involves in a pure penalty with 5 mins, pure penalty represents the natural resistance of the users to the transfer. Thus, the experiment tested the network with a pure penalty from 0 to 15 min, with an interval of 0.25 min. The total mode miles are chosen as the index, which is used to describe the overall usage of the corresponding subnetwork. See Figure 18.Figure 18
Effect of pure penalty on the flow distribution. (a) Effect of pure penalty on the subnetworks. (b) Effect of pure penalty on the number of transfer.
(a)(b)In Figure18(a), when the transfer penalty increases, the miles on car and bicycle subnetwork increases. For the metro and bus subnetworks, which rely more on transfer, the increase of transfer penalty leads to a decrease of their travel miles. In Figure 18(b), we can see that the number of transfers is highly related to the transfer time. Intuitively, the decrease of the number of transfers is good for the transportation system. However, in this test, this is achieved by increasing the transfer time rather than improving the connection of the network. The long-time of waiting pushes the users to the congested links of the single modes, while they cannot use the underutilized links of other subnetworks through combined trip to bypass the congested links. The results suggest that the number of transfers alone may not be a good index to value the integration of the multimodal network.(2) Effect of the Congestion. We further investigated the congestion effect under different demand levels. See Figure 19. The OD demand varies from 5% to 300% of the original demand with an increase of 5%. As the OD demand increases, the users choose road and metro subnetworks at first. When the demand reaches 45% of the original demand, the car subnetwork becomes relatively more congested than other three subnetworks, i.e., bus, bicycle, and car hailing, and their share of travel miles increases while car and metro decreases. As shown in Figure 15, under the original demand, the car and metro subnetworks are already over congested, the decrease speed for these two modes is increasing, and more users choose other three modes.Figure 19
Effect of congestion on the flow distribution.This result suggests that there may be one or two modes that the users will choose at first because their competitions are much better than other modes. Only when these preferred subnetworks are congested, then users will choose the remaining modes. The natural differences of the travel modes are hard to change only via an improvement of the infrastructures. A systematic demand control and operational management are also required to balance the flow distribution.


## 5. Conclusions

This paper proposed a bilevel IMNDP to optimize subnetwork hubs, guide a modal shift, and achieve a balanced flow distribution among all the subnetworks. The upper-level problem had the goal of optimizing various network hubs, and the decision variables were classified as the route-based and location-based optimization. The route-based optimization included the route layout and frequency setting, whereas the location-based consisted of the hub locations. In the lower-level problem, a NL mode-choice structure was devised to describe the mode choice of the contained modes in the path set. A unified VI formulation was developed and transformed into an NLP that could be solved using commercial software. To facilitate the application in a real-size network, a hybrid MSA was further developed to solve the lower-level problem. To solve the proposed bilevel model, an HGSADC along with the initialization, crossover, mutation, repair operators, biased fitness, and diversification mechanisms were introduced in this paper. The solution properties and model performance were verified in three experiments, and a significant modal shift and service level increase could be observed.This study opens several research directions. (1) In this study, the modal shift was valued with a predefined matrix, which was based on the current use of the network and experience, while the balanced use of the network was assessed by the service levels of all the links and hubs. However, how to assess a real balanced use and what is a reasonable modal shift pattern still requires much work. (2) The optimization of different hubs in this paper only considered limited number of constraints. However, the metro, car hailing, and many other modes were not shown in this paper, and the optimization of these modes and their combinations also requires deeper research. (3) As discussed in the literature review, the integration of multimodal network can be achieved by many other factors, such as policy, environment, market needs, fare, and coordinated schedules which can be used as the decision variables as well. (4) The impedance function requires a systematic data collection and calibration methods regarding to the different links, hubs, single, and combined modes. (5) Besides the mode choice, multiple classes of users, link interactions, dynamic assignments, demand uncertainties, exact-method algorithms and many other factors could also apply to the GMS/TAP.The main contributions of this paper are as follows: (1) the extension of the traditional MS/TAP model framework for a wider range of modes. (2) The methodologies for the assessment and optimal design of a hub integration scheme. However, these contributions are more theoretical than practical. Without a calibrated impedance function, the model is not suitable to be adopted into the real-size network with more travel modes. Thus, the priority of our future research is to calibrate a system of impedance functions. Considering there are mature data for the car, metro, and bus subnetworks separately, the research could start from the extraction method of the transfer flow between the three modes using the cellular signaling data, bus card data, and GPS data.



---

*Source: 3442921-2022-12-13.xml*