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.

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 Section 2. 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 ticketing

The 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 [3739], 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 [5254] 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 [5759]. 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 [6466] 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 [7275]. 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 Table 1.

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 Section 3.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 by . Denote the subnetworks as , where each subnetwork corresponds to a basic travel mode, m, where . , , and are the node, subnetwork link, and hub set, respectively. These subnetworks are combined with three link subsets. The set of auxiliary links, , connects the subnetwork nodes and their hubs within an acceptable walking distance. The sets of embarking and alighting links ( and , respectively) connect the OD nodes with the subnetwork nodes and represent the embarking and alighting activities. The set of transfer links, , 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 Figure 1. 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.

To build a NL choice structure from the path set. A given path set of OD pair is denoted as , where can be divided into K + 1 subsets according to the number of transfers, . For the sth path in , . The mode of can be denoted with a mode chain , where is the nth basic mode, . For convenience, a complete mode chain, is denoted as , and an intermediate mode chain with the first n modes, is denoted as .

For the NL mode-choice structure, denote the qth nest in column k and level n as , which can be simplified to . In this study, if the mode sequence in is the same as the mode sequence in , then . If the first n modes in are the same as , then , meaning is the lower-level nest of .

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 , is fixed. In the mode-choice structure, the flow of nests in column k, level −1 is denoted as . Then, C1 is formulated as , , and , and should satisfy the following constraints:where is the equilibrium utility of the nest in column k, level −1, , and 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 and its lower-level nest flow as . C2 is represented as , k = 0∼K and n = 0∼k − 1, and the following constraints should be obeyed:where is the equilibrium utility of nest , is the equilibrium utility set of all nests in level n + 1, and is a binary variable for determining whether .C3: Choice of Path. Paths in each mode nest should follow the UE principle. Denote the cost of path as . C3 is represented as , , , and :

Equation (4) states that only when path cost equals its corresponding mode nest’s equilibrium utility , then can path carry flow; otherwise, the path flow is 0. It should be noted that the condition is to find the corresponding mode nest of the complete mode chain , 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.

Denote as the flow of path ’s corresponding nest in level −1 and as the flow of ’s corresponding nest in level n, which means and . Then, C1 can be reformulated as , , , and :

While C2 can be reformulated as , k = 0∼K, , n = 0∼k – 1, and :

Equations (5) and (6) only replaces and in equations (1) and (2) with and , and the other variables are the same.

3.3.2. Unified Optimization Problem

Denote the equilibrium flow vector of all OD pairs as . 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 solution 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 , , , and ;whereThe solution space of the VI problem isEquation (9) ensures that the total path flow of OD equals demand , 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, , where indicates that flow vector F is an n-dimensional real number vector, , B is the incident matrix of equation (9), denoting as . The VI problem in equation (7) can be transformed into the following NLP.

Proposition 2. Suppose that D is a nonempty polyhedron, . Then, solves the VI problem in equation (7) if the following optimization problem has an optimal value of 0 and , such that is an optimal solution:

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, 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 Table 2.

IMNDP 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

Since it is beyond the work of this paper to cover all hubs of different modes, for simplification, only the five modes in Section 3.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:

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, is denoted as the total flow of mode m among all OD pairs with distance l. The predefined weight parameter is used to value the 1 unit of flow change compared to the original network if 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:

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:

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.

(1)for route r from 1 to Rmax
(2) Whether route r needs to construct a new terminal with probability and
(3) Determine the starting and ending node: if new terminal should be built, select a node from CN as the starting node or ending node 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, is not the , and the route satisfies the route length constraint, bus node repeat constraint, and time < Tmax do
(5)  Find the node set with distance to satisfying the stop distance constraint
(6)  if is empty or all nodes have been visited in
(7)   , t = t-1
(8)  elseif is not empty
(9)   Randomly select a node from that has not been visited
(10)   if
    
   elseif
    Mark has been visited and go back to line 6
   end
(11)if route r is feasible then
  Insert route r into initialization solution
(12)elseif route r is infeasible then
  Go 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 . It is defined as the average distance from S to its closest neighbours in the subpopulation. 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 is the sum of the above three distances.

A dynamic penalty is then introduced to handle the infeasible solutions. Suppose the number of feasible and infeasible solutions are and , respectively. A basic penalty value, , is adjusted for every iterations, and when the number of increases by 1, decreases by , and vice versa. (see Algorithm 2)

(1)Initialization
(2)while number of iterations without improvements < ItNI do
(3)while maximum solution size is not reached do
  select parent solutions S1 and S2, create offspring O using route-based and location-based crossover operator. Insert O into solution set (crossover)
  end while
(4)  select solution S1 and mutate with route-based and location-based mutation operator (mutation)
(5) for each solution S1, if infeasible then
(6) Repair S1 with probability Prep (repair)
  if repairable then
   insert into feasible subpopulation.
  else
   insert into infeasible subpopulation
  end if
(7)elseif S1 feasible then
  insert into feasible subpopulation
end if
(8) calculate objective values with GMS/TAP
(9) calculate the biased fitness value
  Select survivors(selection)
(10)if best solution not improved for Itdiv iterations, then
  Diversify population
(11) Adjust penalty parameters for infeasibility
end while
 Return best feasible solution

Then, the survivors of each iteration are selected based on the rank of objective function and the rank of Hamming distance , 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, , is formulated:where 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.

To adopt the original network into GMS/TAP, a modified network with the same properties was proposed. As shown in Figure 5, 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.

For simplification, the car 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 number K = 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 , , , , , , , and . The path set for is presented in Table 3. A bold number in the path column represents the transfer hub.

The optimal flow, , 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.

In Table 4, 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 modes b+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.

Table 5 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 Figure 6. 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.

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 the x 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.

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 Table 6.

In 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 Table 7.

(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.

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 the mtr + 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.

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.

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 Figure 11 (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.

The network comprises 5 subnetworks, i.e., bicycle, car, car hailing, metro, and bus subnetwork. In Figure 12, 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).

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 number K 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.

(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.

In the optimal scheme (see Figure 15), 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.

(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.

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-mode c 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.

The flow of each mode in the optimal scheme is shown in Figure 17. 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.

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.

In Figure 18(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.

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.

Appendix

A. Proof of Proposition 1

The unified equilibrium conditions are derived based on a similar proof in García and Marín [83].

When C2 holds, referring to García and Marín [83], assuming is the log-sum of its lower-level nests, we haveand it can be derived as

Substitute equation (A.2) in equation (6), then

Change n to k − 1 and substitute in equation (A.3), then

Change n to k − 2 in equation (A.3) and substitute it in equation (A.4) to remove , then repeat the process, and then we have

Like equation (A.3), when C1 holds, assuming is the log-sum of , we have

Substitute equation (A.6) in equation (A.5), then we have

When C3 holds, implying that if , then ; otherwise , . Substitute it in equation (A.7), then we have

For path , , , and are the same, denote as and as . Since , , and are flow related, denote the equilibrium flow vector of all OD pairs is , then we have the unified equilibrium conditions.

Proposition A.1. A vector is an equilibrium flow vector only if there exists the set of value , and , so that

Equation (A.9) can be reformulated as a nonlinear complement problem, which is

According to Nagurney [84]; the nonlinear complement problem and the VI problem (7) have precisely the same solutions, Proposition 1 is proved.

B. Extracts of the Algorithm

The extracts of the algorithm are reported in the following flowchart, see Figure 20.

The major steps include those as follows:Step 1: Initialization. Call the initialization operator to generate the first Pop with n solutions, where one solution consists of the bus route and frequency, the layout of parking lot, and bicycle stations.Step 2: Crossover and Mutation. Call the crossover and mutation operators for Pop to generate Pop1. Obtain Pop2 by combining the Pop and Pop1 together.Step 3: Infeasible Management. Classify each solution in Pop2 into feasible and infeasible subset; for infeasible subset, call the repair operator to fix the solutions.Step 4: GMS/TAP Process. Each solution in the feasible and infeasible subset should go through a GMS/TAP process. Obtain the NL mode-choice structure of each OD pair, then call the Algorithm 3 to output the objective functions of each solution.Step 5: Biased Fitness. Calculate the Hamming distance of feasible and infeasible subset, and obtain the biased fitness by adding the hamming distance and objective function. Select the survivors and set the iteration number t = t + 1;Step 6: Diversification Mechanism. Determine whether the algorithm should end and whether a diversification should be conducted.

(1)Initialization: Generate the path set using hybrid K-shortest path method. Set t = 1 and generate a random initial flow vector F (t);
(2)For all OD pairs, build the mode-choice structure, and .
(3)while
(4) Set as a zero vector. Calculate path utility vector C (t) based on F (t).
(5)  For all OD pairs, for k from 0 to K, find the minimum path utility of paths with k transfers from C (t) and denote as . Calculate according to equation (1).
(6)  For all nests in , find the minimum path utility that belongs to nest from C (t) and denote as . Calculate all the according to equation (2).
(7)  For each mode nest , assign all the flow to the path with minimum path cost, and denote the flow vector as .
(8) Update
(9)t = t + 1.
end while

Data Availability

The data can be achieved by contacting the first author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (no. 51638004, U20A20330, and 52131203).