A Methodology for Increasing Convergence Speed of Tra ﬃ c Assignment Algorithms Based on the Use of a Generalised Averaging Function

: In this paper, we propose a generalisation of the Method of Successive Averages (MSA) for solving tra ﬃ c assignment problems. The generalisation consists in proposing a di ﬀ erent step sequence within the general MSA framework to reduce computing times. The proposed step sequence is based on the modiﬁcation of the classic 1 / k sequence for improving the convergence speed of the algorithm. The reduction in computing times is useful if the assignment problems are subroutines of algorithms for solving Network Design Problems—such algorithms require estimation of the equilibrium tra ﬃ c ﬂows at each iteration, hence, many thousands of times for real-scale cases. The proposed algorithm is tested with di ﬀ erent parameter values and compared with the classic MSA algorithm on a small and on two real-scale networks. A test inside a Network Design Problem is also reported. Numerical results show that the proposed algorithm outperforms the classic MSA with reductions in computing times, reaching up to 79%. Finally, the theoretical properties are studied for stating the convergence of the proposed algorithm.


Introduction
The equilibrium traffic assignment problem has been widely studied since the early 1950s [1], both with respect to theoretical properties and proposing solution methods and algorithms. The Method of Successive Averages (MSA) is one of the most commonly used algorithms for solving the problem; it was proposed by References [2,3], and its convergence was stated by Reference [4]. The MSA proposed by References [2,3] operates on link flows (MSA-FA). On the other hand, Reference [5] proposed the MSA-CA which operates on link costs; while, Reference [6] devised the MSA-ACO based on Ant Colony Optimisation. Solution algorithms in the case of elastic demand were proposed and analysed by References [5,7,8].
The MSA algorithm at each iteration, k, averages the traffic flows calculated by stochastic uncongested network loading, with weight 1/k, with the results of the previous iteration, with weight (k-1)/k. Alternative step sequences (or weighting methods) were proposed, among others, by References [9][10][11]. More recently, Reference [12] proposed three revised versions of the traditional MSA method, namely: (i) Restarting MSA (RMSA), where the iteration counter is re-initialised after a specific number of iterations; (ii) Double RMSA (R2MSA), where the iteration counter is re-initialised to a variable value; (iii) Weighted MSA (WMSA), based on the self-regulated averaging method proposed by Reference [13] for solving the stochastic user equilibrium problem. The above-mentioned MSA In order to solve the fixed-point problems (4) we propose an algorithm based on the MSA general framework. The first formulation of MSA [2] generates a succession of feasible link flow vectors, f k , starting from an initial solution (usually f 0 = 0). At each iteration k, the solution f k is generated by combining traffic flows obtained by a Stochastic Uncongested Network (SUN) assignment, f SUN k = ϕ SUN (c k ), with the previous solution, f k−1 At each iteration k, in order to generate the solution f k , the results of the SUN assignment, f SUN k , are averaged, with weight 1/k, with the results of the previous iteration, f k−1 , with weight (k − 1)/k: that is: Since the algorithm works directly on flows, it is also known as MSA-FA (Flow Averaging). Algorithm 1 reports the code of the traditional MSA-FA algorithm, where ε represents the stop threshold.
It is worth noting that the averages applied to the flows belonging to successive iterations have to be considered an algorithmic trick and should not be confused with the weights that simulate the learning process in the dynamic assignment models. Indeed, the MSA algorithms apply to stationary models, while the simulation of learning processes (comparison between expected and experienced performances) are typical of day-to-day dynamics models [24,25].
Since, as shown in Appendix B, most of the MSA-based algorithms proposed in the literature for solving the traffic assignment problem are based on a modification of the weight formulation (i.e., term 1/k in the traditional MSA algorithm), our proposal consists in providing a further variant of the MSA-FA algorithm substituting the term k in Equation (5) or (6) with the following function: Appl. Sci. 2020, 10, 5698 4 of 23 where η is a parameter between 0 and 1; for η = 1 we obtain the classic MSA algorithm, and for η = 0 we obtain that f k = f SUN k (in this case the algorithm does not consider at each iteration the results of the previous ones). This method tends to muffle the weight given to f k−1 (with η < 1) with respect to the corresponding weight given with the classic MSA algorithm. The weights of f k−1 , w −1 , and f SUN k , w SUN , at each iteration k become: Algorithm 1 Traditional MSA-FA algorithm 1: In Algorithm 2 the code of the generalised MSA-FA algorithm (i.e., GMSA-FA) is reported. Figure 1 reports the values at each iteration k of w SUN (equal to 1 − w −1 ) for values of ηbetween 1 (classic MSA) and 0 and for the first 100 iterations.
It can be noted that the weight of f SUN k decreases less rapidly for lower values of parameter η.
Convergence of the proposed algorithm is stated in Appendix A, where the algorithm is shown to converge for each value of η > 0.

Numerical Results
Numerical tests were conducted on a small test network and on two real-scale networks (an urban and a rural network). Each network, by means of the graph theory, was represented as a set of links and nodes. The links represent the phases of the trips (as for instance, the roads) and the nodes represent the transitions between two consecutive links. Among the nodes, it is possible to identify the centroid nodes that represent the origin and destination of each trip.

Small Network
The small network consisted of 16 oriented links (see Figure 2), whose link performances (i.e., the link travel times) were expressed by means of a BPR function [26], that is: where t l is the travel time associated to link l; L l is the length of link l; v o,l is the free-flow speed associated to link l; f l is the link flow associated to link l; Cap l is the capacity of link l; α and β are function parameters. The main features of the considered network are summarised in Table 1.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 21 It can be noted that the weight of fSUN k decreases less rapidly for lower values of parameter η.
Convergence of the proposed algorithm is stated in Appendix A, where the algorithm is shown to converge for each value of η > 0.

Numerical Results
Numerical tests were conducted on a small test network and on two real-scale networks (an urban and a rural network). Each network, by means of the graph theory, was represented as a set of links and nodes. The links represent the phases of the trips (as for instance, the roads) and the nodes represent the transitions between two consecutive links. Among the nodes, it is possible to identify the centroid nodes that represent the origin and destination of each trip.

Small Network
The small network consisted of 16 oriented links (see Figure 2), whose link performances (i.e., the link travel times) were expressed by means of a BPR function [26], that is: where tl is the travel time associated to link l; Ll is the length of link l; vo,l is the free-flow speed associated to link l; fl is the link flow associated to link l; Capl is the capacity of link l; α and β are function parameters. The main features of the considered network are summarised in Table 1.    Bold links have higher capacity and free-flow speed with respect to the others. We assume that nodes 1, 2, 3 and 4 are also the origins and destinations of the trips, according to the OD vector, d, reported in Table 2. In the tests, we adopted Multinomial Logit as the stochastic path choice model and assumed that the cost of each link was equal to the running time expressed in minutes with parameter θ equal to 0.5 (coefficient of variation about 0.2). We implemented the proposed algorithm in Visual Basic. We tested the proposed algorithm with the following values of parameter η: 1 (classic MSA), 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.01, 0. For each value, we tested the OD vector reported in Table 2 multiplied by the following coefficients: 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, and 2.0. The use of different OD vectors, as well as different multipliers, is to analyse the algorithm's performance for different congestion levels. Indeed, the vectors thus obtained produce average saturation degrees (flow/capacity ratio) on the network between 0.249 and 0.838, as reported in Table 3. The convergence threshold was assumed to be equal to 1% for all tests. In all cases, the proposed algorithm was able to reduce the iterations for the convergence with some values of η lower than 1 (MSA). It can be noted that, upon reducing the value of η, the number of iterations starts to decrease until a minimum value is reached. The number of iterations then begins to increase until the case of η = 0 when, according to also to the theoretical properties examined in Appendix A, the algorithm may not converge.
In Table 4, the number of iterations required for convergence for all cases is reported (the algorithm was stopped in any event when the iterations reached 999); minimum values are written in bold. In the same table, the percentage variations in terms of iteration number are reported. However, it is worth noting that all variants of the algorithm provide the same results. The difference among them is only the number of iterations necessary to achieve the convergence condition. In Figure 3, the results of In bold are reported the best results. Note that the reduction in the number of iterations, and hence, in computing time, provided by the proposed algorithm can be very significant for medium-high congested networks, i.e., up to −61% compared with the traditional version of the MSA algorithm.

Real-Scale Urban Network
The real-scale urban test network is the network of Benevento, a town in the south of Italy with about 61,000 inhabitants. The network, as represented in Figure 4, has 1577 oriented links, comprising about 216 km of roads, 678 nodes and 80 centroids (66 internal and 14 external). The peak-hour OD matrix was made available by studies carried out while drawing up the Urban Traffic Plan. The route choice model is Multinomial Logit with an implicit enumeration of paths according to the algorithm proposed in Reference [27].
In bold are reported the best results.
Note that the reduction in the number of iterations, and hence, in computing time, provided by the proposed algorithm can be very significant for medium-high congested networks, i.e., up to −61% compared with the traditional version of the MSA algorithm.

Real-Scale Urban Network
The real-scale urban test network is the network of Benevento, a town in the south of Italy with about 61,000 inhabitants. The network, as represented in Figure 4, has 1577 oriented links, comprising about 216 km of roads, 678 nodes and 80 centroids (66 internal and 14 external). The peak-hour OD matrix was made available by studies carried out while drawing up the Urban Traffic Plan. The route choice model is Multinomial Logit with an implicit enumeration of paths according to the algorithm proposed in Reference [27]. In order to test the proposed algorithms under several network congestion levels, over the available OD matrix, seven other matrices were generated, multiplying all cells of the OD matrix by the following factors: 0.6, 0.8, 1.2, 1.4, 1.6, 1.8, and 2.0. The corresponding average saturation degrees are summarised in Table 5. On this network, we tested the proposed algorithm with the same values of parameter η as those adopted for the small test network; we tested the algorithm with two different values of the coefficient of variation (standard deviation to mean ratio), C v , equal to 0.10 and 0.05. Moreover, in these cases, the proposed algorithm was able to reduce the iterations for convergence with some values of η lower than 1 (MSA). In this case, the best performances were obtained with values between 0.4 and 0.6 (see Table 6) for C v = 0.10 and between 0.3 and 0.6 (see Table 7) for C v = 0.05; the reductions in the number of iterations reached −76% (C v = 0.10) and −79% (C v = 0.05) with respect to the traditional version of the MSA algorithm. In Figures 5 and 6, the results of Tables 6 and 7 are illustrated. Moreover, in this case, the utility of the proposed algorithm is more significant for medium-high demand levels.
In bold are reported the best results.
In bold are reported the best results.
In bold are reported the best results.

Real-Scale Rural Network
The real-scale rural test network is the regional road network of Campania, an Italian Region with more than 5.8 million inhabitants. The network, as represented in Figure 7, has 764 road links, representing about 8700 km of roads, 262 road nodes and 91 centroids.

Real-Scale Rural Network
The real-scale rural test network is the regional road network of Campania, an Italian Region with more than 5.8 million inhabitants. The network, as represented in Figure 7, has 764 road links, representing about 8700 km of roads, 262 road nodes and 91 centroids.
We tested the proposed approach considering the OD matrix in the morning peak-hour, where the total trips are 165,328 veh/h; we also tested values of demand between 0.6 and 2.0 times the actual value, producing an average saturation degree between 0.4 and 1.25 (see Table 8). Moreover, in this case, we use the Multinomial Logit route choice model with the implicit enumeration of paths [27]. On this network, we tested the proposed algorithm with the same values of parameter η previously adopted, assuming a coefficient of variation equal to 0.3. Moreover, in these cases, the proposed algorithm was able to reduce the iterations for convergence with some values of η lower than 1 (MSA). The best performances were obtained with values between 0.3 and 0.6 (see Table 9); the reductions in the number of iterations reached −69% with respect to the traditional version of the MSA algorithm. In Figure 8, the results of Table 9 are illustrated.
In bold are reported the best results.

Application to a Network Design Problem
The proposed algorithm was tested inside a procedure for solving Network Design Problem (NDP) on the rural real-scale network described in Section 4.3. In particular, we considered a network capacity expansion problem [28] where there were limited resources for improving the performances of a rural road network. In this problem, the current configuration of the road network is known, and no other roads can be added to it; instead, each existing road can be improved, increasing its capacity and free-flow speed. The decision variables of this problem are binary, yi = 0/1, where yi = 1 indicates that road i is improved in the solution, while yi = 0 indicates that road i remains in the starting configuration. The details of NDP and of its model formulation are reported in Appendix C.
In order to compare the results obtained with the proposed muffled algorithm and the classic MSA, we tested the Neighbourhood Search Algorithm (NSA) with steepest descent method; it was used, for instance, in References [28,29]. A brief description of the algorithm is reported in Appendix D.
The tests were conducted under the actual demand (1.0 d), and therefore, assuming the parameter  = 0.5 in the muffled proposed algorithm, since it is the best value for this demand level (see Table 9). The results of the comparison between MSA and muffled algorithm are summarised in

Application to a Network Design Problem
The proposed algorithm was tested inside a procedure for solving Network Design Problem (NDP) on the rural real-scale network described in Section 4.3. In particular, we considered a network capacity expansion problem [28] where there were limited resources for improving the performances of a rural road network. In this problem, the current configuration of the road network is known, and no other roads can be added to it; instead, each existing road can be improved, increasing its capacity and free-flow speed. The decision variables of this problem are binary, y i = 0/1, where y i = 1 indicates that road i is improved in the solution, while y i = 0 indicates that road i remains in the starting configuration. The details of NDP and of its model formulation are reported in Appendix C.
In order to compare the results obtained with the proposed muffled algorithm and the classic MSA, we tested the Neighbourhood Search Algorithm (NSA) with steepest descent method; it was used, for instance, in References [28,29]. A brief description of the algorithm is reported in Appendix D.
The tests were conducted under the actual demand (1.0 d), and therefore, assuming the parameter η = 0.5 in the muffled proposed algorithm, since it is the best value for this demand level (see Table 9). The results of the comparison between MSA and muffled algorithm are summarised in Table 10; note that the computing time saved is about 44% (over 2.5 h). The different numbers of examined neighbourhoods and objective function values are only due to the convergence test; indeed, even if the convergence test is the same, the flows generated at each iteration of the algorithms can be slightly different and may influence the final solution. In Figure 9, the comparison among MSA and muffled algorithm is reported in terms of computing times and objective function values. These results underline the utility of the proposed algorithm inside Network Design Problem solution procedures; indeed, we have to consider also that several meta-heuristic algorithms use NSAs as subroutines (Scatter Search, Tabu Search, etc.), and then, the computing time saved should be multiplied by the number of NSAs performed for reaching the final solution.  Figure 9. Comparison between algorithms in terms of computing times vs. objective function values.

Conclusions and Research Prospects
The results obtained on the small and on the real-scale networks show that the proposed algorithm is able to reduce the number of iterations for convergence, and hence, computing times compared with the classic MSA approach; in particular, the reductions arrive in percentage up to -79% with respect to the classic MSA for the best value of the parameter. On all networks the results are similar, showing that the benefits of the proposed algorithm are significant for medium-high congested demand levels (average saturation degrees over about 0.4); less significant, but not negligible, are the benefits for lower demand levels.
In all tests, the best values of parameter η lie between 0.3 and 0.6. An important characteristic of the proposed algorithm is that there is a parameter that can be opportunely chosen, so to optimise the performances before starting with network design procedure; moreover, examining Figures 3, 5, 6 and 8 it is evident that, especially for medium-high levels of demand, a minimum can be found. In our test, we adopt a 0.1 step, but the parameter η can assume any value between 1 (MSA) and 0 (0 excluded), and therefore, other values can be tested in order to find the one that minimises the computing times.
The reduction in computing times can be not important if the assignment procedure is performed only once or a few times: with current computing capabilities, the savings are lower than a minute on the real-scale network. By contrast, if the assignment is a subroutine of the Network Design Problem, that requires the calculation of the equilibrium traffic flows many thousands of times, the reduction in computing time of an assignment is transferred to the computing time of the network design algorithm, allowing as much as several days to be saved in the calculation. In our test on the rural road network, the proposed algorithm is able to save about 44% of computing time, equivalent to about 2.5 h for a single NSA procedure.
Future research will focus on testing other real-scale networks, testing the algorithm on multimodal networks and under the assumption of elastic demand and proposing a similar algorithm for solving the combined assignment-control problem. Moreover, testing the proposed algorithm within other real-scale network design procedures will be the subject of future studies.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Theoretical Properties
There are three theoretical properties of a mathematical problem and its solution algorithm: existence, uniqueness and convergence. The first two properties are related to the mathematical problem and consist in proving that there is at least one solution (existence), and there is at most one solution (uniqueness). The combination of these two properties implies that the solution to the problem exists and is unique.
Once the existence and uniqueness of the problem are stated, it is necessary to adopt or develop suitable solution algorithms that are able to provide the numerical solution to the problem. Hence, it is necessary to state theoretically the convergence (i.e., the ability to provide the solution in a finite number of iterations) of the proposed algorithms.
The conditions ensuring the existence and the uniqueness of traffic assignment problem can be found in Reference [5]. As for convergence, the same paper proposed an extension of Blum's theorem [30] where it was shown that a convergent solution algorithm for solving the fixed-point problem of a function y = λ(x) which satisfies some suitable conditions could be formulated by the following recursive equation: If the sequence µ k k>0 satisfies the following conditions: In particular, the above-mentioned suitable conditions are: 1.
x ∈ S x and y ∈ S x where S x is a non-empty, compact and convex set; 2. a unique solution to the fixed-point problem x * = λ x * ; 3. a function ϕ(x) ≥ 0 ∀x ∈ S x where ϕ(x) is continuous with first ∇ϕ(x) and second ∇ 2 ϕ(x) derivative continuous; x * ∈ Sx;
In Reference [5], it is showed that if the sequence µ k k>0 satisfies the condition: then the elements of the sequence described by Equation (A1) belong to set S x , which is convex. Moreover, by fixing ∇ϕ(x) = c(x) − c * , with x = f and c * = c f * , all Conditions (1)-(6) are satisfied. Hence, the convergence proof is only related to checking Conditions (A2)-(A4). Finally, Reference [5] showed that the traditional approach of MSA, i.e., ξ(k) = k, satisfies (A2)-(A4) conditions, since: In the case of the proposed generalised function, i.e., ξ(k) = 1 + ((k − 1) · η), condition (A4) is satisfied, and since for k → +∞ we obtain ξ(k) ∼ η · k, it is possible to show that: i.e., all convergence conditions are satisfied for the proposed algorithm for each value of η > 0.

Appendix B. Algorithm Performance Comparison
Most of the MSA-based algorithms proposed in the literature for solving the traffic assignment problem are based on a modification of the weight formulation (i.e., term 1/k in the traditional MSA algorithm) by adopting a 1/ξ(k) function, that is: The differences in formulation among these algorithms are related to the adopted ξ(k) function. Indeed, the traditional MSA algorithm [2,3] assumes: Likewise, our proposal described in Section 3 adopts: Obviously, in the case of η = 1, our proposal provides the traditional MSA formulation, that is (A11) degenerates into (A10).
The comparison among algorithm formulations has been implemented in the case of the urban network of Benevento (described in Section 4.2) by adopting a Cv value equal to 0.05. Numerical results are shown in Table A1 and Figure A1.
Numerical results have shown that the proposed approach (i.e., the muffled method), at least in the case under consideration (the urban network of Benevento with a Cv coefficient equal to 0.05) almost always provides the best results, i.e., it requires the minimum number of iterations to reach convergence.
In bold are reported the best results.
Numerical results have shown that the proposed approach (i.e., the muffled method), at least in the case under consideration (the urban network of Benevento with a Cv coefficient equal to 0.05) almost always provides the best results, i.e., it requires the minimum number of iterations to reach convergence.
It is worth noting that, although the convergence can be demonstrated theoretically (as already shown in Appendix A), a theoretical proof cannot be performed in the case of the convergence speed property which represents an empirical element to be verified by means of numerous and appropriate tests. However, although it is not possible to generalise the optimal performance of the proposed methodology in an absolute way, the results obtained give hope that at least the proposed method can be included among the ones providing optimal performances. Obviously, the confirmation of this statement can only be obtained through a large number of tests in the case of different networks both It is worth noting that, although the convergence can be demonstrated theoretically (as already shown in Appendix A), a theoretical proof cannot be performed in the case of the convergence speed property which represents an empirical element to be verified by means of numerous and appropriate tests.
However, although it is not possible to generalise the optimal performance of the proposed methodology in an absolute way, the results obtained give hope that at least the proposed method can be included among the ones providing optimal performances. Obviously, the confirmation of this statement can only be obtained through a large number of tests in the case of different networks both in terms of context (urban vs. rural), network dimension (urban, provincial or regional network), demand and/or congestion levels.   Figure A2. Identification of decision variables.

Appendix D. Network Design Solution Algorithms
In order to solve the NDP described in Appendix C, we have used a steepest descent Neighbourhood Search Algorithm (NSA); this algorithm was previously used for solving NDPs as a subroutine of meta-heuristic algorithms, among others, by References [28,29].
Given a solution, y, the set of solutions that can be obtained from solution y changing only one value of a variable y i , from 0 to 1 or vice versa, is called neighbourhood of y, N(y); note that other kinds of neighbourhoods can be defined, depending on the kind of problem and on the dimension of neighbourhood (for instance, changing more than one variable per time). The NSA, starting from a solution, y k , generates the following solution, y k+1 , so that: and, with the steepest descent approach, z(y k+1 ) = Min {z(y); ∀ y ∈ N(y k )} (A22) where z(y) formally represents the objective function. The procedure then generates, at each iteration, a solution better than the previous one, choosing, among all solutions belonging to the neighbourhood, the one with the best objective function value. The procedure ends when solution y k is a local optimum, that is when: z(y k ) ≤ z(y) ∀ y∈ N(y k ) (A23) This algorithm is able to find a local optimum, and therefore, is often used inside meta-heuristic algorithms that need local search subroutines. In the test reported in this paper, we use the NSA only one time, starting from the current solution: y 0 = 0.