Optimization of Multimodal Discrete Network Design Problems Based on Super Networks

: In this paper, we investigate the multimodal discrete network design problem that simultaneously optimizes the car, bus, and rail transit network, in which inter-modal transfers are achieved by slow trafﬁc modes including walking and bike-sharing. Speciﬁcally, a super network topology is presented to signify the modal interactions. Then, the generalized cost formulas of each type of links in the super network are deﬁned. And based on the above formulas a bi-objective programming model is proposed to minimize the network operation cost and construction cost with trafﬁc ﬂow equilibrium constraints, investment constraints and expansion constraints. Moreover, a hybrid heuristic algorithm that combines the minimum cost ﬂow algorithm and simulated annealing algorithm is presented to solve the proposed model. Finally, the effectiveness of the proposed model and algorithm is evaluated through two numerical tests: a simple test network and an actual multimodal transport network.


The Multimodal Network Design Problem
There are a few studies about the multimodal DNDP, which mainly focus on two aspects: one is how to optimize the bus network considering congestion interaction with other modes. The other one is how to optimize both the car network and the bus network.
For the first aspect, Wan and Lo studied transit network design problem in the multimodal network and presented a systematic phase-wise methodology that considering both the effect of congestion and integration of modal transfers [30]. Moreover, they modeled inter-route and inter-modal transfers through the state augmented multimodal (SAM) network approach. Uchida et al. established a multimodal network design problem model to solve the problem of optimal frequency design for bus services. They applied sensitivity analysis to define linear approximation functions between the Probit SUE link flows and the design parameters, which are then used as constraints in the sub-problem of the NDP instead of the original stochastic user equilibrium (SUE) condition [31]. Li et al. considered the impact of bike-sharing on transfer, optimized the design problem of the bimodal transit network, and proposed a continuum model to optimize the bimodal transit system and shared bikes in a grid network at the same time [32].
For the second aspect, Gallo et al. optimized both car and bus networks by increasing the capacity of existing roads and the service frequency of bus routes, but the network topology remained unchanged and the added optimization was not considered [33]. Zhang et al. simultaneously optimized the auto network and bus network by establishing a single-level mathematical program with complementarity constraints and solved the model with an active-set algorithm, however, they only considered two transportation modes without considering the modal interactions [34].
Cai et al. reviewed the research status of cycling network planning and design in cold-climate cities in China and worldwide and summarized the relevant methods that can be used in China [35]. Liu et al. investigated the optimal network design problem of bike paths, which are on or adjacent to roadways but are physically separated from motorized traffic within the existing urban network. They proposed a mixed-integer nonlinear nonconvex model for the problem and solved the proposed model with a global optimization method and a matheuristic [36]. In addition, in this paper, based on the optimal scheme of the bike network design problem, we use the known bike routes and walking routes to connect car, bus, rail transit networks as a multimodal transportation network, then study the multimodal transportation DNDP.
In summary, we consider the interactions among various transportation modes and comprehensively optimize the multimodal DNDP of the car, bus, and rail transit, which is meaningful and innovative. The main contributions are listed as follows: Firstly, the slow travel modes in the multimodal transportation network are considered in this paper, which is usually neglected in existing studies. Apart from walking, bike-sharing has become another popular transfer mode. Therefore, we take the bicycle and walking as transfer modes to describe the modal interactions of car, bus, and rail transit. Secondly, the car, bus, and rail transit network are simultaneously studied for the optimization of multimodal DNDP. Most studies only considered the car network and bus network. However, owing to the large capacity, fast speed, and high comfort of rail transit, it occupies a very important position in multimodal travel behavior. Therefore, it is necessary to consider rail transit into multimodal DNDP. Moreover, a multimodal super network topology is presented, which clearly expresses inter-modal transfers of car, bus, and rail transit networks, as well as inter-route transfers of the bus and rail transit network. Besides, a bi-objective programming model of multimodal DNDP is proposed to make decisions about expanding or adding links in the car, bus, or rail transit network. It aims to minimize the network operation cost and construction cost. Lastly, a hybrid heuristic algorithm based on the minimum cost flow algorithm and the simulated annealing algorithm is developed to efficiently and effectively solve the proposed model. The proposed algorithm in this paper and the branch and bound algorithm proposed by reference [14] are all suitable for solving combinatorial optimization problems. However, the branch and bound algorithm needs more branch operations and consumes more time. Compared with the algorithms proposed by references [18,19], they are improved based on the traditional DNDP model and solve the global optimal solution through the Mix Integer Linear Programming (MILP) solver. While we use mathematical planning methods and network flow theory to establish the multimodal DNDP model and apply heuristic algorithms to efficiently solve the optimal solution. They did not make use of the network structure of relevant problems. In addition, the MILP solver cannot solve the problem of a large number of integer variables to the optimal in a reasonable time, which thus hinders the computational efficiency.
The rest of the paper is organized as follows: In Section 2 a multimodal super network topology is presented. In Section 3 a bi-objective programming optimization model of multimodal DNDP is proposed. In Section 4 a hybrid heuristic algorithm to solve the proposed model is developed. In Section 5 two numerical tests are implemented to validate the proposed model and algorithm. Concluding remarks and possible future research directions are presented in the last section.

Multimodal Network Representation
There are 5 kinds of transportation modes considered in this paper: car, bus, rail transit, bicycle, and walking. Multimodal transportation weighted directed graph G = (V, A, W, K). is presented, in which some notations and definitions are shown in Table 1. Table 1. The notations and definitions of the multimodal transportation network.

Notations Definitions
K Set of transportation modes. It includes car, rail transit, and bus, which are indexed by c, r, and p respectively. Thus for any mode k ∈ K = {c, r, p}.
V Set of nodes, node i ∈ V. In a multimodal transportation network, nodes have multiple attributes, which can represent bus stations, rail transit stations, road intersections, and interchange hubs.
A Set of links, link (i, j) ∈ A. In this paper, the links in topological networks are divided into four types: A E ( the set of entering links), A L (the set of leaving links), A D (the set of driving links), A T (the set of transfer links). It is A E,k Set of the entering links of transportation mode k. It means that passengers begin to enter the network.
A L,k Set of leaving links of transportation mode k. It indicates that passengers' leave behavior from the network.
A D,k Set of driving links of transportation mode k. It expresses that transportation mode k is driving in the network.
A T,k Set of links for transferring to transportation mode k. It includes the transfers between the same or different transportation modes.
W Set of weights, weight w ij ∈ W. The weight of links in the multimodal network has a variety of attributes, which can be represented as a generalized cost formula.
The relationship among the elements in the multimodal transportation super network studied in this paper is shown in Figure 1. , is driving in the network. , Set of links for transferring to transportation mode . It includes the transfers between the same or different transportation modes. Set of weights, weight ∈ . The weight of links in the multimodal network has a variety of attributes, which can be represented as a generalized cost formula.
The relationship among the elements in the multimodal transportation super network studied in this paper is shown in Figure 1.

Procedure of Establishing a Multimodal Super Network Topology
In this section, we present a multimodal transportation network that meets two requirements: one expresses the transfer among car, bus and rail transit modes; the other one shows the internal transfer of bus or rail transit mode. As shown in Figure 2, a traveler departs from the origin to the destination . The transfer relationship among the transportation modes are described as follows: Firstly, passengers on bus line 2 can transfer to rail transit line 1 by transfer node 3, which represents the inter-modal transfer between bus and rail transit. Secondly, node 6 is the transfer station between bus line 2 and bus line 3, which stands for the inter-route transfer in the bus network. Moreover, passengers can drive from node 7 to node 10 and transfer to rail transit line 2 through the transfer link (9,10), which indicates the inter-modal transfer between car and rail transit. Besides, by the transfer node 12 passengers can transfer from rail transit line 1 to rail transit line 2, which shows the inter-route transfer in the rail transit network. Finally, all transfer nodes 3→3′, 6→6′, and 12→12′ are achieved by walking; and transfer links (9,10) and (16,17) are achieved by bike-sharing. Next, based on the super network theory, the process of establishing a multimodal transportation super network topology is as described in the following section.

Procedure of Establishing a Multimodal Super Network Topology
In this section, we present a multimodal transportation network that meets two requirements: one expresses the transfer among car, bus and rail transit modes; the other one shows the internal transfer of bus or rail transit mode. As shown in Figure 2, a traveler departs from the origin O to the destination D. The transfer relationship among the transportation modes are described as follows: Firstly, passengers on bus line 2 can transfer to rail transit line 1 by transfer node 3, which represents the inter-modal transfer between bus and rail transit. Secondly, node 6 is the transfer station between bus line 2 and bus line 3, which stands for the inter-route transfer in the bus network. Moreover, passengers can drive from node 7 to node 10 and transfer to rail transit line 2 through the transfer link (9,10), which indicates the inter-modal transfer between car and rail transit. Besides, by the transfer node 12 passengers can transfer from rail transit line 1 to rail transit line 2, which shows the inter-route transfer in the rail transit network. Finally, all transfer nodes 3 → 3 , 6 → 6 , and 12 → 12 are achieved by walking; and transfer links (9,10) and (16,17) are achieved by bike-sharing. Next, based on the super network theory, the process of establishing a multimodal transportation super network topology is as described in the following section.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 6 of 27 Figure 2. The diagram of the multimodal transportation network.

Bus Network
The bus network is illustrated in Figure 3, which includes three bus routes, namely bus line 1, bus line 2, and bus line 3, where the pink node 6 is the transfer station. The information of routes in the bus network is shown in Table 2.  The bus network is illustrated in Figure 3, which includes three bus routes, namely bus line 1, bus line 2, and bus line 3, where the pink node 6 is the transfer station. The information of routes in the bus network is shown in Table 2.

Bus Network
The bus network is illustrated in Figure 3, which includes three bus routes, namely bus line 1, bus line 2, and bus line 3, where the pink node 6 is the transfer station. The information of routes in the bus network is shown in Table 2.  The transfer node 6 in Figure 3 can be extended to nodes 6 and 6′ through super network expansion, so Figure 3 becomes the bus super network topology, which is shown in Figure 4.  The transfer node 6 in Figure 3 can be extended to nodes 6 and 6 through super network expansion, so Figure 3 becomes the bus super network topology, which is shown in Figure 4.

Rail Transit Network
The rail transit network is shown in Figure 5, which contains two rail transit routes, namely rail transit line 1 and rail transit line 2, where the pink node 12 is the transfer station. The information of routes in the rail transit network is shown in Table 3.

Rail Transit Network
The rail transit network is shown in Figure 5, which contains two rail transit routes, namely rail transit line 1 and rail transit line 2, where the pink node 12 is the transfer station. The information of routes in the rail transit network is shown in Table 3.

Rail Transit Network
The rail transit network is shown in Figure 5, which contains two rail transit routes, namely rail transit line 1 and rail transit line 2, where the pink node 12 is the transfer station. The information of routes in the rail transit network is shown in Table 3.  Like transfer node 6 discussed above, the transfer node 12 can be extended to nodes 12 and 12′ by super network expansion. As a result, the rail transit super network topology is as shown in Figure 6.

Car Network
The car network is shown in Figure 7. It contains three car routes, named car route 1, car route 2, and car route 3, respectively. The information of each route in the car network is shown in Table 4. Like transfer node 6 discussed above, the transfer node 12 can be extended to nodes 12 and 12 by super network expansion. As a result, the rail transit super network topology is as shown in Figure 6.
station. The information of routes in the rail transit network is shown in Table 3.  Like transfer node 6 discussed above, the transfer node 12 can be extended to nodes 12 and 12′ by super network expansion. As a result, the rail transit super network topology is as shown in Figure 6.

Car Network
The car network is shown in Figure 7. It contains three car routes, named car route 1, car route 2, and car route 3, respectively. The information of each route in the car network is shown in Table 4.

Car Network
The car network is shown in Figure 7. It contains three car routes, named car route 1, car route 2, and car route 3, respectively. The information of each route in the car network is shown in Table 4.   In general, for car mode passengers or drivers needn't transfer in one trip, so it is unnecessary to expand nodes in the car network, the car super network topology is shown in Figure 8.  In general, for car mode passengers or drivers needn't transfer in one trip, so it is unnecessary to expand nodes in the car network, the car super network topology is shown in Figure 8. In general, for car mode passengers or drivers needn't transfer in one trip, so it is unnecessary to expand nodes in the car network, the car super network topology is shown in Figure 8. In summary, by connecting the same nodes among different modes, a multimodal transportation super network topology ( Figure 9) is obtained, of which three layers are bus, rail transit, and car network topology respectively. In Figure 9, the pink dotted links clearly illustrate the inter-modal transfers among bus, rail transit, and car, as well as the inter-route transfers of bus or rail transit. In summary, by connecting the same nodes among different modes, a multimodal transportation super network topology ( Figure 9) is obtained, of which three layers are bus, rail transit, and car network topology respectively. In Figure 9, the pink dotted links clearly illustrate the inter-modal transfers among bus, rail transit, and car, as well as the inter-route transfers of bus or rail transit.

Formulating the Multimodal Discrete Network Design Problem
In this section, firstly the generalized cost formulas of each type of links are defined. Then, a bi-objective programming optimization model of multimodal DNDP will be introduced.

The Generalized Cost Formula
Many factors affect the impedance of links in the multimodal networks. In this paper, the impact of travel time, monetary cost, and comfort are mainly considered. Among them, travel time includes walking time, driving time, and waiting time. Monetary cost can be divided into two types: mileage-based fare (such as the car mode) and route-based

Formulating the Multimodal Discrete Network Design Problem
In this section, firstly the generalized cost formulas of each type of links are defined. Then, a bi-objective programming optimization model of multimodal DNDP will be introduced.

The Generalized Cost Formula
Many factors affect the impedance of links in the multimodal networks. In this paper, the impact of travel time, monetary cost, and comfort are mainly considered. Among them, travel time includes walking time, driving time, and waiting time. Monetary cost can be divided into two types: mileage-based fare (such as the car mode) and route-based fare (such as the bus and rail transit mode).
The perception of comfort varies with the travel time due to different transportation modes. Besides, we reserve some time to prevent emergencies. Consequently, the generalized cost of links is mainly composed of travel time, monetary cost, and comfort loss, risk reserve time, which is specifically quantized as: where α, β, γ, δ: the weight coefficients and we have α + β + γ + δ = 1, T: the travel time cost, P: the corresponding time cost converted from currency conversion, S: the corresponding time cost converted from the loss of comfort and R: the reserved time for risk.

The Generalized Cost Formulas of Entering Links and Leaving Links
(1) Travel time In the actual transportation system, entering links A E,k , ∀k ∈ {r, p} indicate the process of reaching the stations, which mainly include walking time and waiting time; entering links A E,k , k = {c} express the process of reaching the parking lot, which mostly includes walking time. Therefore, it can be concluded that the travel time of the entering links mainly consists of walking time and waiting time. Consequently, the travel time of the entering link can be formulated as: where t E,k walk,ij : the average walking time of the entering links and t E,k wait,ij : the waiting time of the entering links. In this paper, it is set that t E,k wait,ij = 1 2 f , ∀(i, j) ∈ A E,k , k ∈ { p, r} , f is the departure frequency, and t E,k wait,ij = 0, ∀(i, j) ∈ A E,k , k = { c} . Meanwhile, in the actual transportation system, leaving links A L,k (∀k ∈ K) show the process of arriving at the destination, which mostly contains walking time. Hence, the travel time of the leaving links can be described as: where t L,k walk,ij is the average walking time of leaving links.

(2) Monetary cost
Since there is no monetary cost on the entering links and leaving links, we set the monetary cost to zero, which is formulated as: where P E,k ij : the monetary cost of entering links and P L,k ij : the monetary cost of leaving links.

(3) Comfort loss
Since the trip is very short and the value of comfort loss is very small, we also set the comfort loss to zero on the entering links and leaving links. The equation is as follows: where S E,k ij : the comfort loss of entering links and S L,k ij : the comfort loss of leaving links. (4) Risk reserve time According to travel time, so we set the risk reservation time of the entering links as: where ρ E,k ij : the delay parameter of the entering links, and ρ E,k ij > 1, ∀(i, j) ∈ A E,k , k ∈ K. Similarly, the risk reservation time of the leaving links is formulated as: where ρ L,k ij : the delay parameter of leaving links and ρ L,

The Generalized Cost Formulas of Driving Links
(1) Travel time In the actual transportation system, driving links A D,k , ∀k ∈ {c}, stand for the road between adjacent intersections. In this paper, we set the travel time of car driving links as: where t D,c0 ij : the free-flow travel time, x ij : the newly assigned traffic flow, M D,c ij : the original traffic flow, c D,c ij : the capacity of driving links, ξ: the average number of passengers and ϕ, σ: the undetermined coefficient.
Meanwhile, in the actual transportation system, driving links A D,k , ∀k ∈ {r, p} indicate adjacent stations, and the travel time of the bus and rail transit driving links are set as: where t D,k ij : the average travel time of the transportation mode k.

(2) Monetary cost
Since the car is charged depending on mileages, in this paper, the monetary cost of the car driving links is described as: where p D,c0 a : the monetary cost per unit mileage, L D,c ij : the length of the link (i, j). Moreover, as the price of the ticket is up to mileages, the monetary cost of the bus and rail transit driving links are set as: where p D,k0 ij : the starting fare, L D,k0 : the number of miles corresponding to the starting fare of the mode k, L D,k ij : the length of the link (i, j) and p D,k ij : the increasing ticket price per unit mileage.

(3) Comfort loss
As the comfort loss is proportional to the travel time on the car network, it can be formulated as: where: γ D,c : the degree of comfort loss per unit time of the car mode.
In addition, the comfort loss is mainly determined by travel time and congestion on the bus or rail transit network, so it can be described as: where η D,k 0 : the comfort loss parameter when the transportation mode k is empty, η D,k 1 : the comfort loss parameter when congestion occurs, q D,k ij : the number of passengers, f D,k ij : the departure frequency and o D,k ij : the capacity of vehicles. (4) Risk reserve time Similarly, we set the risk reservation time of the driving links as follows: where ρ D,k ij : the delay parameter of driving links and ρ D,

The Generalized Cost Formulas of Transfer Links
(1) Travel time It has been mentioned above that the transfer modes mainly consist of walking and bike-sharing. Therefore, the travel time of transferring to the transportation mode k can be formulated as: or: where t T,k walk,ij : the average walking time for transfer to the transportation mode k, t T,k cycling,ij : the average cycling time for transfer to the transportation mode k, t T,k wait,ij : the waiting time for transfer to the transportation mode k, and we set t T,k Since there is no monetary cost for walking, and only the monetary cost of bike-sharing exists, the monetary cost of transferring to the mode k can be represented as: or: where p T,k0 ij : the starting fare for cycling, L T,k0 : the number of miles corresponding to the starting fare for cycling, L T,k ij : the length of the link (i, j) and p T,k ij : the increasing fare per unit mileage.

(3) Comfort loss
Since the transfer trip is very short and the comfort loss can be ignored, we set the comfort loss of transferring to the transportation mode k as: Similarly, the risk reservation time of transferring to the transportation mode k is expressed as: where ρ T,k ij : the delay parameter of transferring to the transportation mode k and ρ T,k

Bi-Objective Programming Model
The objective of multimodal DNDP is to decide on an optimal scheme U * from the candidate schemes U = {U c , U r , U p } to minimize the total cost, which includes the network operation cost and the construction cost of the optimization scheme. U is the set of candidate links; U c is the set of added or expanded links for the car network; U r is the set of added or expanded links for the rail transit network, and U p is the set of added or expanded links for the bus network.
(1) Network operation cost According to the generalized cost formulas mentioned above, the network operation cost of the multimodal transportation network can be formulated as follows: Here, x k ij is the traffic flow of the link (i, j).

(2) Construction cost
For DNDP, there are two methods to optimize the network: one is to add a new link; the other one is to expand the original link. In the multimodal super network topology, an added link is to connect two unconnected nodes, while an expanded link is to add a link to the originally connected nodes.
As shown in Figure 10: for the car network, since node 3 and node 14 were not connected originally, the added link (3,14) (the blue dotted link) represents the added optimization. Meanwhile, as node 7 and node 10 were connected originally, so link (7,10) (the blue dotted link) illustrates the expanded optimization. Therefore, we use the 0-1 variable y c ij and z c ij to respectively express whether to add or expand links. If the link (i, j) is added/expanded, y c ij /z c ij is equal to 1, otherwise y c ij /z c ij is equal to 0, ∀(i, j) U c For the bus network, the added method indicates adding a new link on the existed route, like link (16,17) (the green dotted link). The expanded method means to increase the frequency of departure, so it will change the whole route, such as link (1,2), (2,3), (3,4), (4,5), (5,6) (the green dotted link). Hence, if the link (i, j) is added/expanded the 0-1 variable y p ij /z p ij is equal to 1, otherwise y p ij /z p ij is equal to 0, ∀(i, j) U p in the bus network. Similarly, for the rail transit network, the added method expresses to add a new link on the existed route, like link (1,3) (the red dotted link). And the expanded method represents to increase the frequency of departure, so it will change the whole route, such as links (10,12), (12,13), (13,14) (the red dotted link). Thus, if the link (i, j) is added/expanded the 0-1 variable y r ij /z r ij is equal to 1, otherwise y r ij /z r ij is equal to 0, ∀(i, j) U r in the rail transit network.
is the traffic flow of the link ( , ).
(2) Construction cost For DNDP, there are two methods to optimize the network: one is to add a new link; the other one is to expand the original link. In the multimodal super network topology, an added link is to connect two unconnected nodes, while an expanded link is to add a link to the originally connected nodes.
As shown in Figure 10: for the car network, since node 3 and node 14 were not connected originally, the added link (3,14) (the blue dotted link) represents the added optimization. Meanwhile, as node 7 and node 10 were connected originally, so link (7,10) (the blue dotted link) illustrates the expanded optimization. Therefore, we use the 0-1 variable and to respectively express whether to add or expand links. If the link ( , ) is added/expanded, / is equal to 1, otherwise / is equal to 0, ∀( , )ϵ . For the bus network, the added method indicates adding a new link on the existed route, like link (16,17) (the green dotted link). The expanded method means to increase the frequency of departure, so it will change the whole route, such as link (1,2), (2,3), (3,4), (4,5), (5,6) (the green dotted link). Hence, if the link ( , ) is added/expanded the 0-1 variable / is equal to 1, otherwise / is equal to 0, ∀( , ) in the bus network. Similarly, for the rail transit network, the added method expresses to add a new link on the existed route, like link (1,3) (the red dotted link). And the expanded method represents to increase the frequency of departure, so it will change the whole route, such as links (10,12), (12,13), (13,14) (the red dotted link). Thus, if the link ( , ) is added/expanded the 0-1 variable / is equal to 1, otherwise / is equal to 0, ∀( , ) in the rail transit network. Consequently, the construction cost of the optimization scheme is formulated as: where B k1 ij : the cost of the added link (i, j) for the mode k, ∀k ∈ K and B k2 ij : the cost of the expanded link (i, j) for the mode k, ∀k ∈ K.
(3) Bi-objective programming model In summary, a bi-objective programming model for multimodal DNDP is proposed based on Equations (21) and (22) to minimize the network operation cost and construction cost as follows: where θ, τ: the weight coefficients and we have θ + τ = 1, BG k : the budget cost of the transportation mode k, M k ij ,c k ij : the original traffic flow and the original capacity of the link (i, j) respectively, q rs : the traffic demand of OD pair rs, c k1 ij /c k2 ij : the capacity of the added/expanded link (i, j) for the mode k, ω: the binary variable, if the link (i, j) is expanded, ω is equal to 1, otherwise, ω is equal to 0 on the bus or rail transit route and N: the number of expanded links on the bus or rail transit route.
Constraints (24)- (26) indicate that the cost of optimization schemes is less than the budget cost for car, bus, and rail transit network.
Constraint (27) shows that the sum of the assigned traffic flow and the original traffic flow is less than the sum of the original capacity and the increased capacity of the added or expanded optimization. Constraints (28)-(31) are traffic flow equilibrium constraints, where constraint (28) means that for the starting point, the outflow minus the inflow is equal to q rs . Constraint (29) shows that for the destination point, the outflow minus the inflow is equal to −q rs . Constraint (30) expresses that for each node, the outflow and inflow must be equal. Constraint (31) indicates that the traffic flow cannot be less than zero.
Constraint (32) ensures that the expanded method will change the whole existed route for the bus and rail transit network. Specifically, if one link of the route is expanded, then the other links must be expanded.

Solution Algorithm
Since there are two types of decision variables in the proposed optimization model: the assignment of traffic flow and the decision of adding or expanding candidate links, the algorithm is divided into two parts: Firstly, we need to assign the traffic flow to minimize the network operation cost. Since the capacity of each link is considered in the proposed model, it is appropriate that the minimum cost flow algorithm is applied to solve this problem.
Secondly, as multimodal DNDP is an integrated optimization problem, conventional solving algorithms are difficult to solve efficiently, heuristic algorithms are developed in this paper: simulated annealing algorithm to search the near-optimal decision U * from the candidate schemes U = {U c , U r , U p }.

Minimum Cost Flow Algorithm
Let µ k be the minimum cost augmentation chain of the k-th iteration, ε k be the adjustable flow on the augmentation chain µ k , and x k be the feasible flow of the k-th iteration. Based on the shortest path method, a network L(x) is constructed, whose nodes are the nodes of the original network G, specifically, the link (v i , v j ) of the original network G is changed into two links with opposite directions (v i , v j ) and v j , v i . Define the weights of links (v i , v j ) and v j , v i in L(x) as l ij and l ji : Set the matrix D = d ij n×n , d ij represents the shortest length from v i to v j in the original network G; S = s ij n×n , s ij represents the end number of the first link of the Above all, the algorithm for solving the minimum cost feasible flow with the target flow q rs is as follows: Step 1: Set k = 0 Start with a minimum cost initial feasible flow x (0) , x (0) ≤ q rs (generally, x (0) = 0).
Step 2: Generate L x (k) , find the minimum cost path from v r to v s in L x (k) , and use the Floyd algorithm to solve the shortest path.
Step 2.1: Set m = 0 Step 2.2: m = m + 1 Calculate: where: d If m = n go to step 2.2, else go to step 2.3.
Step 2.3: when m = n, the algorithm ends.
ij is the shortest path length from v i to v j ; ij is the end node of the first link in the shortest path from v i to v j .
Step 3: If L x (k) does not have the shortest path from v r to v s , then x (k) is regarded as the maximum flow, there is no minimum cost flow with the target flow of q rs , end the algorithm; otherwise, go to the next step.
Step 4: Find the augmented chain µ k corresponding to the shortest path from v r to v s in L x (k) in the original network G, and calculate ε k = min min c ij − x k ij , min x k ij . If ε k ≤ q rs − x (k) , according to the maximum flow adjustment scheme, adjust the augmented flow ε k on the augmented chain µ k to obtain a new flow x = x (k) + ε k , and let k = k + 1, assign new flow x to x k , go to step 2; otherwise, when ε k > q rs − x (k) , adjust the augmented flow q rs − x (k) on the augmented chain µ k to obtain a new flow x , where x is the minimum feasible flow of the target flow q rs . End the algorithm.

Simulated Annealing Algorithm
(1) Determination of solution space In this problem, all the possible added and expanded candidate schemes U = {U c , U r , U p } constitute the solution space, and a specific candidate scheme can be regarded as a partition of the set U = {U 1 , . . . , U m }, where ∪ m i=1 U i = U, U i ∩ U j = ∅, i = j, and U i or U j is the feasible solution. The solution space is composed of all possible partitions, that According to the discussion presented in Section 3.5, the objective function of the optimization schemes {U 1 , · · · , U m } can be defined as: (

3) The choice of initial solution
The initial solution is the initial point of the simulated annealing algorithm. A large number of experiments show that the final solution of the simulated annealing algorithm does not depend on the selection of the initial solution, so we randomly choose an initial as U i , i ∈ {1, · · · , m}.

(4) The generation of a new solution
For this algorithm, we first need to use a certain strategy to generate a new solution based on the current solution in the space, then decide whether to accept the new solution. When the new solution is accepted, which will be treated as a new current solution, and the value of the objective function is modified correspondingly. The next iteration is started on this basis. If the new solution is not accepted, the next experiment will be continued based on the current solution.
Since the proposed optimization model is a binary integer programming, we can arbitrarily choose u ∈ U i . If u is 0, it will be replaced by 1; otherwise, if it is 1, it will be replaced by 0. Therefore, a new solution U j is generated and the difference of values of the objective function in two iterations can be calculated by the following formula:

(5) Acceptance criterion
An acceptance criterion is set to judge whether the new solution is accepted or not. Generally, the most commonly used acceptance criteria is Metropolis criteria: where t is the control parameter. According to the settings above, the procedure of the simulated annealing algorithm is listed as follows: Step 1: Randomly select an initial solution U i , let the current solution U 0 = U i ; the current number of iteration steps k = 0; the current temperature t k = t max .
Step 2: When reaching the maximum number M of iterations the inner loop, go to step 3; otherwise, randomly select a neighbor U j from the neighborhood N(U i ) and calculate f ji , if f ji ≤ 0 then U i = U j , otherwise, if exp − f ji /t k > random(0, 1) (represents a uniform random number between 0 and 1), then U i = U j , repeat step 2.
Step 3: k = k + 1, t k+1 = d(t k ) (representing the formula of temperature drop), and the classic simulated annealing algorithm cooling formula is d(t k ) = t 0 1+t k . If the end temperature t end , is reached, go to Step 4, otherwise, go to Step 2.
Step 4: Output the calculation result and stop.
As shown in Figure 11, the simulated annealing algorithm includes an inner loop and an outer loop. The inner loop is the second step, which means random search in the neighborhood N(U i ) at the same temperature t. The outer loop mainly includes the temperature drop change t k+1 = d(t k ), the number of iteration steps k = k + 1, and the stop condition in the third step.
Step 4: Output the calculation result and stop.
As shown in Figure 11, the simulated annealing algorithm includes an inner loop and an outer loop. The inner loop is the second step, which means random search in the neighborhood ( ) at the same temperature . The outer loop mainly includes the temperature drop change +1 = ( ), the number of iteration steps = + 1, and the stop condition in the third step.  Consequently, we combine the minimum cost flow algorithm and simulated annealing algorithm as a hybrid heuristic algorithm that can solve the proposed model. The relationship between the minimum cost flow algorithm and the simulated annealing algorithm is shown in Figure 11. The minimum cost flow algorithm will be used to assign traffic flow before solving the objective function in each iteration of the simulated annealing algorithm.

Numerical Examples
In this section, two numerical tests are implemented to validate the proposed model and algorithm: a simple test network (the multimodal transportation network shown in Figure 2) and an actual network (an area of Beijing, China). All experiments are executed in Python software on a Lenovo laptop configured with an Intel Core i5-8250U 3.40 GHz CPU equipped with 8 GB RAM.

Simple Network
Before optimization, the maximum feasible flow of the network is 78, and the traffic demand q rs is 100, so it is necessary to improve the network. According to the passenger demand, some candidate links are given, as shown in Figure 12: the blue dotted links (3,14), (7,10) are the candidate links for the car network, namely U c ={(3,14), (7,10)}; the green dotted link (16,17) is the candidate link for the bus network, namely U p ={ (16,17)}; the red dotted link (1,3) is the candidate link for the rail transit network, namely U r ={(1,3)}. Therefore, it can be concluded that the set of candidate schemes U = {U c , U r , U p }= { (3,14), (7,10), (16,17), (1,3)}, so there are in total 2 4 schemes. Next, we will use the exact algorithm and heuristic algorithm respectively to find the near-optimal scheme U * from the 2 4 schemes. (3,14), (7,10) are the candidate links for the car network, namely ={ (3,14), (7,10)}; t green dotted link (16,17) is the candidate link for the bus network, namely ={(16,17 the red dotted link (1,3) is the candidate link for the rail transit network, name ={(1,3)}. Therefore, it can be concluded that the set of candidate schemes { , , }= {(3,14), (7,10), (16,17), (1,3)}, so there are in total 2 4 schemes. Next, we will u the exact algorithm and heuristic algorithm respectively to find the near-optimal schem * from the 2 4 schemes. (1) Exact algorithm To test the feasibility of the developed heuristic algorithm, an exact algorithm is fi used to solve the model. The specific steps are as follows: Step 1: Calculate the generalized cost of each type of links (1) Exact algorithm To test the feasibility of the developed heuristic algorithm, an exact algorithm is first used to solve the model. The specific steps are as follows: Step 1: Calculate the generalized cost of each type of links Firstly, give the information of links for each network in the multimodal transportation, among which, the car network mainly includes length (km), zero traffic flow driving time (min), existing vehicle (pcu/h), capacity (pcu/h); the bus network and rail transit network mainly include average driving time (min), frequency (number/ hour); transfer network mainly includes length (km), walking time (min) and cycling time (min). Then, calculate the generalized cost of each type of links according to the formulas in Section 3, and the results are shown in Figure 13: in brackets, the first term is the generalized cost of links and the second term is the capacity of links.
Step 2: Traffic flow assignment The minimum cost flow algorithm is coded by PYTHON based on the steps in Section 4.1 to assign the traffic flow. The traffic flow assignment and the composition of transportation mode for each shortest path are obtained, as shown in Table 5, in which all candidate links are taken. Firstly, give the information of links for each network in the multimodal transportation, among which, the car network mainly includes length (km), zero traffic flow driving time (min), existing vehicle (pcu/h), capacity (pcu/h); the bus network and rail transit network mainly include average driving time (min), frequency (number/ hour); transfer network mainly includes length (km), walking time (min) and cycling time (min). Then, calculate the generalized cost of each type of links according to the formulas in Section 3, and the results are shown in Figure 13: in brackets, the first term is the generalized cost of links and the second term is the capacity of links. Figure 13. The generalized cost and capacity of each link in the super network topology.
Step 2: Traffic flow assignment The minimum cost flow algorithm is coded by PYTHON based on the steps in Section 4.1 to assign the traffic flow. The traffic flow assignment and the composition of transportation mode for each shortest path are obtained, as shown in Table 5, in which all candidate links are taken.  Step 3: Calculate the value of the objective function The optimal result of traffic flow assignment in step 2 is used to calculate the value of the network operation cost by Equation (21); then based on the optimization cost and the budget cost of each candidate link to calculate the value of construction cost by Equation (22); finally, the weighted sum of the operation cost and the construction cost by Equation (23) is the value of the objective function. When the generalized cost formula coefficients α = 0.25, β = 0.25, γ = 0.25, δ =0.25 in Equation (1), and the objective function weight coefficients θ = 0.5, τ = 0.5 in Equation (23), the calculation results are illustrated in Table 6 and Figure 14.
It can be concluded from Table 6 and Figure 14 that the optimal solution of the exact algorithm is U * = (1,3), and the corresponding operation cost is 4114, construction cost is 3000, and total cost is 3557. Some schemes are without values of construction cost B and objective function Z in Table 6. Because the maximum feasible flow of the schemes does not reach q rs = 100, so the optimization schemes are not adopted, and there is no need to calculate the objective function.  0  0  0  0  78  3429  --1  0  0  0  90  4041  --0  1  0  0  78  3428  --0  0  1  0  100  4114  3000  3557  0  0  0  1  78  3365  --1  1  0  0  90  4028  --1  0  1  0  100  3958  4000  3979  1  0  0  1  90  3977  --0  of the network operation cost by Equation (21); then based on the optimization cost and the budget cost of each candidate link to calculate the value of construction cost by Equation (22); finally, the weighted sum of the operation cost and the construction cost by Equation (23) is the value of the objective function. When the generalized cost formula coefficients = 0.25, = 0.25, = 0.25, = 0.25 in Equation (1), and the objective function weight coefficients = 0.5, = 0.5 in Equation (23), the calculation results are illustrated in Table 6 and Figure 14.  0  0  0  0  78  3429  --1  0  0  0  90  4041  --0  1  0  0  78  3428  --0  0  1  0  100  4114  3000  3557  0  0  0  1  78  3365  --1  1  0  0  90  4028  --1  0  1  0  100  3958  4000  It can be concluded from Table 6 and Figure 14 that the optimal solution of the exact algorithm is * = (1,3), and the corresponding operation cost is 4114, construction cost is 3000, and total cost is 3557. Some schemes are without values of construction cost and objective function in Table 6. Because the maximum feasible flow of the schemes does (2) Heuristic algorithm Similarly, the developed heuristic algorithm is also coded by PYTHON based on the steps in Section 4.2 to solve the model, and we set the initial temperature t max = 500, the end temperature t end = 100, the maximum number of iterations M = 50. Besides, according to Equation (43), the relative error between the heuristic algorithm and the exact algorithm is calculated. The results are shown in Table 7.
GAP= (near-optimal solution − exact solution)/exact solution × 100% (43) Clearly, the solution calculated by the heuristic algorithm is the same as the exact algorithm, which verifies the effectiveness and feasibility of the heuristic algorithm. To further test the heuristic algorithm, then we carry out sensitivity analysis of parameters in Equations (1) and (23).

(3) Sensitivity analysis
In this section, the sensitivity analysis is carried out on the weight coefficients of the generalized cost formula Equation (1) α, β, γ, δ and the weight coefficients of the objective formula Equation (23) θ, τ. The specific results are listed in Table 8.
It can be seen from Table 8 that the results obtained by the heuristic algorithm are basically consistent with the exact solution, and only two calculation results deviate from the exact value. Moreover, the relative errors of the two results obtained by the two algorithms are within 1.5%, which validates the effectiveness and reliability of the heuristic algorithm within the allowable error range. In addition, the exact algorithm takes several hours to obtain the optimal solution, while the heuristic algorithm takes only within 56 s. Therefore, the heuristic algorithm is superior to the exact solution algorithm in terms of calculation time.

Actual Network
The actual transportation network in this case is an area from Xizhimen to Dongdan in Beijing, China. It is one of the busiest traffic areas in Beijing. The network of this area includes multiple bus, rail transit lines and integrates various travel modes, such as car, rail transit, bus, bicycle, and walking. It has many important transfer hub stations. So it is a typical multimodal transportation network. Moreover, there are many facilities for sharing bicycles near each transfer station, which facilitates passengers' transfer behavior. That is consistent with the model proposed in this paper, which takes bike-sharing as a transfer mode. In addition, due to the huge passenger demand from Xizhimen to Dongdan and the bottleneck of the transportation network, the network capacity is insufficient. Thus, it is necessary to expand or add some links for the car, bus, and rail networks respectively. This is also an important reason why this regional network is selected as an actual case. The multimodal transportation network diagram is as follows, which combined car, bus, and rail transit modes.
As shown in Figure 15: the origin is Xizhimen while the destination is Dongdan. The transfer relationship among modes is as follows, among them, the pink nodes are the transfer stations:  Besides, the inter-modal transfer between car and rail transit includes 1 route: passengers drive from the origin Xizhimen to node 13 and transfer to rail line 3 by node 13.
includes multiple bus, rail transit lines and integrates various travel modes, such as car, rail transit, bus, bicycle, and walking. It has many important transfer hub stations. So it is a typical multimodal transportation network. Moreover, there are many facilities for sharing bicycles near each transfer station, which facilitates passengers' transfer behavior. That is consistent with the model proposed in this paper, which takes bike-sharing as a transfer mode. In addition, due to the huge passenger demand from Xizhimen to Dongdan and the bottleneck of the transportation network, the network capacity is insufficient. Thus, it is necessary to expand or add some links for the car, bus, and rail networks respectively. This is also an important reason why this regional network is selected as an actual case. The multimodal transportation network diagram is as follows, which combined car, bus, and rail transit modes.
As shown in Figure 15: the origin is Xizhimen while the destination is Dongdan. The transfer relationship among modes is as follows, among them, the pink nodes are the transfer stations: Besides, the inter-modal transfer between car and rail transit includes 1 route: passengers drive from the origin Xizhimen to node 13 and transfer to rail line 3 by node 13. Before optimization, the maximum feasible flow of the network is 1524, and the traffic demand is 2000, so it is necessary to improve the network. According to the passenger demand, the candidate links are given, as shown by the dotted line in Figure 15: for the car network, the blue dotted links (1,2), (2,3) represent the expansion optimization, and the blue dotted link (3,7) is the candidate link for the adding optimization, namely ={(1,2), (2,3), (3,7)}; for the bus network, the green dotted links (12,13), (15,16), (17,20) are the candidate links for adding optimization, namely ={ (12,13), (15,16), (17,20)  The multimodal transportation super network topology is built and the link information of each network is given to calculate the generalized cost of each type of links, as shown in Figure 16: in brackets, the first term is the generalized cost of links and the second term is capacity of links.
Finally, the proposed heuristic algorithm will be coded by PYTHON for the numerical test based on Section 4, which combines the minimum cost flow algorithm and simulated annealing algorithm. And we set the initial temperature t max = 1000, the end temperature t end = 100, the maximum number of iterations M = 200. To analyze the optimal schemes under different weight coefficients, 12 groups of numerical experiments are carried out. The near-optimal schemes under different weight coefficients are shown in the table below. totally 2 7 schemes. Next, we use the heuristic algorithm to find the near-optimal scheme from the 2 7 schemes.
The multimodal transportation super network topology is built and the link information of each network is given to calculate the generalized cost of each type of links, as shown in Figure 16: in brackets, the first term is the generalized cost of links and the second term is capacity of links. Finally, the proposed heuristic algorithm will be coded by PYTHON for the numerical test based on Section 4, which combines the minimum cost flow algorithm and simulated annealing algorithm. And we set the initial temperature = 1000, the end temperature = 100, the maximum number of iterations = 200. To analyze the optimal schemes under different weight coefficients, 12 groups of numerical experiments are carried out. The near-optimal schemes under different weight coefficients are shown in the table below.
It can be concluded from Table 9 that the choice of the near-optimal scheme is affected by the weight coefficients: (1) When the value of is larger than in Equation (1), the travel time is the main priority, and passengers are more willing to choose the car mode that can reach the destination quickly. In this case, when = 0.6 (group 8,9,10), the car's candidate links (1,2), (2,3) are all adopted in various groups, which are consistent with the choice of passengers. On the contrary, the monetary cost is focused on, and passengers prefer to choose the transportation mode that can save a lot of money. While = 0.6 (group 11,12,13), the bus's candidate links (12,13), (15,16), (17,20) are all adopted in various groups, which are also tally with the actual situation. It can be concluded from Table 9 that the choice of the near-optimal scheme is affected by the weight coefficients: (1) When the value of ff is larger than fi in Equation (1), the travel time is the main priority, and passengers are more willing to choose the car mode that can reach the destination quickly. In this case, when ff = 0.6 (group 8,9,10), the car's candidate links (1,2), (2,3) are all adopted in various groups, which are consistent with the choice of passengers. On the contrary, the monetary cost is focused on, and passengers prefer to choose the transportation mode that can save a lot of money. While fi = 0.6 (group 11,12,13), the bus's candidate links (12,13), (15,16), (17,20) are all adopted in various groups, which are also tally with the actual situation (2) When the value of`is larger than ø in Equation (23), the decision-maker pays more attention to the network operation cost. In this case, when`= 0.9 (group 2,5,8,11), the car's candidate links (1,2), (2,3) are the principal part of the near-optimal scheme under different groups. Conversely, the construction cost is mainly considered. While ø = 0.9 (group 3,6,9,12), the bus's candidate links (12,13), (15,16), (17,20) are the main component of the near-optimal scheme in various groups. These are consistent with the actual situation, which proves that the developed heuristic algorithm can efficiently find the near-optimal solution with good stability and reliability.
Therefore, in an actual operation situation, the value of each parameter can be set according to the attention of the decision-maker, for example, when the travel time is more important, the value of α is larger than β, otherwise, the monetary cost is mainly considered; when decision-maker pays more attention to the operation cost, the value of θ is larger than τ, conversely, the construction cost is principal considered.

Conclusions
In this paper, we comprehensively optimize a multimodal transportation DNDP that is composed of car, bus, rail transit, bicycle, and walking. Firstly, a multimodal super network topology that includes four types of links (entering links, leaving links, driving links, transfer links) is proposed, which clearly expresses inter-modal transfers among car, bus, and rail transit networks, as well as inter-route transfers in the bus or rail transit network. And to complete the above transfers, walking and bike-sharing are considered as transfer modes. Specifically, in Section 2 we present a simple multimodal transportation network that meets two types of the above transfer relationships and we also discuss the process of establishing a multimodal transportation super network topology. Secondly, the generalized cost formulas of each type of links in the super network topology are defined. Among them, the generalized cost is mainly composed of four types, namely travel time, monetary cost, and comfort loss, risk reserve time. Based on the above formulas, a biobjective programming model to minimize the network operation cost and construction cost is proposed to decide whether to expand or add links for the car, bus, and rail network respectively. Moreover, based on the minimum cost flow algorithm and simulated annealing algorithm, we develop a hybrid heuristic algorithm to solve the proposed model, which can efficiently and effectively get the near-optimal solution. Finally, two numerical tests of a simple test network (designed in Section 2) and an actual network (an area from Xizhimen to Dongdan in Beijing, China) are implemented: (1) For the simple network, to test the feasibility of the developed heuristic algorithm, an exact algorithm is first used to solve the model. And the results show that the proposed algorithm is effective and feasible within the allowable error range of 1.5%. (2) For the actual network, 12 groups of numerical experiments are carried out to analyze the optimal schemes under different weight coefficients. The results indicate that the choice of the near-optimal scheme is affected by the weight coefficients, and the results are consistent with the actual situation, Consequently, the proposed model and developed algorithm can promote the planning of multimodal transportation systems in the actual transportation network.
In this paper, the proposed model and developed algorithm are flexible, and the objective function and constraints can be modified. In our future works, we can consider the network accessibility, environmental pollution, and other aspects, and modify the objective function and constraints correspondingly. In addition, the algorithm can be combined with other heuristic algorithms, such as tabu search algorithm, genetic algorithm, and so on, to get the near-optimal solution more quickly.
The study of multimodal transportation DNDP is conductive to guide passengers' multimodal travel behavior and relieve congestion. The practical implications are as follows: (1) By studying the multimodal transportation DNDP, the optimal investment scheme of transportation network planning and construction can be obtained to provide comparison and reference for traffic planning departments, decision-makers, and researchers. (2) In the complicated urban transportation network formed by various modes, grasping the mutual influence between the operation of the multimodal transportation networks as a whole is helpful to improve the level of transportation management measures, and also conducive to giving full play to the advantages of various modes and coordinate the combination of various modes. These will bring great improvement to the transportation congestion, thus reducing vehicle delays and transportation costs.
Certain disadvantages of the proposed methods need to be mentioned here. In this paper, only single OD is considered, and the passenger demand is determined. In future work, we will study the multimodal DNDP of multi-OD under uncertain demand. In addition, this paper only studies five representative travel modes: car, bus, and rail transit, walking, and bike-sharing. However, they can be subdivided into many types in real life. For example, the car model includes private cars, taxis, online car-hailing, etc. And electromobile, motorcycle and so on are also the travel modes that many people choose. Therefore, we can take these into account to improve the proposed model in future work.