# Computational Methods for Calculating Multimodal Multiclass Traffic Network Equilibrium: Simulation Benchmark on a Large-Scale Test Case

Authors: Mostafa Ameli; Jean-Patrick Lebacque; Ludovic Leclercq
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

This study reviews existing computational methods to calculate simulation-based dynamic network equilibrium. We consider a trip-based multimodal approach for the dynamic network loading. Mode and path choices are carried out at the same level; therefore, travel times depend on the travel path and the mode attributes of travelers. This study develops a multiclass model with several parameters per class. Two different categories of algorithms (heuristic and metaheuristic) are considered in order to solve the discrete dynamic traffic assignment (DTA) problem. Finally, we analyze the equilibrium in a large-scale multimodal DTA test case (Lyon 6th + Villeurbanne) in order to investigate the performance of different optimization approaches to solve trip-based DTA. The results show that, in a multimodal and heterogeneous setting, the metaheuristic methods provide better solutions than the heuristic methods in terms of optimality and computation time. These improvements are even more significant than in a homogeneous setting.


---

## Body

## 1. Introduction

Network user equilibrium is the situation where all users’ travel costs in all routes used are equal and lower than those related to the not chosen route [1]. Dynamic traffic assignment (DTA) is a process to assign a given set of trips to a possible path (set of routes and modes) following the network user equilibrium by considering the dynamics of the transportation system. In a large-scale network, when we have numerous paths per origin-destination (OD) pairs and also large numbers of origins and destinations, the problem becomes almost intractable analytically because of multiple flow exchanges at nodes. The simulation-based approach, introduced by [2], is an effective way to address the question of DTA in large-scale networks [3]. Simulation-based DTA models rely on a network performance module called dynamic network loading (DNL). In DNL, the dynamic traffic simulator does the network loading to calculate network states and travel times; then the traffic assignment component determines the route and path flow distribution [4]. This combination, for large-scale networks, allows us to consider sophisticated traffic models in order to incorporate disaggregate information into the estimates of travel costs [5, 6].Traffic simulators can be divided into two classes: flow-based models, which consider the flow of each path, and trip-based models, which define how many travelers take each path. Macroscopic traffic flow models fall into the first category, while microscopic models belong to the second. The macroscopic approach and flow-based models are usually fast in equilibrium calculation as the path flow discipline is more flexible (flows are not necessarily equivalent to traveler units). However, they are less realistic for OD pairs without adding integrality constraints with low demand as vehicles are split into small particles in practice. In this study, we decide to focus on the trip-based approach in which each traveler is reproduced individually. The trip-based approach can be considered as a multiagent system, wherein each trip impacts the system by its route (including path and mode(s)) choice. Microscopic traffic simulators are now widely used for operational studies, and we have chosen to focus on DTA performance for this kind of model. Trip-based DNL attempts to assign particle-discretized time-dependent origin/destination flows in a dynamic network equilibrium framework [7]. The goal of this study is to design and solve a large-scale trip-based traffic network equilibrium model.The user equilibrium situation, theoretically, is the solution of a fixed-point problem [8]. From a practical point of view, it means that each user perceives no delay or extra cost compared to other users with the same origin and destination [9]. This refers to minimizing the total gap between users’ travel cost of each OD pair and the corresponding shortest path (the path(s) with minimum travel cost) travel cost [10]. In the context of DTA, the Method of Successive Average (MSA) algorithm remains by far the most widely used solution method in simulation-based DTA [11]. The simplicity of the MSA implementation and the nonrequirement of derivative information are the main reasons for its widespread use [12]. However, Yang and Jayakrishnan [13] showed that the convergence is not guaranteed to be the user equilibrium solution with the traditional MSA algorithm. Based on this indication, several heuristic approaches are proposed to improve the determination of the step size (e.g., [14–16]). Ameli et al. [17] did a benchmark on most of the heuristic methods in the literature. Here, in addition to heuristic algorithms, we also considered metaheuristic approaches (proposed by [18]) in order to compare different methods to solve the DTA problem.The travel cost is the criterion for each user to select the travel route and mode(s). It is important to consider generalized travel costs in simulation-based DTA in order to include different features of travelers’ behavior. Generally, travelers consider different criteria in selecting their optimal properties of travel, including travel route and travel mode(s). Thus, the total demand of the system can be clustered into different user classes based on the definition of their travel cost function. Many attempts have been made on modeling travelers’ behavior in the traffic network. Most multiclass flow models classify travelers by identifying traveler attributes. In the literature, this classification process identifies eight categories of information about travelers: (1) choice model, (2) cost function, (3) economic attributes, (4) generalized cost function, (5) knowledge level of the network, (6) risk taking attitude, (7) social class, and (8) vehicle class. Moreover, some research uses hybrid classification (see, e.g., [19–21]). In this study, we attempt to define a generic cost function for users in order to represent their behavior accurately. For the convenience of the reader, we discuss all categories and identify the different approaches in Appendix A (Table 1) in order to clarify our contribution compared to the previously mentioned papers.Table 1
List of notations used in this article.
Equilibrium model:WOD pairs, subset of origin× destination nodes, W⊂N×N.ISet of user classes.Pw,τSet of paths for OD pairw in departure time τ.Pw,τi∗Set of shortest paths for OD pairw and class i in departure time interval τ.iIndex of user class,i∈I.τIndex of departure time interval,τ∈T.wIndex of origin-destination (OD) pair,w∈W.pIndex of path,p∈Pw,τ.pi∗Index of shortest path for classi, pi∗∈Pw,τi∗.πw,p,τiNumber of users in classi from OD pair w, assigned to path p in departure time τ.πw,p,τi∗.Optimal number of users in classi from OD pair w, assigned to path p in departure time τ.GCp,ταi,βiGeneralized cost of pathp for user class i in departure time τ.GCw,τi∗Minimum generalized cost of OD pairw for user class i in departure time τ.Optimization algorithms:jOuter-loop iteration index.kInner-loop iteration index.jmaxMaximum number of outer-loop iterations.kmaxMaximum number of inner-loop iterations.πpi,kThe number of users in classi on path p in iteration k.Pi∗The set of shortest paths in iterationk for class i.zjThe path flow distribution of the outer-loop iterationj.zkThe path flow distribution of the inner-loop iterationk.AGapsThe travel cost gap per user of the path flow distributions.TGapsThe total travel cost gap of the path flow distributions.TGapwsThe total travel cost gap of the path flow distribution of OD pairw in solution s.GapαThe threshold for the exploration rate of the inner loop.GapβThe threshold for the convergence rate of the outer loop.Multiclass DTA models are widely investigated in order to take into account the heterogeneity of the users [22]. While the model becomes more accurate for representing the users’ travel cost function, finding the equilibrium becomes more difficult computationally. There are many studies about developing algorithms to find user equilibrium [9]. This study aims to find out which category of solution algorithms is more efficient for the multiclass DTA problem. There are two main categories of solution methods to solve simulation-based DTA with computational methods: heuristic methods based on fixed-point theory and metaheuristic methods that have been recently introduced [18]. In the literature, these methods are applied individually to networks with different sizes and different levels of details about users and network dynamics [23]. In this study, we apply solution methods with the multiclass model to a real test case (Lyon 6th + Villeurbanne: about 17km2) in order to benchmark all the solution methods. The goal is to answer which solution method is the most efficient in generalized and advanced numerical simulation settings.In the next section, the multiclass model is presented. The Wardrop condition for the network equilibrium model is discussed for trip-based DTA inSection 2. The benchmarked algorithms that find the network equilibrium are explained in Section 3. In Section 4, the numerical experiment and the simulator are presented. Section 5 presents the numerical results and discusses them based on several indicators. Finally, we conclude this study in Section 6.


## 2. Dynamic Traffic Assignment Model

### 2.1. Multiattribute Generalized Travel Cost

In order to describe the traveler’s behavior, a generalized cost function that takes into account four types of information mentioned inSection 1: (2), (3), (4), and (7), is used. Typically, the generalized cost GC function will integrate travel time TT and travel cost TC. Also, GC will take into account the class-dependent value of time (VOT). Zhan and Ukkusuri [24] denoted αi for class i and a vector βi of economic and social class parameters. Therefore, we need to define, for class i and path p (p∈P: set of paths), the path generalized cost GCpαi,βi.The vector of parametersβi includes the traveler mode(s) and social attributes. For instance, one class in the network can be students who pay less to use the transit network and who have the specific value of time. The demand of OD pair w (w∈W: set of OD pairs) is fixed and(1)dwi=∑p∈Pwπpi,(2)dwi≥0where i is the user class index i∈1,2,…,I, Pw is a set of paths for w, and πpi is the number of users of class i on path p. After assigning the user to the path between OD pair w, πa (the number of users on link a) is(3)πai=∑p∈Pπpiδap(4)πa=∑i=1Iπai,where a is the link index (a∈A: set of links in the network), πai is the number of users of class i on link a, and δap=1 if path p contains link a and 0 otherwise. The travel time and cost functions depend on the congestion of the link. Therefore, both functions are calculated by the following formulas:(5)TTa=TTaπa,TCaiβi=TCaiβ,πa,and, consequently, the link generalized cost function is calculated as follows:(6)GCaαi,βi=TCaβ,πa+α.TTaπa.Equation (6) calculates the generalized cost of each link for each class of users. Finally, the cost of path p is calculated by the following equation:(7)GCpαi,βi=∑a∈AGCaα,βδap.Now, we can build a DTA model based on the generalized cost defined by equation (7).


### 2.2. Equilibrium Conditions

Consider a networkGN,A with a finite set of nodes N and a finite set of directed links A. The demand is time-dependent. The period of interest (planning horizon) of duration H is discretized into a set of small time intervals indexed by τ (τ∈T=τ0,τ0+σ,τ0+2σ,…,τ0+Mσ and τ0+Mσ=H). σ is the duration of the time intervals. In an interval τ, travel time and traffic conditions do not change. The important notations to introduce the dynamic equilibrium model and optimization algorithms are as follows.For each classi, for all OD pairs w∈W and for all paths p∈Pw, the dynamic traffic network equilibrium condition with fixed travel demands for the aforementioned traffic network equilibrium problem is(8)GCp,ταi,βi−GCw,τi∗≥0,πp,τiGCp,ταi,βi−GCw,τi∗=0,πp,τi≥0.Lu et al. [10] reformulated the problem as a nonlinear problem in order to minimize the gap function. The gap function is defined as the gap between average path travel time and the shortest path travel time. Therefore, the solution of this problem is equivalent to finding the solution to the following variational inequality:(9)∑w∈W∑τ=1T∑i=1I∑p∈Pw,τGCw,p,τi∗πw,p,τi−πw,p,τi∗≥0,∀πw,p,τi,πw,p,τi∗∈C,where C denotes the flow constraints based on dwi. Note that πw,p,τi, a feasible assignment, is the optimal assignment, if Equation (9) holds for all feasible πw,p,τi. In Equation (9), both πw,p,τi and πw,p,τi∗ are decision variables and hence the gap function is a function of both variables. Existence of solution(s) to Equation (9) requires the continuity of the generalized cost functions of paths with respect to the numbers of users on paths. In this case, the unicity of the solution requires the strict monotonicity of the OD paths generalized cost functions with respect to OD path flows (i.e., the number of cars) [25]. In general urban transportation networks, these conditions do not hold [22].


## 2.1. Multiattribute Generalized Travel Cost

In order to describe the traveler’s behavior, a generalized cost function that takes into account four types of information mentioned inSection 1: (2), (3), (4), and (7), is used. Typically, the generalized cost GC function will integrate travel time TT and travel cost TC. Also, GC will take into account the class-dependent value of time (VOT). Zhan and Ukkusuri [24] denoted αi for class i and a vector βi of economic and social class parameters. Therefore, we need to define, for class i and path p (p∈P: set of paths), the path generalized cost GCpαi,βi.The vector of parametersβi includes the traveler mode(s) and social attributes. For instance, one class in the network can be students who pay less to use the transit network and who have the specific value of time. The demand of OD pair w (w∈W: set of OD pairs) is fixed and(1)dwi=∑p∈Pwπpi,(2)dwi≥0where i is the user class index i∈1,2,…,I, Pw is a set of paths for w, and πpi is the number of users of class i on path p. After assigning the user to the path between OD pair w, πa (the number of users on link a) is(3)πai=∑p∈Pπpiδap(4)πa=∑i=1Iπai,where a is the link index (a∈A: set of links in the network), πai is the number of users of class i on link a, and δap=1 if path p contains link a and 0 otherwise. The travel time and cost functions depend on the congestion of the link. Therefore, both functions are calculated by the following formulas:(5)TTa=TTaπa,TCaiβi=TCaiβ,πa,and, consequently, the link generalized cost function is calculated as follows:(6)GCaαi,βi=TCaβ,πa+α.TTaπa.Equation (6) calculates the generalized cost of each link for each class of users. Finally, the cost of path p is calculated by the following equation:(7)GCpαi,βi=∑a∈AGCaα,βδap.Now, we can build a DTA model based on the generalized cost defined by equation (7).


## 2.2. Equilibrium Conditions

Consider a networkGN,A with a finite set of nodes N and a finite set of directed links A. The demand is time-dependent. The period of interest (planning horizon) of duration H is discretized into a set of small time intervals indexed by τ (τ∈T=τ0,τ0+σ,τ0+2σ,…,τ0+Mσ and τ0+Mσ=H). σ is the duration of the time intervals. In an interval τ, travel time and traffic conditions do not change. The important notations to introduce the dynamic equilibrium model and optimization algorithms are as follows.For each classi, for all OD pairs w∈W and for all paths p∈Pw, the dynamic traffic network equilibrium condition with fixed travel demands for the aforementioned traffic network equilibrium problem is(8)GCp,ταi,βi−GCw,τi∗≥0,πp,τiGCp,ταi,βi−GCw,τi∗=0,πp,τi≥0.Lu et al. [10] reformulated the problem as a nonlinear problem in order to minimize the gap function. The gap function is defined as the gap between average path travel time and the shortest path travel time. Therefore, the solution of this problem is equivalent to finding the solution to the following variational inequality:(9)∑w∈W∑τ=1T∑i=1I∑p∈Pw,τGCw,p,τi∗πw,p,τi−πw,p,τi∗≥0,∀πw,p,τi,πw,p,τi∗∈C,where C denotes the flow constraints based on dwi. Note that πw,p,τi, a feasible assignment, is the optimal assignment, if Equation (9) holds for all feasible πw,p,τi. In Equation (9), both πw,p,τi and πw,p,τi∗ are decision variables and hence the gap function is a function of both variables. Existence of solution(s) to Equation (9) requires the continuity of the generalized cost functions of paths with respect to the numbers of users on paths. In this case, the unicity of the solution requires the strict monotonicity of the OD paths generalized cost functions with respect to OD path flows (i.e., the number of cars) [25]. In general urban transportation networks, these conditions do not hold [22].


## 3. Solution Algorithm to Solve Trip-Based Dynamic Traffic Assignment

The development of efficient algorithms for the numerical computation of equilibria is a topic as important as the qualitative analysis of equilibria. The complexity of equilibrium problems, coupled with their increasing scale, is precluding their resolution via closed-form analytics. Note that we will use a trip-based simulation-based algorithm to look for large-scale user equilibrium. As mentioned before, the solution space is then discrete rather than continuous. The optimizer (algorithm) is embedded in the optimization framework of the DTA problem. First, we present the optimization framework, which comes from the literature. Then we focus on the optimizer part, wherein the path assignment is updated by heuristic and metaheuristic methods.Finding an equilibrium through simulation (when no closed-form analytical solution is available) typically involves a solution scheme that relies on an iterative procedure. In large-scale application, this iterative procedure contains three costly steps: traffic simulation (network loading), time-dependent shortest path calculation, and path flow optimization. In order to focus on the path flow optimization, we choose column generation approach.Figure1 presents the optimization framework based on the column generation approach. The framework separates the three costly steps in order to minimize the number of shortest path calculations which is very costly in large-scale traffic networks. In other words, the idea of the column generation approach is based on the generation of the set of paths to which users may be assigned. Classical approach starts from an empty set and augments it in every iteration as required. The column generation approach contains two loops: outer loop and inner loop. The outer loop is responsible for path discovery, while the inner loops implement the path flow optimization for a given path set. Note that the classic approach executes both steps in a single loop. In large-scale network problems, it is extremely costly to keep the data of all the possible paths between each OD pair, so, the column generation approach definitely reduces memory requirements [26]. Lu et al. [10] implement and examine the column generation approach in DTA context and show that it not only reduces memory requirements but also outperforms other algorithms in convergence with a designed swapping algorithm embedded in the inner loop.Figure 1
Simulation-based DTA process: column generation approach.This approach provides us with a framework to compare different path flow optimization algorithms to address which of them is more efficient for the multiclass multimodal DTA problem. Therefore, we are going to focus on the green modules in Figure1. Before presenting the framework in detail, we need to define some solution quality indicators to measure the distance between each solution and optimal solution (UE path flow distribution).The total gap as a quality indicator in time intervalτ is(10)TGapπ,GC∗=∑w∈W∑i=1k∑p∈Pw,τπw,piGCw,pi−GCw,pi∗.The quality of the solution is evaluated by two indicators. The first indicator that calculates the travel cost gap per user in seconds is(11)AGapπ,GC∗=∑w∈W∑τ=1T∑i=1I∑p∈Pw,τπw,p,τiGCw,p,τi−GCw,p,τi∗∑w∈W∑τ=1T∑i=1I∑p∈Pw,τπw,p,τi.The total gap (Equation (10)) and the first quality indicator (Equation (11)) are calculated for each solution with the updated shortest path. The perfect user equilibrium state has zero AGap; thus, the best optimization algorithm obtains minimum AGap. In other words, AGap measures the distance between the algorithms’ solution and UE. The second indicator is the assignment violation. The violation indicator is defined in order to count users assigned on nonoptimal path(s). The following steps calculate the violation indicator:(1)
Calculate the user violation: It is defined by considering the gap of each user. The following functionUVwi defines the user violation:(12)UVwi=1;ifETTwi−TTw∗TTw∗≥0.1,0;o.w.,whereETTi is the experienced travel time of the user i who travels for OD pair w and TTw∗ is the shortest path of OD pair w. In other words, if the gap between the user perceived travel time and shortest path travel time is bigger than 10% of the shortest path travel time, the user is in violation.(2)
Compute the OD violation: The OD pairw is in violation when there are more than 10% of the users on w in violation. Particularly, function ODVw defines the OD violation:(13)ODVw=1;if,∑i∈IwUVwiIw≥0.1,0;o.w.,whereIw is the set of users who travel for OD pair w.(3)
The violation indicator of networkG is the share of ODs which are in violation. Therefore, the second indicator for the quality of the solution is defined as follows:(14)VG=∑w∈WODVwW.As the first indicator evaluates the overall quality of the solution, the assignment violation evaluates the solution’s quality at the OD level. Based onAGap and VG, we can measure the quality of the solution at each iteration. The AGap indicator is used for checking the convergence. The steps of the column generation for simulation-based DTA are detailed in Algorithm 1.Algorithm 1: Simulation-based DTA: equilibrium calculation.
Result: UE path flow distributionInitialization: Set outer loop iteration index to onej=1, zj= All-or-Nothing assignment (i.e., everyone is placed on the shortest path); Execute time-dependant shortest path algorithm and update the path set for all ODs;CalculateAGapzj based on new shortest path(s);Determinejmax and kmax;Whilej<jmaxand (AGapzj≥Gapαor There is at least one new shortest path for one OD) doSet inner loop iteration index to onek=1;Setzk=zj and AGapz0=AGapzj+2×Gapβ;Whilek<kmaxandAGapzk−AGapzk−1>GapβdoExecuteoptimization algorithm to calculate zk+1;Execute traffic simulation withzk+1;Identify the shortest path for all ODs based on the simulation results;UpdateAGapzk+1 based on updated shortest path(s);Keep the best solution of the current inner loop;Setk=k+1;endProduce the best solution of the inner loop aszk;Setzj+1=zk;Execute time-dependant shortest path algorithm based onzj+1;Update the path set for all ODs;CalculateAGapzj+1 based on new/updated shortest paths;Setj=j+1;end.The framework has two convergence tests to stop the calculation process in order to save computation time. In both steps, first, the iteration index is checked with the maximum number of loops (jmax and kmax) that we fix for each loop. Suppose that it is not reached (j<jmax or k<kmax); then we check the solution’s quality at the current loop, but there are three differences between the outer-loop convergence test and inner-loop convergence test. First, in the outer loop, we calculate the quality indicator (Equation (11)) based on the new shortest path(s), but the path set is fixed during the inner-loop operations, and only the shortest path is updated in each iteration. Second, in the inner loop, we calculate the relative total gap (i.e., the difference of the total gap between two iterations) to check the improvement of the process and compare it with Gapβ. Gapβ is a threshold for checking the exploration rate of the optimization algorithm in the inner loop. But, in the outer loop, the goal is to monitor the quality of the solution. Therefore, AGap is calculated and compared with the fixed value Gapα. If we have a solution with AGap equal to or less than threshold Gapα, we terminate the algorithm. Third, in the outer loop, we have one more condition other than checking the solution’s quality. The condition is the shortest path checking. If we have no new shortest path per OD and good quality of the solution, then the process is converged. It means we will continue the process when the quality of the solution is not good until we no longer have a new shortest path.The optimization algorithm is embedded in the inner loop of the process (Figure1) which updates the assignment plan in each iteration. The heuristic and metaheuristic algorithms are defined in this step.


### 3.1. Heuristic Methods

Generally, the heuristic methods are path-swapping descent direction methods that decompose the problem into many (origin-destination, departure time interval, and user class) subproblems. They solve each of the subproblems by adjusting time-varying OD flows between all nonshortest paths and the shortest path(s). We did a cross-comparison of heuristic algorithms in [17] with homogeneous users. Here we present the best algorithm for a large scale, which may have a good performance in multiclass and multimodal contexts.


#### 3.1.1. Method of Successive Average with Ranking Method

The first method is an extension of the Method of Successive Average (MSA). The MSA algorithm is the classical method to solve the traffic assignment problem [27] and it is widely used in theory and application for DTA problems [28]. Equation (15) presents the swapping formula of MSA algorithm [29]:(15)zk=kk+1zk−1+1k+1yk.where yk is extracted from the auxiliary path assignments obtained by the all-or-nothing procedure; that is, everyone is placed on the shortest path. Note that the step size in this method is 1/k+1. It means that, during iteration k for this method, a fraction 1/k+1 of users on each nonshortest path swaps to the one of the shortest path in the current iteration. MSA ranking adds a ranking method to choose the users from nonshortest paths to the shortest path(s). In other words, it ranks all users by their generalized cost and swaps the users with higher GC to the shortest path(s). The step size is the same as the MSA method equal to 1/k+1 of all users on nonshortest paths.


#### 3.1.2. Gap-Based Method

To enable the deployment of large-scale (simulation-based) DTA models, [10] proposes a heuristic approach, which we called “gap-based” method in this study. This method updates the path flow in optimization iteration k based on gap between the path p generalized cost of user class iGCpi and that OD shortest path generalized cost of user class iGCpi∗:(16)πpi,k+1=max0,πpi,k.1−ρkGCpi−GCpi∗GCpi;ifp∉Pi∗πpi,k+1=πpi,k+∑p∈Pw−Pi∗πpi,k.ρkGCpi−GCpi∗/GCpiPi∗;ifp∈Pi∗,ρk=1k+1;ifj=0,1;o.w.,where j denotes the outer-loop iteration index, k denotes the inner-loop iteration index, and Pi∗ is the number of the shortest paths in iteration k. This path assignment updating scheme implies a natural path flow adjustment mechanism: flows on the nonshortest paths are moved to the shortest path and the volume that moved out from a nonshortest path is proportional to the gap between the nonshortest path and shortest path over nonshortest path cost.MSA ranking and gap-based method have two main drawbacks which are also the same for all heuristic methods that come from the flow-based traffic assignment:Using the step size guarantees the algorithm to converge but decreases the flexibility of the method to swap more users to the shortest path(s). Moreover, the final solution is obtained based on the initialization (all-or-nothing) because the algorithm converges (by decreasing the step size) after a while even we use a not good initial assignment which is not really obvious in the large scale.Many of the methods in the literature (see, e.g., [8, 11, 12, 15, 16, 28, 30–32]) are using the average travel time of each path to calculate the indicators and to distinguish the shortest path(s) and other paths per OD. This means they do not consider users’ experienced travel time directly. Moreover, the optimizer is applied for each OD and swap users based on the different values of the paths’ average travel time. Note that users in the same departure time interval and on the same path can perceive different travel time values because of the network dynamic evolution.These drawbacks show that we always look for the local optimum based on initialization and our step size. Here we present a probabilistic approach to help the heuristic optimization framework in order to find the global optimum for trip-based traffic assignment.


#### 3.1.3. Probabilistic Method

The probabilistic method is a trip-based approach. This method has no step size for optimization [17]. For each user, the probability of swap to the shortest path(s) is calculated by the following formula:(17)Swap indicator=EGCp−GCpi∗GCp,where EGCp is the experienced generalized cost of the user. Indeed, the path-swapping approach is applied to every user, not paths (number of users on path p). Afterwards, we use a random number uniformly distributed between 0 and 1 for taking a decision about swapping the user or not. If the random number is lower than the swap indicator, the user swaps to the shortest path; otherwise, the user keeps his/her path for the current iteration. The probabilistic process enables the algorithm to explore more the solution space without any step size. This heuristic method converges intuitively for a large-scale network because the swap indicator decreases by swapping based on random numbers. Note that, at the end of each outer loop, the final assignment plan is the best assignment in terms of average gap between different inner-loops solutions. It is obvious that, because of the swap indicator, the algorithm converges by random numbers to the solution with a minimum number of possible swaps, which indicates the small value of gap for each user. In other words, with every swapping decision, users try to optimize the GC in order to reduce the swap indicator (the possibility of swapping) for the next inner loop.


#### 3.1.4. Hybrid Methods

Verbas et al. [33] applied an extension of the gap-based method to transit assignment. The method uses Equation (16) to determine the number of swapping users for each path and then the probabilistic method is applied to choose users from each path for swapping. This method is called “gap-based prob.” We examined this method in DTA test cases, and it fully dominates the gap-based method in different monoclass test cases [17]. Therefore, we consider this method in our numerical experiments.The second hybrid method is similar to the probabilistic method with a different swap indicator:(18)Swap indicator=1k+1EGCp−GCpi∗GCp.This method is called “step-size prob.” in this study. It has a good performance in large-scale and congested traffic networks [20]. We developed a smart step-size selection procedure for this method in which we specialize the step size for each OD pair and decrease the step size when the OD gap is not improved compared to the previous iteration [17]. The method is called “smart step-size prob.” All the methods that are presented in this section will apply to multiclass multimodal DTA problem in order to evaluate their performance and also compare them with metaheuristic methods.


### 3.2. Metaheuristic Methods

As discussed in the previous section, in iterationk, the classical (heuristic) algorithms try to improve the path flow distribution by swapping a fraction 1/k+1 of users to the shortest path(s) from each nonshortest path. Then, one simulation is launched based on the updated path flow distribution. We have to wait until the simulation run is finished to know the new link travel times and adjust the path flow distribution to be tested in the next iteration accordingly. This is the serial process of the classical algorithms, which limits the solution space exploration and computational process. We developed two metaheuristic algorithms in [18] in order to overcome the problem of the serial approach. Two parallelized metaheuristic approaches are designed: the first derived from the simulated annealing framework and the second derived from that of the genetic algorithm.


#### 3.2.1. Simulated Annealing Algorithm

The algorithm starts with an initial solution generated randomly; that is, users choose their path from the OD path set randomly. For solutionz, the total gap TGap. between the users’ travel cost and the shortest path travel cost (Equation (10)) is considered as the energy of the solution. A set of neighbor solutions are generated with respect to the current one based on the temperature T of the current iteration. A neighbor solution is defined as a candidate to replace the current solution. The current phase of the iteration depends on the temperature of the process. Inspired by the physics of matter, this study distinguishes three different methods to generate a neighbor solution, gas, liquid, and solid; these methods represent the states of matter in nature. When the temperature is high (T>λ where λ denotes the boiling temperature), the gas method is applied. When running the SA algorithm, by decreasing the temperature, the algorithm enters the liquid phase (λ>T>λ′ where λ′ denotes the melting temperature) and then the liquid method is applied. When the temperature is quite low T<λ′, the solid method is applied. The algorithm is presented in Algorithm 2.Algorithm 2: Simulated annealing algorithm.
Result: Inner loop path flow distribution of iteration kInitialization: SetT=T0, where T0 denotes the initial temperature;Set inner loop iteration index to onek=1 and AGapz0=AGapz−1=2AGapzj;Whilek<kmaxandAGapzk≠AGapzk−1≠AGapzk−2doifT>λthenGas Phase:Generate a solution candidate by Randomization method.Randomization: users choose their path randomly from the path set;else ifT>λ′thenLiquid Phase:Generate the first solution candidate by Randomization method; Generate the second and third solution candidates by applying MSA and Gap-based methods tozi;elseSolid Phase:Generate the first and second solution candidates by applying MSA and Gap-based methods tozi;endExecute traffic simulation in parallel for all candidates;Identify the shortest path for all ODs based on the simulation results for each candidate;UpdateAGap. based on updated shortest path(s) for all candidates;Calculate the solution acceptance probability for each candidate by the following equation:PSaccepted′=1=e−∇E/Twheres′ denotes a candidate, Saccepted′ denotes the binary decision variable and ∇E=TGaps′−TGaps;Take decision about accepting each candidate to determinezk+1;Decrease the temperature by the following formula:T=T0/lnk+1;Setk=k+1;endProduce the final solution for the outer loop.Algorithm 3: Genetic algorithm.
Result: Inner loop path flow distribution of iteration kInitialization: Generate the initial solution population based on Randomization;Ifj>1thenAdd the solution comes from the outer loop.Set inner loop iteration index to onek=1;Mp0=Mp−1=0, where Mpk denotes the mean average gap of the population of iteration kpk;whilek<kmaxandMpk≠Mpk−1≠Mpk−2do1. GA-operators: Apply the GA operators to each pair of selected solutions.1.1. Selection: Random selection based on the crossover rate Cr and population size PS to compute the number of selected solutions for the crossover process: SS=PS×Cr (Note that based on [36], we set Cr=0.5 and the PS is set to 10 individuals in this study.); 1.2. Crossover: Apply a non-uniform crossover by using a bit-vector mask method [37]. We select two different solutions (parents) from the set of selected solutions. Then apply the crossover between each pair of solutions. Two new solutions are generated.; 1.3. Mutation: Apply the mutation operator by replacing one OD assignment of the solution by another corresponding OD assignment from another solution. The possibility of the mutation for one OD assignment is calculated based on the quality of the OD assignment:MPw=TGapws/TGapswhereTGapws denotes the total gap of the OD assignment w in solution s (note that in application of each mutation operator, two new solutions are generated).;2. GA Inner: Apply the GA Inner operators to every selected solutionsAdaptive random method [18];Execute traffic simulation in parallel for all new solutions obtained from the previous steps;Identify the shortest path for all ODs based on the simulation results for every new solutions;UpdateAGap. based on updated shortest path(s);Apply the selection operator to keep the size of population fix by a replacement strategy.;Setk=k+1;endProduce the final solution (a solution with the minimumAGap in the population pk) for the outer loop.


#### 3.2.2. Genetic Algorithm

The genetic algorithm (GA) (first formalized by [34]) is inspired by natural selection and genetic variation. The genetic algorithm includes three main operators: selection, crossover, and mutation. The description of the operators is well presented in [35]. The GA framework is adapted to the traffic assignment problem in [18]. The algorithm is designed as a two-layer GA process to search solution space by changing the path flows in GA Inner and overcome the drawback of OD impact by considering different combinations of OD assignments in GA-operators. In other words, the heuristic algorithms plus a random method are applied in GA Inner, and the GA operators in one upper level generate different combinations of OD. GA Inner is a new layer in the optimization process compared to the classical GA. It has been added to account for correlations between OD assignments. Moreover, GA considers a set of solutions instead of a single solution at every iteration. The steps of the GA are detailed in Algorithm 3.


## 3.1. Heuristic Methods

Generally, the heuristic methods are path-swapping descent direction methods that decompose the problem into many (origin-destination, departure time interval, and user class) subproblems. They solve each of the subproblems by adjusting time-varying OD flows between all nonshortest paths and the shortest path(s). We did a cross-comparison of heuristic algorithms in [17] with homogeneous users. Here we present the best algorithm for a large scale, which may have a good performance in multiclass and multimodal contexts.


### 3.1.1. Method of Successive Average with Ranking Method

The first method is an extension of the Method of Successive Average (MSA). The MSA algorithm is the classical method to solve the traffic assignment problem [27] and it is widely used in theory and application for DTA problems [28]. Equation (15) presents the swapping formula of MSA algorithm [29]:(15)zk=kk+1zk−1+1k+1yk.where yk is extracted from the auxiliary path assignments obtained by the all-or-nothing procedure; that is, everyone is placed on the shortest path. Note that the step size in this method is 1/k+1. It means that, during iteration k for this method, a fraction 1/k+1 of users on each nonshortest path swaps to the one of the shortest path in the current iteration. MSA ranking adds a ranking method to choose the users from nonshortest paths to the shortest path(s). In other words, it ranks all users by their generalized cost and swaps the users with higher GC to the shortest path(s). The step size is the same as the MSA method equal to 1/k+1 of all users on nonshortest paths.


### 3.1.2. Gap-Based Method

To enable the deployment of large-scale (simulation-based) DTA models, [10] proposes a heuristic approach, which we called “gap-based” method in this study. This method updates the path flow in optimization iteration k based on gap between the path p generalized cost of user class iGCpi and that OD shortest path generalized cost of user class iGCpi∗:(16)πpi,k+1=max0,πpi,k.1−ρkGCpi−GCpi∗GCpi;ifp∉Pi∗πpi,k+1=πpi,k+∑p∈Pw−Pi∗πpi,k.ρkGCpi−GCpi∗/GCpiPi∗;ifp∈Pi∗,ρk=1k+1;ifj=0,1;o.w.,where j denotes the outer-loop iteration index, k denotes the inner-loop iteration index, and Pi∗ is the number of the shortest paths in iteration k. This path assignment updating scheme implies a natural path flow adjustment mechanism: flows on the nonshortest paths are moved to the shortest path and the volume that moved out from a nonshortest path is proportional to the gap between the nonshortest path and shortest path over nonshortest path cost.MSA ranking and gap-based method have two main drawbacks which are also the same for all heuristic methods that come from the flow-based traffic assignment:Using the step size guarantees the algorithm to converge but decreases the flexibility of the method to swap more users to the shortest path(s). Moreover, the final solution is obtained based on the initialization (all-or-nothing) because the algorithm converges (by decreasing the step size) after a while even we use a not good initial assignment which is not really obvious in the large scale.Many of the methods in the literature (see, e.g., [8, 11, 12, 15, 16, 28, 30–32]) are using the average travel time of each path to calculate the indicators and to distinguish the shortest path(s) and other paths per OD. This means they do not consider users’ experienced travel time directly. Moreover, the optimizer is applied for each OD and swap users based on the different values of the paths’ average travel time. Note that users in the same departure time interval and on the same path can perceive different travel time values because of the network dynamic evolution.These drawbacks show that we always look for the local optimum based on initialization and our step size. Here we present a probabilistic approach to help the heuristic optimization framework in order to find the global optimum for trip-based traffic assignment.


### 3.1.3. Probabilistic Method

The probabilistic method is a trip-based approach. This method has no step size for optimization [17]. For each user, the probability of swap to the shortest path(s) is calculated by the following formula:(17)Swap indicator=EGCp−GCpi∗GCp,where EGCp is the experienced generalized cost of the user. Indeed, the path-swapping approach is applied to every user, not paths (number of users on path p). Afterwards, we use a random number uniformly distributed between 0 and 1 for taking a decision about swapping the user or not. If the random number is lower than the swap indicator, the user swaps to the shortest path; otherwise, the user keeps his/her path for the current iteration. The probabilistic process enables the algorithm to explore more the solution space without any step size. This heuristic method converges intuitively for a large-scale network because the swap indicator decreases by swapping based on random numbers. Note that, at the end of each outer loop, the final assignment plan is the best assignment in terms of average gap between different inner-loops solutions. It is obvious that, because of the swap indicator, the algorithm converges by random numbers to the solution with a minimum number of possible swaps, which indicates the small value of gap for each user. In other words, with every swapping decision, users try to optimize the GC in order to reduce the swap indicator (the possibility of swapping) for the next inner loop.


### 3.1.4. Hybrid Methods

Verbas et al. [33] applied an extension of the gap-based method to transit assignment. The method uses Equation (16) to determine the number of swapping users for each path and then the probabilistic method is applied to choose users from each path for swapping. This method is called “gap-based prob.” We examined this method in DTA test cases, and it fully dominates the gap-based method in different monoclass test cases [17]. Therefore, we consider this method in our numerical experiments.The second hybrid method is similar to the probabilistic method with a different swap indicator:(18)Swap indicator=1k+1EGCp−GCpi∗GCp.This method is called “step-size prob.” in this study. It has a good performance in large-scale and congested traffic networks [20]. We developed a smart step-size selection procedure for this method in which we specialize the step size for each OD pair and decrease the step size when the OD gap is not improved compared to the previous iteration [17]. The method is called “smart step-size prob.” All the methods that are presented in this section will apply to multiclass multimodal DTA problem in order to evaluate their performance and also compare them with metaheuristic methods.


## 3.1.1. Method of Successive Average with Ranking Method

The first method is an extension of the Method of Successive Average (MSA). The MSA algorithm is the classical method to solve the traffic assignment problem [27] and it is widely used in theory and application for DTA problems [28]. Equation (15) presents the swapping formula of MSA algorithm [29]:(15)zk=kk+1zk−1+1k+1yk.where yk is extracted from the auxiliary path assignments obtained by the all-or-nothing procedure; that is, everyone is placed on the shortest path. Note that the step size in this method is 1/k+1. It means that, during iteration k for this method, a fraction 1/k+1 of users on each nonshortest path swaps to the one of the shortest path in the current iteration. MSA ranking adds a ranking method to choose the users from nonshortest paths to the shortest path(s). In other words, it ranks all users by their generalized cost and swaps the users with higher GC to the shortest path(s). The step size is the same as the MSA method equal to 1/k+1 of all users on nonshortest paths.


## 3.1.2. Gap-Based Method

To enable the deployment of large-scale (simulation-based) DTA models, [10] proposes a heuristic approach, which we called “gap-based” method in this study. This method updates the path flow in optimization iteration k based on gap between the path p generalized cost of user class iGCpi and that OD shortest path generalized cost of user class iGCpi∗:(16)πpi,k+1=max0,πpi,k.1−ρkGCpi−GCpi∗GCpi;ifp∉Pi∗πpi,k+1=πpi,k+∑p∈Pw−Pi∗πpi,k.ρkGCpi−GCpi∗/GCpiPi∗;ifp∈Pi∗,ρk=1k+1;ifj=0,1;o.w.,where j denotes the outer-loop iteration index, k denotes the inner-loop iteration index, and Pi∗ is the number of the shortest paths in iteration k. This path assignment updating scheme implies a natural path flow adjustment mechanism: flows on the nonshortest paths are moved to the shortest path and the volume that moved out from a nonshortest path is proportional to the gap between the nonshortest path and shortest path over nonshortest path cost.MSA ranking and gap-based method have two main drawbacks which are also the same for all heuristic methods that come from the flow-based traffic assignment:Using the step size guarantees the algorithm to converge but decreases the flexibility of the method to swap more users to the shortest path(s). Moreover, the final solution is obtained based on the initialization (all-or-nothing) because the algorithm converges (by decreasing the step size) after a while even we use a not good initial assignment which is not really obvious in the large scale.Many of the methods in the literature (see, e.g., [8, 11, 12, 15, 16, 28, 30–32]) are using the average travel time of each path to calculate the indicators and to distinguish the shortest path(s) and other paths per OD. This means they do not consider users’ experienced travel time directly. Moreover, the optimizer is applied for each OD and swap users based on the different values of the paths’ average travel time. Note that users in the same departure time interval and on the same path can perceive different travel time values because of the network dynamic evolution.These drawbacks show that we always look for the local optimum based on initialization and our step size. Here we present a probabilistic approach to help the heuristic optimization framework in order to find the global optimum for trip-based traffic assignment.


## 3.1.3. Probabilistic Method

The probabilistic method is a trip-based approach. This method has no step size for optimization [17]. For each user, the probability of swap to the shortest path(s) is calculated by the following formula:(17)Swap indicator=EGCp−GCpi∗GCp,where EGCp is the experienced generalized cost of the user. Indeed, the path-swapping approach is applied to every user, not paths (number of users on path p). Afterwards, we use a random number uniformly distributed between 0 and 1 for taking a decision about swapping the user or not. If the random number is lower than the swap indicator, the user swaps to the shortest path; otherwise, the user keeps his/her path for the current iteration. The probabilistic process enables the algorithm to explore more the solution space without any step size. This heuristic method converges intuitively for a large-scale network because the swap indicator decreases by swapping based on random numbers. Note that, at the end of each outer loop, the final assignment plan is the best assignment in terms of average gap between different inner-loops solutions. It is obvious that, because of the swap indicator, the algorithm converges by random numbers to the solution with a minimum number of possible swaps, which indicates the small value of gap for each user. In other words, with every swapping decision, users try to optimize the GC in order to reduce the swap indicator (the possibility of swapping) for the next inner loop.


## 3.1.4. Hybrid Methods

Verbas et al. [33] applied an extension of the gap-based method to transit assignment. The method uses Equation (16) to determine the number of swapping users for each path and then the probabilistic method is applied to choose users from each path for swapping. This method is called “gap-based prob.” We examined this method in DTA test cases, and it fully dominates the gap-based method in different monoclass test cases [17]. Therefore, we consider this method in our numerical experiments.The second hybrid method is similar to the probabilistic method with a different swap indicator:(18)Swap indicator=1k+1EGCp−GCpi∗GCp.This method is called “step-size prob.” in this study. It has a good performance in large-scale and congested traffic networks [20]. We developed a smart step-size selection procedure for this method in which we specialize the step size for each OD pair and decrease the step size when the OD gap is not improved compared to the previous iteration [17]. The method is called “smart step-size prob.” All the methods that are presented in this section will apply to multiclass multimodal DTA problem in order to evaluate their performance and also compare them with metaheuristic methods.


## 3.2. Metaheuristic Methods

As discussed in the previous section, in iterationk, the classical (heuristic) algorithms try to improve the path flow distribution by swapping a fraction 1/k+1 of users to the shortest path(s) from each nonshortest path. Then, one simulation is launched based on the updated path flow distribution. We have to wait until the simulation run is finished to know the new link travel times and adjust the path flow distribution to be tested in the next iteration accordingly. This is the serial process of the classical algorithms, which limits the solution space exploration and computational process. We developed two metaheuristic algorithms in [18] in order to overcome the problem of the serial approach. Two parallelized metaheuristic approaches are designed: the first derived from the simulated annealing framework and the second derived from that of the genetic algorithm.


### 3.2.1. Simulated Annealing Algorithm

The algorithm starts with an initial solution generated randomly; that is, users choose their path from the OD path set randomly. For solutionz, the total gap TGap. between the users’ travel cost and the shortest path travel cost (Equation (10)) is considered as the energy of the solution. A set of neighbor solutions are generated with respect to the current one based on the temperature T of the current iteration. A neighbor solution is defined as a candidate to replace the current solution. The current phase of the iteration depends on the temperature of the process. Inspired by the physics of matter, this study distinguishes three different methods to generate a neighbor solution, gas, liquid, and solid; these methods represent the states of matter in nature. When the temperature is high (T>λ where λ denotes the boiling temperature), the gas method is applied. When running the SA algorithm, by decreasing the temperature, the algorithm enters the liquid phase (λ>T>λ′ where λ′ denotes the melting temperature) and then the liquid method is applied. When the temperature is quite low T<λ′, the solid method is applied. The algorithm is presented in Algorithm 2.Algorithm 2: Simulated annealing algorithm.
Result: Inner loop path flow distribution of iteration kInitialization: SetT=T0, where T0 denotes the initial temperature;Set inner loop iteration index to onek=1 and AGapz0=AGapz−1=2AGapzj;Whilek<kmaxandAGapzk≠AGapzk−1≠AGapzk−2doifT>λthenGas Phase:Generate a solution candidate by Randomization method.Randomization: users choose their path randomly from the path set;else ifT>λ′thenLiquid Phase:Generate the first solution candidate by Randomization method; Generate the second and third solution candidates by applying MSA and Gap-based methods tozi;elseSolid Phase:Generate the first and second solution candidates by applying MSA and Gap-based methods tozi;endExecute traffic simulation in parallel for all candidates;Identify the shortest path for all ODs based on the simulation results for each candidate;UpdateAGap. based on updated shortest path(s) for all candidates;Calculate the solution acceptance probability for each candidate by the following equation:PSaccepted′=1=e−∇E/Twheres′ denotes a candidate, Saccepted′ denotes the binary decision variable and ∇E=TGaps′−TGaps;Take decision about accepting each candidate to determinezk+1;Decrease the temperature by the following formula:T=T0/lnk+1;Setk=k+1;endProduce the final solution for the outer loop.Algorithm 3: Genetic algorithm.
Result: Inner loop path flow distribution of iteration kInitialization: Generate the initial solution population based on Randomization;Ifj>1thenAdd the solution comes from the outer loop.Set inner loop iteration index to onek=1;Mp0=Mp−1=0, where Mpk denotes the mean average gap of the population of iteration kpk;whilek<kmaxandMpk≠Mpk−1≠Mpk−2do1. GA-operators: Apply the GA operators to each pair of selected solutions.1.1. Selection: Random selection based on the crossover rate Cr and population size PS to compute the number of selected solutions for the crossover process: SS=PS×Cr (Note that based on [36], we set Cr=0.5 and the PS is set to 10 individuals in this study.); 1.2. Crossover: Apply a non-uniform crossover by using a bit-vector mask method [37]. We select two different solutions (parents) from the set of selected solutions. Then apply the crossover between each pair of solutions. Two new solutions are generated.; 1.3. Mutation: Apply the mutation operator by replacing one OD assignment of the solution by another corresponding OD assignment from another solution. The possibility of the mutation for one OD assignment is calculated based on the quality of the OD assignment:MPw=TGapws/TGapswhereTGapws denotes the total gap of the OD assignment w in solution s (note that in application of each mutation operator, two new solutions are generated).;2. GA Inner: Apply the GA Inner operators to every selected solutionsAdaptive random method [18];Execute traffic simulation in parallel for all new solutions obtained from the previous steps;Identify the shortest path for all ODs based on the simulation results for every new solutions;UpdateAGap. based on updated shortest path(s);Apply the selection operator to keep the size of population fix by a replacement strategy.;Setk=k+1;endProduce the final solution (a solution with the minimumAGap in the population pk) for the outer loop.


### 3.2.2. Genetic Algorithm

The genetic algorithm (GA) (first formalized by [34]) is inspired by natural selection and genetic variation. The genetic algorithm includes three main operators: selection, crossover, and mutation. The description of the operators is well presented in [35]. The GA framework is adapted to the traffic assignment problem in [18]. The algorithm is designed as a two-layer GA process to search solution space by changing the path flows in GA Inner and overcome the drawback of OD impact by considering different combinations of OD assignments in GA-operators. In other words, the heuristic algorithms plus a random method are applied in GA Inner, and the GA operators in one upper level generate different combinations of OD. GA Inner is a new layer in the optimization process compared to the classical GA. It has been added to account for correlations between OD assignments. Moreover, GA considers a set of solutions instead of a single solution at every iteration. The steps of the GA are detailed in Algorithm 3.


## 3.2.1. Simulated Annealing Algorithm

The algorithm starts with an initial solution generated randomly; that is, users choose their path from the OD path set randomly. For solutionz, the total gap TGap. between the users’ travel cost and the shortest path travel cost (Equation (10)) is considered as the energy of the solution. A set of neighbor solutions are generated with respect to the current one based on the temperature T of the current iteration. A neighbor solution is defined as a candidate to replace the current solution. The current phase of the iteration depends on the temperature of the process. Inspired by the physics of matter, this study distinguishes three different methods to generate a neighbor solution, gas, liquid, and solid; these methods represent the states of matter in nature. When the temperature is high (T>λ where λ denotes the boiling temperature), the gas method is applied. When running the SA algorithm, by decreasing the temperature, the algorithm enters the liquid phase (λ>T>λ′ where λ′ denotes the melting temperature) and then the liquid method is applied. When the temperature is quite low T<λ′, the solid method is applied. The algorithm is presented in Algorithm 2.Algorithm 2: Simulated annealing algorithm.
Result: Inner loop path flow distribution of iteration kInitialization: SetT=T0, where T0 denotes the initial temperature;Set inner loop iteration index to onek=1 and AGapz0=AGapz−1=2AGapzj;Whilek<kmaxandAGapzk≠AGapzk−1≠AGapzk−2doifT>λthenGas Phase:Generate a solution candidate by Randomization method.Randomization: users choose their path randomly from the path set;else ifT>λ′thenLiquid Phase:Generate the first solution candidate by Randomization method; Generate the second and third solution candidates by applying MSA and Gap-based methods tozi;elseSolid Phase:Generate the first and second solution candidates by applying MSA and Gap-based methods tozi;endExecute traffic simulation in parallel for all candidates;Identify the shortest path for all ODs based on the simulation results for each candidate;UpdateAGap. based on updated shortest path(s) for all candidates;Calculate the solution acceptance probability for each candidate by the following equation:PSaccepted′=1=e−∇E/Twheres′ denotes a candidate, Saccepted′ denotes the binary decision variable and ∇E=TGaps′−TGaps;Take decision about accepting each candidate to determinezk+1;Decrease the temperature by the following formula:T=T0/lnk+1;Setk=k+1;endProduce the final solution for the outer loop.Algorithm 3: Genetic algorithm.
Result: Inner loop path flow distribution of iteration kInitialization: Generate the initial solution population based on Randomization;Ifj>1thenAdd the solution comes from the outer loop.Set inner loop iteration index to onek=1;Mp0=Mp−1=0, where Mpk denotes the mean average gap of the population of iteration kpk;whilek<kmaxandMpk≠Mpk−1≠Mpk−2do1. GA-operators: Apply the GA operators to each pair of selected solutions.1.1. Selection: Random selection based on the crossover rate Cr and population size PS to compute the number of selected solutions for the crossover process: SS=PS×Cr (Note that based on [36], we set Cr=0.5 and the PS is set to 10 individuals in this study.); 1.2. Crossover: Apply a non-uniform crossover by using a bit-vector mask method [37]. We select two different solutions (parents) from the set of selected solutions. Then apply the crossover between each pair of solutions. Two new solutions are generated.; 1.3. Mutation: Apply the mutation operator by replacing one OD assignment of the solution by another corresponding OD assignment from another solution. The possibility of the mutation for one OD assignment is calculated based on the quality of the OD assignment:MPw=TGapws/TGapswhereTGapws denotes the total gap of the OD assignment w in solution s (note that in application of each mutation operator, two new solutions are generated).;2. GA Inner: Apply the GA Inner operators to every selected solutionsAdaptive random method [18];Execute traffic simulation in parallel for all new solutions obtained from the previous steps;Identify the shortest path for all ODs based on the simulation results for every new solutions;UpdateAGap. based on updated shortest path(s);Apply the selection operator to keep the size of population fix by a replacement strategy.;Setk=k+1;endProduce the final solution (a solution with the minimumAGap in the population pk) for the outer loop.


## 3.2.2. Genetic Algorithm

The genetic algorithm (GA) (first formalized by [34]) is inspired by natural selection and genetic variation. The genetic algorithm includes three main operators: selection, crossover, and mutation. The description of the operators is well presented in [35]. The GA framework is adapted to the traffic assignment problem in [18]. The algorithm is designed as a two-layer GA process to search solution space by changing the path flows in GA Inner and overcome the drawback of OD impact by considering different combinations of OD assignments in GA-operators. In other words, the heuristic algorithms plus a random method are applied in GA Inner, and the GA operators in one upper level generate different combinations of OD. GA Inner is a new layer in the optimization process compared to the classical GA. It has been added to account for correlations between OD assignments. Moreover, GA considers a set of solutions instead of a single solution at every iteration. The steps of the GA are detailed in Algorithm 3.


## 4. Numerical Experiments

The multiclass model is applied to a large-scale multimodal network of Lyon 6th + Villeurbanne. This network that covers 25km2 and includes 1,883 nodes, 3,383 links, 94 origins, and 227 destinations is located in Lyon Metropolis and is described in Figure 2. All the signalized intersections in the real field have been implemented in the simulator with their actual signal timing. In this work, we use Symuvia (an open source simulator (https://github.com/Ifsttar/Open-SymuVia)) as a trip-based simulator for calculating the needed variables in the network. Symuvia is a microscopic simulator based on a Lagrangian discretization of the LWR model [38]. We get the simulation time-step to 1 second and collect the link traffic information (travel times) every 1 minute. Note that the dynamic simulator uses Equation (6) to calculate the generalized travel cost on a link at each time interval. There are buses (with/without a dedicated line), trams, and metro lines considering public transportations. In the simulation process, all these patterns, in addition to the private car and walking option, are simulated. For instance, for using a metro line, a user walks (walking time) to the station and then waits for the next train (waiting time). When the train arrives at the station, if the train has free capacity, the user can take the train (boarding time) and then take the next patterns until arriving at the destination point. All of these actions are considered for the travel time calculation. The users’ routes are determined by the DTA model and the rolling horizon technique [39], which determine the path flow distribution based on a prediction period of 30 minutes and an assignment period of 15 minutes. The multimodal time-dependent shortest path algorithm is implemented based on the multiorigin Dijkstra algorithm [40] with considering hyperpaths based on [41] in order to find minimum cost path(s) in terms of travel time.Figure 2
Lyon 6th + Villeurbanne: (a) mapping data ©Google 2020; (b) the traffic network considered in Symuvia.
(a)(b)The network is loaded with 54,190 travelers of all ODs with given departure times in order to represent the 2.5 morning peak hours of the network between 6 : 30 and 9 : 00. The demand profile comes from the study of [42]. Note that, here, we fixed the departure time in order to focus on travelers’ path choices. The main reason is that this study aims to compare the performances of DTA solvers. By fixing the departure time, the input of the problem will be the same for all algorithms. Therefore, the solution space is fixed and given for the route (path and mode) choice problem, and we can compare the performances of different solution algorithms.Figure3(a) presents the demand profile of the numerical experiment. The figure presents the inflow rate of the trips to the traffic network. The total length of the Lyon 6th + Villeurbanne network is 496.73 kilometers. Note that each OD has a time-dependent demand profile based on the real data of Lyon city [42]. The dynamic loading represents the saturation state of the network. To show a quick and synthetic overview of the network state, we plot the Macroscopic Fundamental Diagrams (MFD) of the network in Figure 3(b). This diagram represents the overall evolution of the traffic conditions in the network by plotting the total travel distance versus the total travel time [43]. Note that total travel distance is proportional to the mean network flow, while total travel time is equivalent to vehicle accumulation.Figure 3
The demand scenario of Lyon 6th + Villeurbanne. (a) The demand profile. (b) The macroscopic fundamental diagram from 7 : 00 to 9 : 00.
(a)(b)The MFD shows the network state for every 5 minutes between 7 : 00 and 9 : 00 during the full simulation. Before 7 : 00, the demand is low, and the network is almost empty. As such, we filtered out the first 30 minutes in order to present the critical part of the evolution of the network traffic conditions. For instance, point number 15 of the MFD (Figure3(b)) shows the total travel distance and total travel time of all the travelers in the network between 8 : 10 and 8 : 15. First, the MFD curve increases from (0, 0). The traffic states remain undersaturated (point numbers 1 to 8) when demand is light, in this case from 6 : 30 to 7 : 40. Afterwards, travel production, which is equivalent to the total travel distance for a given period of time, stabilizes, while the accumulation (or total travel time) continues to increase (point numbers 9 to 13). This corresponds to the saturation level occurring from 7 : 40 to 8 : 05. The decrease in travel production and accumulation (point numbers 14 to 24) shows that the network starts to exit the saturation level at the end of the simulation period and slowly returns to the undersaturation level from 8 : 05 to 9 : 00.The trips are divided into nine classes that are presented in Table2. In this table, the user class’s fraction of the total demand is presented by “percentage of trips,” and, for each class, we use a specific VOT. The demand level and VOT values are extracted from the real data of the full Lyon networks [42]. The optimization goal is to find the user equilibrium based on the generalized cost functions.Table 2
User classes and their value of time for the Lyon 6th + Villeurbanne network.
Percentage of trips (%)Class name (job)Value of time (euro/hour)3.20Worker category 152.566.40Worker category 239.2435.1Worker category 318.7226.70Student7.920.70Worker-student10.8018.80Retired9.363.2At home6.845.9Unemployed and other15.36The optimization process is executed for a maximum of 10 outer loops for five assignment periods. At each outer-loop iteration, new shortest path(s) can be discovered; therefore, if at least a new shortest path is found per outer-loop iteration, we have 11 paths per OD. In each outer loop, a maximum of 40 inner loops have been executed to find the best assignment pattern. Note that, in heuristic and SA algorithms, there is one more simulation per outer loop for keeping the best assignment pattern if the process is converged by the maximum number of the inner loops. This means that if the convergence does not happen during all inner loops, we take the best solution and simulate the network to prepare the result for the outer-loop convergence test. The convergence tests of two loops are fixed for all algorithms. In the inner-loop convergence test, the optimization will be stopped when we have less than 1% improvementGapβ=0.01. Improvement is related to the average gap between the current and previous inner loops. For the outer loop, we will stop the process when we do not find a new shortest path for all the ODs and also have a small value of gap or have no improvement in terms of total gap Gapα=0. In the next section, we will compare the performances of heuristic and metaheuristic algorithms in the multiclass multimodal DTA problem based on the numerical results.


## 5. Numerical Results and Discussion

The full optimization process is done for all methods. The aim of the experiments performed on the Lyon 6th + Villeurbanne network is to examine the algorithms’ convergence pattern of the solutions for a large-scale network case. Thus, the UE is calculated for the considered network seven times with the heuristic and metaheuristic algorithms defined inSection 3. The AGap indicator is used to evaluate the quality of the solution.The results for the performance indicators of all the algorithms are presented in Table3. GA leads to ∼10× more simulations, and SA leads to only ∼1.3× simulations compared to the ∼550 simulations required in heuristic algorithms. However, the computation times are significantly lower because of the parallel simulation framework. Moreover, the solutions obtained by the metaheuristic algorithms are significantly closer to the optimal UE than other methods. The metaheuristic algorithms dominate all heuristic methods in the final solution AGap; indeed, they are significantly better than the gap-based prob. algorithm, which is the best heuristic algorithm in this experiment. The GA and SA algorithms manage to reduce the UE AGap of the gap-based prob. method by more than 50% (Table 3). The violation indicator also shows that metaheuristic algorithms work much better than the heuristics methods (reduction of 54%).Table 3
Solution quality and performance indicators [CT: computation time].
Indicator/methodNumber of simulationsIncomplete travels (%)Violation (%),VGAGap (euros)Improvement toAGap compared to GBPCT (hours)Improvement to CT compared to SSPMax. number of cores usedMSA ranking5665,74%13.771.38−24.32%132.26−9.74%1Gap-based prob. (GBP)5944.08%9.031.11—137.52−14.10%1Probabilistic (prob.)5984.66%11.641.23−10.81%143.19−18.81%1Step-size prob. (SSP)5149.11%14.772.61−135.14%120.52—1Smart step-size prob.5915.16%10.221.27−14.41%139.59−15.81%1Simulated annealing7783.01%4.820.3271.17%91.3024.25%3Genetic algorithm58663.79%6.090.5451.35%83.7830.49%18Regarding the percentage of incomplete trips, also the metaheuristic algorithms dominate the others. The incomplete trips indicator denotes the share of trips, which is not finished by the end of the simulation in the final path flow distribution of each algorithm. A lower number of incomplete trips means a lower total travel time spent in the system over the simulation period. The SA algorithm finds the closest solution to UE (minimumAGap) in this study. Moreover, the final solution of SA has the best value for other quality indicators that are highlighted in bold in the SA row in Table 3.Moreover, we measure the computation time (CT) of algorithms to evaluate the performance of them. The CTs of SA algorithm and GA are significantly better than heuristic methods. In particular, it is significantly better than the SSP method (bold values in GA row inside Table3). Note also that the DTA process over the full simulation period requires considerable computational resources because of the network size. Therefore, the computational improvement obtained by the metaheuristic methods is high, 47 hours on average (two days) for the SA and GA, compared to the heuristic algorithms. The GA reduces the CT by one day and a half compared to SSP as the fastest heuristic method.Regarding the computation resources, every simulation usually needs one central processing unit (core) of the processor. The heuristic algorithms, as mentioned before, are run in series, so they need only one core per iteration. The SA algorithm generates a maximum of three new path flow patterns in itsliquid method per iteration. Therefore, we need a maximum of three cores for the simulation process. For the GA, according to population size, crossover rate, and the mutation rate, this value is maximum 30 new solutions from the GA in this study. Because all the experiments are conducted on a 64-bit personal computer with 18 cores, the number of cores is limited to 18. If the number of new solutions in GA is bigger than 18, the algorithm executes in two successive phases: the first 18 simulations in the first phase and the remaining solutions in the second phase.Figure4 presents the convergence pattern for all seven algorithms in the last assignment period. The convergence pattern is drawn based on the average delay of all users. By the best heuristic methods, the final user average gap is around 30 sec, but this value for the metaheuristic method is around 10 sec. It means that the metaheuristic methods also obtain the better solutions in terms of delay. In Figure 4, the convergence pattern of all heuristic algorithms first decreases and then variates in a small range for the last iterations, where they seem to stick to the local optimum. Metaheuristic algorithms’ convergence pattern is always decreasing, which means that the algorithms are exploring and looking for the global optimum. Therefore, both metaheuristic algorithms produce not only better solutions in terms of optimality but also better (i.e., not increasing) convergence patterns than classical algorithms.Figure 4
Convergence patterns for the outer loops.In order to consider the complexity of the problem on the performance of algorithms, a monoclass demand scenario is also designed. The volume of the demand and the traffic network are the same as those in the multiclass scenario. For monoclass users, we consider the weighted average VOT of the multiclass scenario (15.84 euros per hour), extracted from Table2, and the same monetary cost for public transportation. The setting for the optimization framework, including the maximum inner- and outer-loops’ iterations, is fixed to the same values as the multiclass scenario. Table 4 presents the monoclass scenario results with the multiclass scenario. Same as the previous experiment, metaheuristic algorithms dominate heuristic algorithms. The results show that heuristic algorithms cannot guarantee the user equilibrium for the large-scale network.Table 4
Comparison between the performances of algorithms with homogeneous and heterogeneous users.
Demand profile; quality indicator/methodMonoclass usersMulticlass usersAGap (euros)Improvement toAGapcompared to prob.CT (hours)Improvement to CTcompared to smart SSPAGap (euros)Improvement toAGapcompared to GBPCT (hours)Improvement to CT compared to SSPMSA ranking0.47−40.03%105.62−9.28%1.38−24.32%132.26−9.74%Gap-based prob. (GBP)0.42−23.88%101.81−5.33%1.11—137.52−14.10%Probabilistic (prob.)0.34—98.12−1.53%1.23−10.81%143.19−18.81%Step-size prob. (SSP)0.55−64.78%108.71−12.48%2.61−135.14%120.52—Smart step-size prob.0.37−10.45%96.65—1.27−14.41%139.59−15.81%Simulated annealing0.2136.36%85.6211.41%0.3271.17%91.3024.25%Genetic algorithm0.2718.48%72.8424.64%0.5451.35%83.7830.49%The results forAGap and CT of the multiclass setting are worse than those of the monoclass setting because of the VOT of different classes that change the scale of AGap and also require more CT to explore the solution space of the multiclass setting. However, regarding AGap, it is even more beneficial to use metaheuristic methods over heuristic methods in a multiclass environment, as the improvements in AGap are significantly better than those in the monoclass problem (Table 4). The same conclusion can be drawn for CT. The GA algorithm is the fastest algorithm in both scenarios. The performance comparison of solution algorithms in two demand scenarios shows that metaheuristic algorithms are more efficient when the problem’s complexity is increased.


## 6. Conclusion

This study develops the multiclass model by designing the hybrid generalized cost function VOT and cost function, including mode and economic classes. The model is extended for DTA, and variational inequality formulation is presented for the dynamic user equilibrium problem. This problem is also challenging from the computation point of view in the large-scale application. We compared two different assignment frameworks (heuristic and metaheuristic) to execute the model on the large-scale network in order to investigate which algorithm is more efficient. Seven optimization methods are discussed to calculate the user equilibrium. The model and assignment framework are applied to the Lyon 6th + Villeurbanne network, and the two optimization approaches are compared on the large-scale network. The multiattribute multiclass model and optimization process have been proven to be successful when applied to our multimodal large-scale network. The results show that the metaheuristic methods (SA and GA) with parallel computation work significantly better than the heuristic methods. The average user’s gap in generalized cost for metaheuristic algorithms is about 4.8 times in average better than the heuristic methods. Furthermore, the metaheuristic algorithms save a minimum of 36 hours in computation time. GA is the fastest algorithm, and SA provides the best solution for the network user equilibrium problem. Moreover, we can conclude that all heuristic methods do not guarantee the user equilibrium solution in large-scale DTA problems.As metaheuristic algorithms are new, there is room to improve them and also apply other kinds of metaheuristic algorithms. In addition to this, the authors propose the following topics for future works.1. Consider departure time choice models. One possibility can be adding another layer of optimization for the departure time adjustment process. It can also be interesting to use the proposed parallel framework to optimize departure time and route choices simultaneously. 2. Extend the multiattribute, multiclass, trip-based model to activity-based models in order to consider travel with activity chains. 3. Compare and analyze the performance of the metaheuristic method in other large-scale networks with different typologies and the other types of optimization methods. 4. Combine and investigate the metaheuristic methods with other kinds of traffic simulators and shortest path calculation methods.



---

*Source: 8815653-2021-02-19.xml*