Efficient Algorithm for the Nonadditive Traffic Assignment Problem With Link Capacity Constraints

Abstract

This paper presents an insightful examination of the modeling and efficient solution algorithm for the link capacitated nonadditive traffic assignment problem (CNaTAP) to provide highly accurate flow solutions for large-scale networks. Despite the increasing significance of the CNaTAP, the ability to efficiently solve it for satisfactory accuracy in practical applications remains inadequate. Given that existing CNaTAP models and algorithms are typically limited to small experimental networks, the CNaTAP model is formulated as a variational inequality (VI) problem in this paper. This formulation is decomposed into two VI subproblems that involve equilibrium and capacity constraints, utilizing the Karush–Kuhn–Tucker (KKT) conditions. The Lagrangian multipliers for the capacity constraints are treated as fixed costs for the links in the equilibrium subproblem, ensuring the stability of the Cartesian product structure within the feasible set. This approach facilitates the decomposition of OD pairs, enabling the efficient solution of CNaTAP in large-scale networks. In addition, an algorithmic framework is developed that incorporates high-frequency updates of these Lagrangian multipliers, along with an adaptive Barzilai–Borwein (ABB) step-size calculation method applied to expedite convergence in the equilibrium subproblem. Extensive numerical experiments confirm the efficacy of the proposed algorithm in efficiently solving large-scale networks with high convergence accuracy.

This work is licensed under http://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

1. Introduction

The traffic assignment problem (TAP) entails allocating predicted Origin-Destination (OD) demand across the transportation network based on specific rules to derive flow information for each link. As a crucial aspect of transportation planning, the TAP finds extensive application in predicting network traffic distribution patterns, assessing the potential impacts of projects and policies, and analyzing the utilization of existing roadway networks to manage congestion, emissions, and toll revenues. The traditional TAP models assume that travelers adhere to the user equilibrium (UE) principle, i.e., select the path that minimizes travel costs. The path travel cost is typically assumed to be the simple summation of link travel costs along the path. This assumption of path cost additivity enables the formulation of TAP as a convex program without incorporating path flow variables in the objective function. Consequently, while path flow variables remain in the constraint set, explicit storage of path information is unnecessary when solving TAP, providing significant computational advantages for large-scale networks. Numerous efficient algorithms have been developed to solve large-scale TAP [1, 2].

However, in many practical scenarios, the assumption of additivity does not hold [3, 4]. Certain path attributes, such as path-specific tolls [5] or travel time reliability [6], cannot be simply represented as the sum of corresponding link attributes. Moreover, travel decisions often involve a comprehensive assessment of multiple cost factors, including time, money, distance, safety, complexity, and reliability. A common modeling approach integrates these attributes into the path cost measure in a nonlinear manner [7, 8], resulting in nonadditive cost characteristics. The nonadditive traffic assignment problem (NaTAP) is typically formulated as a variational inequality (VI) problem, a nonlinear complementarity problem, or a nonconvex mathematical programming model based on a gap function. Among these approaches, VI-based modeling has garnered significant attention, leading to the development of various solution algorithms applicable to real-world transportation networks [9–11].

While NaTAP improves the realism of the traditional TAP, it still fails to fully capture various control policies and traffic constraints in urban transportation networks. For instance, capacity constraints are a fundamental limitation in transportation systems and cannot be overlooked. A direct approach to enhancing the realism of NaTAP is to incorporate link capacity constraints, which prevent link flows from reaching impractical levels and provide a more accurate representation of queuing and congestion effects.

The feasible sets of the TAP and the NaTAP exhibit a Cartesian product structure, facilitating efficient resolution through OD pair decomposition. In fact, nearly all algorithms for large-scale TAP and NaTAP exploit this property [11, 12]. In contrast, TAPs with side constraints present greater computational challenges, as these constraints disrupt the Cartesian product structure of the feasible set. This disruption prevents the problem from being decomposed into independent subproblems solvable for each OD pair separately, necessitating a joint solution approach. Moreover, the shortest path problem in such cases transforms into a multicommodity minimum-cost flow problem, significantly increasing computational complexity. Table 1 summarizes existing research on TAPs with side constraints.

Table 1

Summary of existing studies on TAPs with various side constraints.

Path travel cost components	Reference	Application scenarios	Methodology	Computed large-scale networks
Additive	Daganzo [13, 14]	Capacity constraints	Asymptotic link performance function	No
	Yang, and Yagar [15, 16]	Capacity constraints	IPF	No
	Prashker, and Toledo [17]	Capacity constraints	IPF	No
	Shahpar et al. [18]	Capacity constraints	IPF	No
	Larsson and Patriksson [19, 20]	Capacity constraints	ALM	No
	Lawphongpanich [21]	Capacity constraints	ALM	No
	Larsson et al. [22]	Generalized constraints	ALM	No
	Nie et al. [23]	Capacity constraints	ALM and IPF	Yes
	Feng et al. [24]	Capacity constraints	ALM	Yes
Nonadditive	Chen et al. [25]	Capacity and environmental constraints	VI (decomposition algorithm)	No
	Xu et al. [26]	Environmental constraints	NCP (steepest descent method)	No

Currently, the mainstream approaches for solving link capacitated traffic assignment problems (CTAPs) fall into two categories: (1) the asymptotic link performance function approach. (2) The unconstrained reformulation approach via a penalty/dual scheme. Daganzo [13, 14] was the first to implicitly account for link capacity constraints through asymptotic link travel-time functions. While this approach yielded feasible solutions, modifying the travel-time function was impractical. Boyce et al. [27] highlighted that adopting an asymptotic link performance function could result in excessively high travel times, convoluted rerouting of trips, and numerical difficulties in the neighborhood of the capacity. In contrast, the second method offers greater flexibility, including the interior penalty function (IPF) and the augmented Lagrangian multiplier (ALM) methods. The IPF method establishes barriers on the boundary of the feasible set to ensure constraint adherence during the solution process. The ALM method, which integrates dual methods with external penalty techniques, relaxes constraints to formulate an approximate convex programming. Nie et al. [23] compared the IPF and ALM methods, employing the gradient projection (GP) algorithm as the subproblem solver. Their findings indicate that the ALM method offers superior computational efficiency and stability. Furthermore, Feng et al. [24] integrated the greedy algorithm into the ALM method framework for solving the CTAP, with numerical experiments demonstrating its effectiveness and high convergence accuracy on large-scale networks. However, these studies continue to assume additivity in path costs, restricting their applicability to the link capacitated nonadditive traffic assignment problem (CNaTAP).

Previous studies have investigated various modeling and solution approaches for general TAP, including NaTAP, with multiple side constraints. For instance, Chen et al. [25] formulated a generalized side-constrained traffic equilibrium model within a VI framework, incorporating both capacity and environmental constraints. They proposed a decomposition algorithm to solve the VI problem under relatively mild assumptions. Similarly, Xu et al. [26] employed a gap function approach to reformulate the general TAP with environmental constraints as an equivalent unconstrained optimization problem, which was subsequently solved using the steepest descent method. However, it is important to note that the methods proposed by Chen et al. [25] and Xu et al. [26] were tested only on small-scale networks, which are considered simple toy networks by current standards.

While several studies have introduced solution methods for CTAPs, most suffer from significant limitations: they either fail to effectively address the nonadditivity of path costs or are restricted to small-scale networks, making them difficult to extend to real-world transportation systems. Therefore, this paper aims to explore an efficient and accurate algorithm for solving CNaTAP in large-scale networks. The key contributions of this study are as follows:

• CNaTAP modeling: the CNaTAP is formulated within a VI framework and decomposed into an equilibrium subproblem and a capacity constraints subproblem based on the Karush–Kuhn–Tucker (KKT) conditions.

• Solution method optimization: in solving the equilibrium subproblem, the Lagrange multipliers associated with capacity constraints are treated as fixed link costs. This approach restores the Cartesian product structure of the feasible set, enabling OD pair decomposition and enhancing computational efficiency for large-scale networks.

• Efficient algorithmic framework: an optimized algorithmic framework is designed to increase the update frequency of Lagrange multipliers, thereby accelerating algorithm convergence.

The remainder of the paper is organized as follows: Section 2 introduces the VI formulation of the CNaTAP and its fundamental properties. Subsequently, Section 3 presents an efficient CNaTAP solution method suitable for large-scale networks. Section 4 outlines the specific implementation details for solving the CNaTAP. Section 5 verifies the properties of the CNaTAP solution and demonstrates the efficiency and high accuracy of our algorithm through a series of numerical experiments. Finally, Section 6 provides a comprehensive paper summary and proposes potential future research directions.

2. Formulation of the CNaTAP

This paper extends the conventional UE TAP model to a more realistic CNaTAP model, incorporating complex cost structures and link capacity side constraints. There are three critical general assumptions: (1) traffic demands are predetermined and fixed, and the traffic conditions in the network are deterministic; (2) travelers possess comprehensive knowledge of the network and select paths with minimize travel costs; and (3) the entire network maintains adequate capacity to fulfill traffic demand. Assumption (1) may be relaxed by introducing elastic demand, while assumption (2) can be expanded by integrating stochastic traffic equilibrium models. Assumption (3) ensures that a solution for the CNaTAP model exists.

This study considers a fixed OD travel demand in a directed traffic network $GN,A$, where $N$ and $A$ denote the sets of nodes and links, respectively. Travel demand $q^{rs}$ is generated from origin $r\in R$ to destinations $s\in S$, where $R$ and $S$ are a subset of $N$. A set of paths $P^{rs}$ connects origin $r$ to destination $s$, with flow $f_p^{rs}$ on each path $p\in P^{rs}$. Let $v_a$ represent the flow on the link $a\in A,$ and $c_a$ represent the capacity of the link $a\in A.$ Let $t_a$ represent the time of the link $a\in A.$ Let ${\delta }_{pa}^{rs}=1$ if the path $p$ from origin $r$ to destination $s$ through link $a$, and 0 otherwise.

Let ${\eta }_p^{rs}$ denote the cost of the path $p$ between OD $rs.$ A general nonadditive path cost function can be formulated as follows [11]:$$\begin{array}{}\text{(1)}& {\eta }_p^{rs}={\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}{\rho }_a+{\displaystyle \sum }_{a\in A}\sigma {\delta }_{pa}^{rs}{t}_a+h_p{\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}{t}_a,\forall a\in A,p\in P^{rs},r\in R,s\in S,\end{array}$$where ${\rho }_a$ represents the direct financial cost (e.g., tolls) on the link $a;$$\sigma$ denotes the operating cost per unit of travel time (e.g., fuel consumption and vehicle rental), and $h_p$ is a function of time value (converting travel time into monetary cost). $h_p$ and $t_a$ continuously and monotonically increase for path travel time and link flow, respectively.

If each traveler tries to minimize their travel cost, the CNaTAP can be formulated as the following VI problem:$$\begin{array}{}\text{(2)}& {\mathbf{\eta }{\mathbf{f}}^{\mathrm{\ast }}}^{\mathrm{\Gamma }}\mathbf{f}-{\mathbf{f}}^{\mathrm{\ast }}\ge 0,\forall \mathbf{f}\in \mathrm{\Phi },\end{array}$$where $\mathbf{f}$ denotes the vector of $\dots ,f_p^{rs},\dots ,$$\mathbf{\eta }$ denotes the path cost vector $\dots ,{\eta }_p^{rs},\dots ,$$\mathbf{\eta }{\mathbf{f}}^{\mathrm{\ast }}$ is the path cost under the path flow pattern ${\mathbf{f}}^{\mathrm{\ast }},$ and $\mathrm{\Phi }$ is described by equations (3)–(6).$$\begin{array}{}\text{(3)}& {\displaystyle \sum }_{p\in P^{rs}}{f}_p^{rs}-q^{rs}=0,\forall r\in R,s\in S,\text{(4)}& v_a={\displaystyle \sum }_{r\in R}{\displaystyle \sum }_{s\in S}{\displaystyle \sum }_{p\in P^{rs}}{f}_p^{rs}{\delta }_{pa}^{rs},\forall a\in A,\text{(5)}& f_p^{rs}\ge 0,\forall p\in P^{rs},r\in R,s\in S,\text{(6)}& g_a=v_a-c_a\le 0,\forall a\in A.\end{array}$$

Equation (3) represents a conservation constraint, while equation (4) is a defining constraint that summarizes all path flows through a given link. Equation (5) imposes a non-negativity constraint on path flow, and equation (6) establishes a capacity constraint on the link.

The abovementioned VI formulation for the CNaTAP demonstrates extensive applicability. It can manage diverse types of side constraint structures and path cost configurations. Consequently, this formulation is not limited to addressing solely nonadditive path costs; it can also accommodate additive path costs. Moreover, the capacity side constraints delineated in the formulation can be substituted with alternative side constraints, such as various complex environmental constraints. Furthermore, side constraints can be extended to different spatial levels, including the link, node, and region. Let ${\pi }^{rs}$ and ${\lambda }_a$ represent the Lagrange multiplier associated with the flow conservation constraints (3) and the link capacity side constraints (6), respectively. We derive the KKT conditions of the VI (2) as follows:$$\begin{array}{}\text{(7)}& \begin{array}{c}{f}_p^{rs}{\eta }_p^{rs}+{\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}{\lambda }_a\frac{\mathrm{\partial }{g}_a}{\mathrm{\partial }{v}_a}-{\pi }^{rs}=0\\ {\eta }_p^{rs}+{\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}{\lambda }_a\frac{\mathrm{\partial }{g}_a}{\mathrm{\partial }{v}_a}-{\pi }^{rs}\ge 0\\ f_p^{rs}\ge 0\end{array},\forall p\in P^{rs},r\in R,s\in S,\text{(8)}& \begin{array}{c}{\pi }^{rs}{\displaystyle \sum }_{p\in P^{rs}}{f}_p^{rs}-q^{rs}=0\\ {\displaystyle \sum }_{p\in P^{rs}}{f}_p^{rs}-q^{rs}\ge 0\\ {\pi }^{rs}\ge 0\end{array},\forall r\in R,s\in S,\text{(9)}& \begin{array}{c}{\lambda }_a-g_a=0\\ -g_a\ge 0\\ {\lambda }_a\ge 0\end{array},\forall a\in A.\end{array}$$

The KKT condition outlined above comprises three subsystems of nonlinear complementary conditions. Equation (7) explicitly states the UE condition for the generalized travel costs of the path, asserting that the generalized travel cost of all utilized paths does not exceed that of unused paths, where ${\pi }^{rs}$ can be understood as the minimum generalized path travel cost of OD $rs.$ Let ${\tilde{\eta }}_p^{rs}={\eta }_p^{rs}+{{\displaystyle \sum }}_{a\in A}{\delta }_{pa}^{rs}{\lambda }_a,$ where ${\tilde{\eta }}_p^{rs}$ denotes the generalized travel cost of path $p$ and $\tilde{\mathbf{\eta }}$ represents its vector form. It follows that equation (7) is equivalent to the following VI:$$\begin{array}{}\text{(10)}& \tilde{\mathbf{\eta }}{{\mathbf{f}}^{\mathrm{\ast }}}^{\mathrm{T}}\mathbf{f}-{\mathbf{f}}^{\mathrm{\ast }}\ge 0,\forall \mathbf{f}\in \mathrm{\Omega },\end{array}$$where $\mathrm{\Omega }$ is described by equations (3)–(5).

Equation (8) represents the equilibrium condition associated with the conservation constraint. Equation (9) encapsulates the equilibrium condition associated with the link capacity constraint. In order to grasp the physical implications of ${\lambda }_a,$ equation (9) can be summarized as follows:$$\begin{array}{}\text{(11)}& {\lambda }_a>0⟹g_a=0,\forall a\in A,\text{(12)}& {\lambda }_a=0⟹g_a\le 0,\forall a\in A.\end{array}$$

A specific mapping correlation exists between ${\lambda }_a$ and $g_a,$ represented as $g_a=g{\lambda }_a.$ Given that ${\lambda }_a$ is positive solely when the link $a$ is saturated, it is rational to attribute it to a negative utility associated with constrained capacity. In this context, ${{\displaystyle \sum }}_{a\in A}{\delta }_{pa}^{rs}{\lambda }_a$ can be understood as the additional travel cost incurred by travelers attempting to use a saturated path. Denoting the vector $\dots ,{\lambda }_a,\dots$ as $\mathbf{\lambda }$ and the vector $\dots ,g_a,\dots$ as $\mathbf{g}$, equation (9) is equivalent to the following VI:$$\begin{array}{}\text{(13)}& -\mathbf{g}{{\mathbf{\lambda }}^{\mathrm{\ast }}}^{\mathrm{T}}\mathbf{\lambda }-{\mathbf{\lambda }}^{\mathrm{\ast }}\ge 0,\forall \mathbf{\lambda }\in \mathrm{\Lambda },\end{array}$$where $\mathbf{g}{\mathbf{\lambda }}^{\mathrm{\ast }}$ is the link capacity side constraint under ${\mathbf{\lambda }}^{\mathrm{\ast }}$, and $\mathrm{\Lambda }$ is described by ${\lambda }_a\ge 0,\forall a\in A$.

In summary, VI (2) can be decomposed into two VIs, represented by equations (10) and (13), under the condition of given and fixed traffic demand. An efficient solution to the CNaTAP can be achieved by jointly addressing VI (10) and VI (13).

3. The Solution Algorithm for the CNaTAP

The efficient OD decomposition algorithm for solving NaTAP cannot be directly applied to CNaTAP, as the introduction of side constraints disrupts the Cartesian product structure of the feasible set. In this section, we propose an approach to enable OD decomposition. Specifically, in solving the equilibrium subproblem (10), treating the Lagrangian multipliers $\mathbf{\lambda }$ as a fixed cost for links temporarily circumvents capacity constraints, ensuring the stability of the Cartesian product structure and facilitating OD pair decomposition. Sequential solving of the decomposed equilibrium subproblems allows focused attention on individual OD:$$\begin{array}{}\text{(14)}& \tilde{\mathbf{\eta }}{{\mathbf{f}}^{\mathrm{\ast }}}^{\mathrm{T}}\mathbf{f}-{\mathbf{f}}^{\mathrm{\ast }}\ge 0,\forall \mathbf{f}\in {\mathrm{\Omega }}^{rs},\end{array}$$where ${\mathrm{\Omega }}^{rs}=f_p^{rs}\in {\mathbb{R}}_{+}^n:{{\displaystyle \sum }}_{p\in P^{rs}}{f}_p^{rs}=q^{rs},n=P^{rs}.$

Due to the nonnegative property of the flow feasible set ${\mathrm{\Omega }}^{rs}$, several methods are available for computing the projection on ${\mathrm{\Omega }}^{rs}$. Similarly, $\mathrm{\Lambda }$ represents a simple nonnegative space, where projections can be computed with minimal effort. Based on these properties, this section employs a projection method to alternately solve VI (13) and VI (14), ultimately obtaining the solution to VI (2). Specifically, in solving the equilibrium subproblem (14), the Lagrangian multipliers are held constant, while in the capacity constraints subproblem (13), the path flow is fixed. These two steps alternate iteratively to solve the original problem (2).

From equations (11) and (12), it follows that the mapping $g_a=g{\lambda }_a$ is monotonic. When the Lagrangian multipliers are treated as the fixed costs of the links, the generalized travel cost increases strictly monotonically with respect to path flow. Furthermore, both feasible sets ${\mathrm{\Omega }}^{rs}$ and $\mathrm{\Lambda }$ are closed convex sets, which guarantees the convergence of the projection algorithm for solving the equilibrium and capacity constraints subproblems.

Since both $g_a$ and ${\tilde{\eta }}_p^{rs}$ maintain monotonicity, and the alternating process preserves this property, the alternating projection method guarantees that each iteration progressively converges to the solution of the original problem. Moreover, the projection operator is nonexpansive on closed convex sets, ensuring that the error does not increase with each iteration. As a result, the iterative sequence generated by the alternating projection method remains bounded, with the error gradually decreasing until convergence is achieved. However, it is important to note that ${\lambda }_a$ may not always be uniquely determined at the optimum solution [19].

3.1. Method for Solving the Equilibrium Subproblem

Goldstein [28] and Levitin et al. [29] initially proposed the projection method, commonly called the GLP projection algorithm. The iterative scheme of equation (14) based on the GLP projection is summarized as follows:$$\begin{array}{}\text{(15)}& {\mathbf{f}}^{k+1}={\text{Proj}}_{{\mathrm{\Omega }}^{rs}}{\mathbf{f}}^k-{\beta }^k\tilde{\mathbf{\eta }}{\mathbf{f}}^k,k=1,2\dots ,\end{array}$$where ${\text{Proj}}_{{\mathrm{\Omega }}^{rs}}\mathbf{x}$ denotes the orthogonal projection from $\mathbf{x}$ to ${\mathrm{\Omega }}^{rs}$, and ${\beta }^k$ represents the step-size of the $k\text{th}$ iteration. A larger step-size may result in divergence, while a smaller step-size may decelerate the convergence rate. Recently, Tan et al. [11] have made significant advancements in determining suitable step-size. They have applied the Barzilai–Borwein (BB) step-size, originally proposed by Barzilai and Borwein [30], to enhance the performance of GP algorithms for solving NaTAP. The computation of the BB step-size involves intelligently utilizing information from previous iterations, thereby avoiding frequent projection operations and mapping evaluations, which makes it particularly effective for solving large-scale problems. The following is the formula for computing the BB step-size in equation (15):$$\begin{array}{}\text{(16)}& {\beta }_{\text{BB}1}^k=\frac{{{\mathbf{f}}^k-{\mathbf{f}}^{k-1}}^2}{{{\mathbf{f}}^k-{\mathbf{f}}^{k-1}}^{\mathrm{T}}\tilde{\mathbf{\eta }}{\mathbf{f}}^k-\tilde{\mathbf{\eta }}{\mathbf{f}}^k},\end{array}$$or$$\begin{array}{}\text{(17)}& {\beta }_{\text{BB}2}^k=\frac{{{\mathbf{f}}^k-{\mathbf{f}}^{k-1}}^{\mathrm{T}}\tilde{\mathbf{\eta }}{\mathbf{f}}^k-\tilde{\mathbf{\eta }}{\mathbf{f}}^{k-1}}{{\tilde{\mathbf{\eta }}{\mathbf{f}}^k-\tilde{\mathbf{\eta }}{\mathbf{f}}^{k-1}}^2}.\end{array}$$

The ${\beta }_{\text{BB}1}^k$ and the ${\beta }_{\text{BB}2}^k$ are denoted as the long BB step-size and short BB step-size, respectively. Drawing inspiration from the superlinear behavior of the BB method in two-dimensional quadratics, Zhou et al. [31] proposed an adaptive Barzilai–Borwein (ABB) step-size calculation method capable of automatically selecting a suitable small or large step-size in each iteration. In this approach, the small step-size can provide a favorable descent direction for the next iteration, while the large step-size can yield a significant reduction. When ${\beta }_{\text{BB}2}^k/{\beta }_{\text{BB}1}^k$ is small, they believe that the gradient norms have not taken a satisfactory descent in the previous iteration and, therefore, select a smaller step-size, ${\beta }_{\text{BB}2}^k$. Conversely, they believe that the gradient norms have achieved or can achieve a satisfactory descent in the prior iteration and thus choose a larger step-size, ${\beta }_{\text{BB}1}^k$. The specific expression is presented in the following:$$\begin{array}{}\text{(18)}& {\beta }^k=\begin{array}{cc}{\beta }_{\text{BB}1}^k,& \text{if}\frac{{\beta }_{\text{BB}2}^k}{{\beta }_{\text{BB}1}^k}<\kappa ,\\ {\beta }_{\text{BB}2}^k,& \text{otherwise},\end{array}\end{array}$$where $0<\kappa <1.$ If a degenerate case ${{\mathbf{f}}^k-{\mathbf{f}}^{k-1}}^{\mathrm{T}}\tilde{\mathbf{\eta }}{\mathbf{f}}^k-\tilde{\mathbf{\eta }}{\mathbf{f}}^{k-1}\le 0$ occurs, the step-size can be set to a predefined parameter, for instance, ${\beta }^0/10$, where ${\beta }^0$ denotes the initial step-size and satisfies $0<{\beta }^0<1$. Numerical experiments conducted in this paper demonstrate that utilizing the ABB method yields superior performance compared to the regular BB method (refer to Section 5.3).

3.2. Method for Solving the Capacity-Constrained Subproblem

The iterative scheme of equation (13) based on GLP projection is outlined in the following:$$\begin{array}{}\text{(19)}& {\mathbf{\lambda }}^{k+1}={\text{Proj}}_{\mathrm{\Lambda }}{\mathbf{\lambda }}^k+{\alpha }^k{\mathbf{g}}^k,k=1,2\dots ,\end{array}$$where ${\alpha }^k$ represents the step-size in the $k\text{th}$ iteration. Generally, ${\alpha }^k$ is initialized with a moderate initial value ($0<{\alpha }^0<1$) and is adjusted to increase or decrease only when there is evidence that it is too small or too large to allow the solution to progressively satisfy the side constraints at a reasonable rate. However, this provokes the determination of whether the current solution satisfies the side constraints.

So far, an intuitive, practical, and uniform judgment criterion has yet to emerge in the literature dealing with the CTAP/CNaTAP. In this paper, we consider that equation (9) describes an NCP, and the NCP can be transformed into an equivalent unconstrained optimization problem through the Fischer–Burmeister function as a gap function as follows [5]:$$\begin{array}{c}\text{(20)}& \min G\mathbf{h}={\displaystyle \sum }_a\phi {\lambda }_a,-g_a=\frac{1}{2}{{\displaystyle \sum }_a\sqrt{{{\lambda }_a}^2+{-g_a}^2}-{\lambda }_a-g_a}^2,\end{array}$$where $\mathbf{h}={\mathbf{\lambda },\mathbf{g}}^{\mathrm{T}}$, $\min G\mathbf{h}=0$ is a global minimum, so the distance of $G\mathbf{h}$ from 0 can be evaluated to determine whether the current solution satisfies the side constraints.

In this paper, ${\alpha }^k$ is updated as follows:$$\begin{array}{}\text{(21)}& {\alpha }^k=\begin{array}{cc}\begin{array}{cc}\frac{{\alpha }^{k-1}}{{\gamma }_1},& \text{if\hspace{0.17em}}G{\mathbf{h}}^k>\sigma G{\mathbf{h}}^{k-1},\\ \frac{{\alpha }^{k-1}}{{\gamma }_2},& \text{otherwise},\end{array}& \text{if\hspace{0.17em}}G{\mathbf{h}}^k>G{\mathbf{h}}^{k-1},\\ \begin{array}{cc}{\alpha }^{k-1}\times {\gamma }_1,& \text{if\hspace{0.17em}}0<\delta <1\times {10}^{-5}\text{and\hspace{0.17em}MaxDiff}<1\times {10}^{-2},\\ {\alpha }^{k-1}\times {\gamma }_2,& \text{if\hspace{0.17em}}1\times {10}^{-5}<\delta <1\times {10}^{-2}\text{and\hspace{0.17em}MaxDiff}<1\times {10}^{-2},\end{array}& \text{if\hspace{0.17em}}G{\mathbf{h}}^k\le G{\mathbf{h}}^{k-1},\end{array}\end{array}$$where $\sigma >1,{\gamma }_1>{\gamma }_2>1,$$\delta =G{\mathbf{h}}^{k-1}-G{\mathbf{h}}^k/G{\mathbf{h}}^{k-1},$ and $\text{MaxDiff}$ represents the equilibrium convergence criterion (refer to Section 4.3). The conditions for increasing ${\alpha }^k$ are quite strict to avoid unnecessary oscillations.

3.3. Algorithm Details

Efficiently solving CNaTAP necessitates the design of an appropriate algorithmic framework. Figure 1 depicts three potential schemes: Plan A involves updating the Lagrange multiplier $\mathbf{\lambda }$ once, followed by multiple column generations and flow shifts until the generalized equilibrium achieves an adequate accuracy under the current $\mathbf{\lambda }$ or the maximum number of iterations is reached before $\mathbf{\lambda }$ is updated again. Essentially, this approach integrates a $\mathbf{\lambda }$-updating step at the outermost level of the NaTAP solution algorithmic framework, akin to the framework adopted by Feng et al. [24]. Plan B entails simultaneous $\mathbf{\lambda }$ updates and column generation, followed by multiple traffic shifts until the generalized equilibrium reaches an adequate accuracy or the maximum number of iterations is attained, after which $\mathbf{\lambda }$ updates and column generation resume. Plan C involves performing one column generation, followed by multiple $\mathbf{\lambda }$ updates and traffic shifts until the generalized equilibrium meets an adequate accuracy or reaches the maximum number of iterations before proceeding to the next column generation. According to the results of our numerical experiments (refer to Section 5.2), Plan C generally outperforms the other plans and is thus implemented in our approach. The specific details of solving the CNaTAP are outlined in Algorithm 1.

[figure(s) omitted; refer to PDF]

Algorithm 1: The solution algorithm for the CNaTAP.

1.  Initialize: lines 2–11.

2.  Set parameters: ${\epsilon }_1>0,{\epsilon }_2>0,{\alpha }^0>0,v_a=0,{\lambda }_a=0.$

3.  For each OD pair do

4.   Find the path with minimum generalized cost ${\overline{P}}^{rs}$ from $r$ to $s$.

5.   Assign all $q^{rs}$ to ${\overline{P}}^{rs}$: $f_{{\overline{P}}^{rs}}^{rs}=q^{rs}.$

6.   Push ${\overline{P}}^{rs}$ into initial working path set $P^{rs}$.

7.   Update $v_a,g_a,$ and compute $G\mathbf{h}.$

8.  Compute $\mathbf{\lambda }$ according to (19), and update ${\tilde{\eta }}_p^{rs}.$

9.  Set main loop counter: $M=0.$

10.  Initialize path flow solution vector: ${\mathbf{f}}^{\mathbf{M}}=\dots f_{{\overline{P}}^{rs}}^{rs}\dots .$

11.  Set ${\beta }^0>0$ for each OD pair.

12.  Main Loop: lines 13–32.

13.  While $\text{MaxDiff}\ge {\epsilon }_1\text{\hspace{0.17em}or\hspace{0.17em}}G\mathbf{h}\ge {\epsilon }_2$ do

14.   $M=M+1.$

15.   Remove unused paths.

16.   For each OD pair do

17.    Find the path with minimum generalized cost ${\overline{P}}^{rs}$ from $r$ to $s$.

18.    Push ${\overline{P}}^{rs}$ into $P^{rs}\text{\hspace{0.17em}and\hspace{0.17em}}{f}_{{\overline{P}}^{rs}}^{rs}=0$ if ${\overline{P}}^{rs}\notin P^{rs};$ otherwise delete ${\overline{P}}^{rs}.$

19.   Compute convergence criterion $\text{MaxDiff}.$

20.   Inner Loop: lines 21–32.

21.   Set inner loop counter: $I=0.$

22.   While ${\text{MaxDiff}}_{rs}\ge \text{MaxDiff}/100\text{\hspace{0.17em}or\hspace{0.17em}}I<\text{MaxInnerIter}$ do

23.   $I=I+1.$

24.   For each OD pair do

25.    Shift path flow based on ${\mathbf{f}}^{I+1}={\text{Proj}}_{{\mathrm{\Omega }}^{rs}}{\mathbf{f}}^I-{\beta }^I\tilde{\mathbf{\eta }}{\mathbf{f}}^I,$

26.    Compute the step-size ${\beta }^{I+1}$ based on (18).

27.   Update $v_a,g_a.$

28.   Compute convergence criterion $G\mathbf{h}.$

29.   Compute the step-size ${\alpha }^I$ based on (21).

30.   Compute ${\mathbf{\lambda }}^{I+1}$ according to (19).

31.   Update ${\tilde{n}}_p^{rs}.$

32.   Compute ${\text{MaxDiff}}_{rs}$.

4. Implementation of the Algorithm

The meticulous implementation of algorithm is vital to ensure optimal performance. This section provides methods to identify the path with minimum generalized cost, establish convergence criteria, and project onto feasible sets.

4.1. The Algorithm for Finding the Path With Minimum Generalized Cost

The nonadditive path cost function utilized in this paper is defined as follows:$$\begin{array}{}\text{(22)}& {\eta }_p^{rs}=5\times \sqrt{{\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}{t}_a}+{\displaystyle \sum }_{a\in A}\frac{{\delta }_{pa}^{rs}{t}_a}{2}+{\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}\text{tol}{\mathrm{l}}_a,\end{array}$$where the link travel time is assumed to follow the standard BPR function:$$\begin{array}{}\text{(23)}& t_a=t_{0}1+0.15{\frac{v_a}{C_a}}^4;\end{array}$$subsequently, the generalized cost of the path can be expressed as$$\begin{array}{}\text{(24)}& {\tilde{\eta }}_p^{rs}=5\times \sqrt{{\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}{t}_a}+{\displaystyle \sum }_{a\in A}\frac{{\delta }_{pa}^{rs}{t}_a}{2}+{\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}\text{tol}{\mathrm{l}}_a+{\displaystyle \sum }_{a\in A}{\delta }_{pa}^{rs}{\lambda }_a.\end{array}$$

Determining the path with the minimum generalized cost poses a significant challenge, necessitating solving the multicommodity minimum-cost flow problem. As per equation (24), the path’s generalized cost encompasses three factors: travel time, path tolls, and additional travel cost. Ideally, if a path offers lower travel time, path tolls, and extra cost compared with other paths, it would have the minimum generalized travel cost. However, this scenario seldom occurs. A common strategy involves enumerating all nondominated paths and selecting the path with the minimum generalized cost from among them. This method, known as finding the goal path through undominated paths, is widely employed in addressing the multiobjective shortest path problem (MSP) [11, 32, 33].

Consider two paths, $p_1\text{\hspace{0.17em}and\hspace{0.17em}}{p}_2,$ with corresponding attribute vectors $A_1,B_1,C_1$ and $A_2,B_2,C_2$. Path $p_1$ dominates $p_2$ if the inequalities $A_1\le A_2,$$B_1\le B_2,$ and $C_1\le C_2$ are simultaneously satisfied and at least one of the inequalities holds strictly. According to the definition of domination, it can be reasonably deduced that the dominated path cannot have the minimum generalized travel cost. Therefore, it is imperative to devise an algorithm that identifies all nondominated paths to determine the path with the minimum generalized travel cost. Drawing inspiration from identifying the set of nondominated paths in the MSP, an exact label setup method [34, 35] is employed to pinpoint the path with the minimum generalized travel cost. Although this paper focuses exclusively on three cost attributes, the method can be readily extended to cases with multiple attributes. Algorithm 2 outlines the procedure for finding the path with the minimum generalized travel cost from a given source to all other nodes in the network.

Algorithm 2: Identify the path with the minimum generalized cost.

1.  Initialize: lines 2–4.

2.  Given origin $r\in R$, set $i=r$.

3.  Set $CT=\varnothing \text{\hspace{0.17em}and\hspace{0.17em}}N{T}_i=\varnothing \text{\hspace{0.17em}for\hspace{0.17em}all\hspace{0.17em}nodes\hspace{0.17em}}i.$

4.  Set label $l=i,0,0,0,0,0,0,\text{null},\text{null}$ and add $l\text{\hspace{0.17em}to\hspace{0.17em}}CT\text{\hspace{0.17em}and\hspace{0.17em}}N{T}_i.$

5.  Update$\mathbf{\text{CT\hspace{0.17em}and\hspace{0.17em}}}\mathbf{N}{\mathbf{T}}_{\mathbf{i}},$finding nondominated labels for all nodes: lines 6–10.

6.  While $CT\ne \varnothing ,$ do

7.   Select the label with the minimum value for the seventh element from $CT.$ Record the node $i$ to which this label $l_i$ belongs. Remove this label from $CT.$

8.   For each succeeding node $j\text{\hspace{0.17em}of\hspace{0.17em}}i,$ do

9.    Record the link between $j\text{\hspace{0.17em}and\hspace{0.17em}}i$ (denoted as $L_{ij}$) and assign a new label, $l_j=j,l_{i2}+t_{L_{ij}},5\times \sqrt{l_{j2}},l_{j2}/2,l_{i5}+tol{l}_{L_{ij}},l_{i6}+{\lambda }_{L_{ij}},l_{j3}+l_{j4}+l_{j5}+l_{j6},L_{ij},l_i.$

10.    If any labels in $N{T}_i$ do not dominate $l_j,$ add $l_j\text{\hspace{0.17em}to\hspace{0.17em}}CT\text{\hspace{0.17em}and\hspace{0.17em}}N{T}_i.$ Moreover, if the $l_j$ dominates any label in $N{T}_i$, remove the dominant label from $N{T}_i.$

11.  Obtain the generalized cost minimum paths: lines 12-13.

12.  For each node $i,i\ne r,$ do

13.  Examine each label in $N{T}_i,$ identify the label with the minimum value in its seventh element, and by backtracking through this label, obtain the path with minimum generalized travel cost from $r\text{\hspace{0.17em}to\hspace{0.17em}}i.$

An empty set ($N{T}_i$) is initially assigned to each node $i$ to store the associated labels. The initial label, the origin $r,$ is added to $N{T}_r$ and the temporary set $CT.$ Each label (denoted by $l_i=i,{l_i}_2,{l_i}_3,{l_i}_4,{l_i}_5,{l_i}_6,{l_i}_7,{l_i}_8,{l_i}_9$) comprises nine elements, representing the corresponding node $i,$ the time from the origin $r$ to node $i,$ the value of time from the origin $r$ to node $i,$ the operating cost from the origin $r$ to node $i,$ the tolls from the origin $r$ to node $i,$ the additional travel cost from the origin $r$ to node $i,$ the generalized travel cost from the origin $r$ to node $i,$ the link connecting node $i$ to the previous node, and the label of the previous node. Consider two labels, $l_1\text{\hspace{0.17em}and\hspace{0.17em}}{l}_2;$$l_1$ dominates $l_2$ if the inequalities $l_{12}\le l_{22},$$l_{15}\le l_{25}$ and $l_{16}\le l_{26}$ hold, and at least one of the inequalities is strictly true. If $CT$ is null, all the nondominated labels have been identified, and the label set at each node is entirely nondominated. Backtracking through these nondominated labels reveals all nondominated paths. Examining the seventh element in the label set of each node $i$ identifies the path with minimum generalized travel cost from $r\text{\hspace{0.17em}to\hspace{0.17em}}i.$

4.2. Computational Procedures for Projecting Onto Feasible Sets

Jayakrishnan et al. [36] introduced a method to simplify the projection operation (15) called the GP method. This approach selects the shortest path as the reference path and relaxes the constraint of demand conservation based on the flow of nonshortest paths. Consequently, only nonnegative constraints remain, greatly simplifying the complexity of projection computation. The GP algorithm has been widely adopted for projection computation in large-scale TAP [11, 37]. Notably, the feasible set ${\mathrm{\Omega }}^{rs}=f_p^{rs}\in {\mathbb{R}}_{+}^n:{{\displaystyle \sum }}_{p\in P^{rs}}{f}_p^{rs}=q^{rs}$ is a simplex, which can be computed using a more concise and efficient method [38]. We consider Algorithm 3, which computes the projection of vector $\mathbf{x}$ onto the simplex ${\mathrm{\Omega }}_{rs},$ where $\mathbf{x}$ denotes the vector $\dots ,x_i,\dots ,$$x_i=f_p^{rs}-\beta \tilde{\eta }{f}_p^{rs}$.

Algorithm 3: Simplex projection algorithm.

1.  Set loop counter: $k=0.$

2.  Set $y_i^0=x_i+1/n{q}^{rs}-{{\displaystyle \sum }}_{i=1}^n{x}_i$ for $i=1,\dots ,n.$${\mathbf{Y}}^k=\dots ,y_i^k,\dots .$

3.  If ${\mathbf{Y}}^k<0$ do

4.   Compute $W=i:y_i^k>0.$

5.  $y_i^{k+1}=\begin{array}{cc}0& \text{if\hspace{0.17em}}i\notin W\\ y_i^k+1/W{q}^{rs}-{{\displaystyle \sum }}_{i\in W}{y}_i^k& \text{if\hspace{0.17em}}i\in W\end{array}$

6.   $k=k+1.$

The computation of Algorithm 3 is straightforward and can be completed in at most $n$ iterations. Panicucci et al. [38] successfully employed this method to solve the asymmetric TAP. Since $\mathrm{\Lambda }$ is a simple nonnegative orthant, projection (19) is tractable to compute, i.e., ${\text{Proj}}_{\mathrm{\Lambda }}{\mathbf{\lambda }}^k+{\alpha }^k{\mathbf{g}}^k=\max 0,{\mathbf{\lambda }}^k+{\alpha }^k{\mathbf{g}}^k.$

4.3. Convergence Criterion

An essential aspect of algorithmic execution lies in adopting an adequate convergence criterion. In this paper, we employ the maximum–minimum cost difference bound proposed by Dial [39] as a measure of equilibrium convergence, denoted as$$\begin{array}{}\text{(25)}& \text{MaxDiff}=\underset{r\in R,s\in S}{\max }{\tilde{\eta }}_{rs}^{\max }-{\pi }^{rs},\end{array}$$where ${\tilde{\eta }}_{rs}^{\max }={\max }_{p\in P^{rs}}{\tilde{\eta }}_p^{rs}$ is the maximum generalized travel cost contained in set $P^{rs}$.

Perederieieva et al. [33] have demonstrated the effectiveness of $\text{MaxDiff}$ as a convergence criterion through numerous numerical experiments. Compared with the Relative Gap (RGAP), a commonly used convergence metric in the NaTAP solutions reveals that $\text{MaxDiff}$ only reaches $1\times {10}^{-5}$ while RGAP can achieve $1\times {10}^{-10}$ or even higher accuracies. Thus, achieving $\text{MaxDiff}<1\times {10}^{-5}$ is adequate to ensure the equilibrium convergence accuracy of the CNaTAP. Furthermore, the convergence measure for capacity side constraints can utilize the $G\mathbf{h}$ mentioned in Section 3.2. Here, the proximity of $G\mathbf{h}$ to 0 serves as an indicator of whether the current solution satisfies the capacity constraints. The capacity constraint is nearly satisfied when $G\mathbf{h}<1\times {10}^{-10}.$

In order to avoid unnecessary innerloop iterations in the later stages of the calculation, we establish both MaxInnerIter and a stopping criterion for the inner loop, ${\text{MaxDiff}}_{rs}$, which measures the precision of the equilibrium within the current path set. ${\text{MaxDiff}}_{rs}$ is defined as follows:$$\begin{array}{}\text{(26)}& {\text{MaxDiff}}_{rs}=\underset{r\in R,s\in S}{\max }{\tilde{\eta }}_{rs}^{\max }-{\tilde{\eta }}_{rs}^{\min },\end{array}$$where ${\tilde{\eta }}_{rs}^{\min }={\min }_{p\in P^{rs}}{\tilde{\eta }}_p^{rs}$ is the minimum generalized travel cost contained in set $P^{rs}$.

5. Numerical Experiments

This section utilizes four varying-sized real traffic networks to demonstrate the proposed CNaTAP mode and assess the performance of the developed algorithms. Initially, we compare the assignment results of the CNaTAP and the NaTAP on the Sioux Falls network to illustrate the capability of the CNaTAP to generate realistic results. Subsequently, we evaluate the relative performance of the three algorithmic frameworks introduced in Section 3.3 across networks of different sizes, providing a comprehensive analysis of their operational efficiency. Furthermore, we examined the impact of employing the ABB and BB step-size computation methods on the overall convergence rate of the algorithm, utilizing both the Winnipeg and Chicago Sketch networks. Finally, a sensitivity analysis is conducted to explore the impact of MaxInnerIter on algorithm efficiency.

All numerical outcomes presented in this section were generated on a 64 bit Windows 11 operating system with a Gen Intel Central Processing Unit (CPU) i5-11320H @ 3.20 GHz CPU and 16G RAM. In the implementation, we set parameters as follows: ${\epsilon }_1=1\times {10}^{-5},{\epsilon }_2=1\times {10}^{-10},$${\alpha }^0=1\times {10}^{-4},{\beta }^0=0.1,\sigma =10,{\gamma }_1=1.2,{\gamma }_2=1.02,\text{MaxInnerIter}=15.$ The four testing networks were sourced from https://github.com/bstabler/TransportationNetworks, and their characteristics are provided in Table 2. All four test networks are toll-free networks. To guarantee the existence of solutions for the CNATAP, we set OD demands equal to half the total for the Sioux Falls and Anaheim networks, and the Chicago Sketch OD demand at 0.4 multiples of the total. The capacity of each link in the Winnipeg network is 3000 units.

Table 2

Detail of test networks.

Network	Nodes	Links	Zones	OD pairs
Sioux Falls	24	76	24	528
Anaheim	416	914	38	1406
Winnipeg	1067	2836	148	4345
Chicago-Sketch	933	2950	387	93,513

5.1. Assignment Results on the Sioux Falls Network

In this section, we utilize the Sioux Falls network to validate the CNaTAP model and ascertain the efficacy of the proposed algorithm. Figures 2 and 3 represent the disparities between the NaTAP and the CNaTAP solutions. The NaTAP is solved using the algorithm developed by Tan et al. [11], resulting in the oversaturation of 26 links, with flow primarily congesting internal links. However, in practical scenarios, limited capacity prompts flow to redistribute onto uncongested paths. This phenomenon is evident in the CNaTAP solution: optional paths are fully utilized, excess flow efficiently redistributes to links 12, 17, 20, 23, 25, 26, 27, 28, 32, 41, 43, 44, 56, 59, and 60, and traffic distribution across all network links adheres to capacity constraints. This example underscores how introducing capacity constraints can significantly alter traffic distribution patterns, offering more realistic and accurate predictions.

[figure(s) omitted; refer to PDF]

Analyzing the two path equilibrium flow patterns represented by the O-D pairs (10, 14), (7, 10). Tables 3 and 4 further validate that the NaTAP and the CNaTAP adhere to the UE principle and the generalized UE principle, respectively, in their successful solutions. The number of paths and corresponding flows depicted in Tables 3 and 4 indicate that the CNaTAP either generates additional paths to accommodate OD demand or distributes flows more evenly, effectively ensuring that flows on all network links comply with capacity constraints. In addition, Table 5 summarizes the 23 active constraints derived from traffic assignment results for the links, providing corresponding Lagrange multipliers for clarity.

Table 3

The NaTAP path flows for OD pairs.

(O, D)	Path	Link sequence	$f_p^{rs}$	${\eta }_p^{rs}$
(10, 14)	1	[27, 34]	1047.058	22.715
	2	[28, 44]	2.942	22.715
(7, 10)	1	[18, 55, 48]	950.000	23.065

Table 4

The CNaTAP path flows for OD pairs.

(O, D)	Path	Link sequence	$f_p^{rs}$	${\tilde{\eta }}_p^{rs}$	${{\displaystyle \sum }}_{a\in p}{\lambda }_a$
(10, 14)	1	[27, 34]	624.871	23.760	2.716
	2	[28, 44]	425.129	23.760	0.000
(7, 10)	1	[18, 55, 48]	454.918	29.340	9.040
	2	[17, 21, 25]	495.082	29.340	0.000

Table 5

Equilibrium link flows and multipliers for the Sioux Falls network.

Link	Capacity	Flow	${\lambda }_a$
1	25900.201	1550	0
2	23403.473	3350	0
3	25900.201	1550	0
4	4958.181	2950	0
5	23403.473	3350	0
6	17110.524	6099.917	0
7	23403.473	3550.083	0
8	17110.524	6099.917	0
9	17782.794	8949.917	0
10	4908.827	2750	0
11	17782.794	8949.917	0
12	4947.995	4686.108	0
13	10,000	8651.329	0
14	4958.181	2950	0
15	4947.995	4686.108	0
16	4898.588	4898.588	9.822
17	7841.811	6035.643	0
18	23403.473	6445.479	0
19	4898.588	4898.588	9.89
20	7841.811	5935.643	0
21	5050.193	4243.015	0
22	5045.823	3655.959	0
23	10,000	8651.329	0
24	5050.193	4143.015	0
25	13915.788	12548.48	0
26	13915.788	12498.48	0
27	10,000	9382.766	0
28	13512.002	13191.68	0
29	4854.918	4854.918	8.948
30	4993.511	4194.041	0
31	4908.827	2800	0
32	10,000	9282.766	0
33	4908.827	4908.827	1.926
34	4876.508	4876.508	2.716
35	23403.473	3550.083	0
36	4908.827	4908.827	1.957
37	25900.201	5608.91	0
38	25900.201	5658.91	0
39	5091.256	5091.256	7.557
40	4876.508	4876.508	2.72
41	5127.526	4681.96	0
42	4924.791	4924.791	0.218
43	13512.002	13191.68	0
44	5127.526	4731.96	0
45	14564.753	10626.56	0
46	9599.181	9599.181	1.415
47	5045.823	3655.959	0
48	4854.918	4854.918	9.04
49	5229.91	5229.91	3.86
50	19679.897	6230.967	0
51	4993.511	4194.041	0
52	5229.91	5229.91	3.787
53	4823.951	4823.951	6.193
54	23403.473	6545.479	0
55	19679.897	6230.967	0
56	23403.473	8766.61	0
57	14564.753	10626.56	0
58	4823.951	4823.951	6.181
59	5002.608	5002.608	2.15
60	23403.473	8816.61	0
61	5002.608	5002.608	2.239
62	5059.912	2950.757	0
63	5075.697	3382.685	0
64	5059.912	2949.804	0
65	5229.91	4598.574	0
66	4885.358	4885.358	2.45
67	9599.181	9599.181	1.499
68	5075.697	3383.638	0
69	5229.91	4597.621	0
70	5000	5000	1.163
71	4924.791	4874.791	0
72	5000	5000	0.809
73	5078.508	4092.878	0
74	5091.256	5091.256	7.276
75	4885.358	4885.358	2.09
76	5078.508	4042.878	0

Note: The 23 activity constraints, with the corresponding Lagrange multipliers highlighted in bold to emphasize the key constraints.

5.2. Comparison of Three Algorithmic Frameworks

In this section, we utilize the four transportation networks outlined in Table 2 to evaluate the computational efficiency of the three algorithmic frameworks illustrated in Figure 1. To ensure a fair comparison, we maintain consistent parameter settings and employ common code whenever feasible. The convergence performance of the three algorithmic frameworks on these test networks is depicted in Figures 4 and 5.

[figure(s) omitted; refer to PDF]

Figures 4 and 5 vividly illustrate the substantial differences in convergence rates among the three algorithmic frameworks. Plan C consistently demonstrates significantly higher convergence rates than Plan B and Plan A across all test networks. Comparatively, Plan B typically requires several times longer to achieve the same convergence accuracy as Plan C and Plan A often fail to attain predetermined accuracy within a short time. As depicted in Figure 4, equilibrium convergence curves generated by Plan B typically exhibit pronounced oscillations on all test networks, necessitating longer running times to converge. Equilibrium convergence curves produced by Plan A show overlapping due to excessive amplitude and frequency of oscillations. Each column generation entails an equilibrium convergence metric calculation, with Plan A involving the most column generations, resulting in overlapping convergence curves. As shown in Figure 5, the convergence curve for the side constraints in Plan B progresses slowly over an extended period, eventually reaching the target accuracy at a slower rate. In contrast, the convergence curve for Plan A exhibits a trailing effect and fails to achieve the target accuracy. Therefore, Plan A will be excluded from further comparisons with Plan B and Plan C in subsequent numerical experiments.

In order to thoroughly investigate why Plan B requires more running time than Plan C, we decompose the total running time into four parts: (1) the time required for network data input and initialization (T1), (2) the time used for removing zero flow paths, performing column generation, and conducting convergence checks (T2), (3) the execution time of the inner loop (T3), and (4) the time for outputting results (T4). We compared the performance of Plan B and Plan C based on these running times across four test networks and plotted the results in Figure 6.

[figure(s) omitted; refer to PDF]

The results depicted in Figure 6 reveal that the impact of T1 and T4 is generally negligible compared with T2 and T3. In certain instances, T1 and T4 are minuscule and indiscernible in the images due to their tiny contribution to the overall running time. Across the four tested networks, the time required for T2 and T3 in Plan B consistently exceeds that of Plan C. The sole distinction between Plan B and Plan C lies in $\mathbf{\lambda }$ updates and column generation frequency. In Plan B, $\mathbf{\lambda }$ updating coincides with column generation, leading to a higher frequency of column generation than in Plan C. Each column generation involves identifying the shortest path tree with generalized cost, a time-consuming computation that escalates with network size. Conversely, Plan C conducts multiple $\mathbf{\lambda }$ updates after each column generation, reducing the frequency of column generation while ensuring timely $\mathbf{\lambda }$ updates based on flow shift results. Extensive research highlights the importance of timely updates of network information in accelerating algorithm convergence [1, 37]. Consequently, the total running time of Plan C significantly outperforms that of Plan B.

5.3. Effects of the ABB Method

This section analyzes the effect of the ABB step-size computation method on the overall performance of solving CNaTAP within Plan C, based on experiments conducted on the Winnipeg and Chicago-Sketch networks. Figure 7 illustrates the impact of employing BB (BB1 and BB2) and ABB step-size computation methods in the equilibrium subproblem on the overall convergence rate of the algorithm.

[figure(s) omitted; refer to PDF]

Figure 7 indicates that the ABB step-size computation method notably accelerates the overall convergence rate of the CNaTAP solving algorithm compared with BB1 and BB2. Across both the Winnipeg and Chicago-Sketch networks, the overall trends of the convergence curves generated by ABB, BB1, and BB2 are similar. However, BB1 and BB2 consistently exhibit significant delays compared with ABB.

5.4. Choosing Best MaxInnerIter for Algorithm

This section evaluates the impact of multiple MaxInnerIter values on the overall convergence rate of the algorithm within Plan C, based on tests using the Chicago-Sketch network. The aim is to determine the optimal MaxInnerIter parameter. The experimental results are shown in Figure 8.

[figure(s) omitted; refer to PDF]

The results in Figure 8 reveal that when set MaxInnerIter is 1 or 100, the convergence curves exhibit pronounced oscillations, challenging achieving the expected convergence accuracy. When the set MaxInnerIter is 50, the convergence rate of the algorithm is notably slower compared with other values. This discrepancy arises because an excessively large MaxInnerIter may result in a low frequency of column generation, thereby impeding the timely incorporation of feasible paths into the working path set or removing zero flow paths from the working path set, potentially leading to algorithm divergence. Conversely, a too-small MaxInnerIter may increase the frequency of column generation, prolonging the total running time and reducing subproblem equilibrium accuracy, risking algorithm divergence. Therefore, selecting an appropriate MaxInnerIter value is pivotal for algorithm performance.

Based on the results in Figure 8, the algorithm demonstrates rapid convergence when MaxInnerIter falls within the range of 10–20. Therefore, we set MaxInnerIter to 15 in this study, as it consistently facilitates rapid convergence across all four test networks, as evidenced in Figures 4 and 5.

6. Conclusions

This paper introduces capacity side constraints to extend the NaTAP to more general scenarios, resulting in the CNaTAP model. This model is formulated using the VI formulation and decomposed into two subproblems related to equilibrium and capacity constraints based on the KKT condition. Due to the specific characteristics of the CNaTAP feasible set, existing efficient algorithms for solving the NaTAP are not directly applicable to the CNaTAP. Initially, we treat the Lagrange multiplier of capacity constraints as the fixed cost of the link when solving the equilibrium subproblem. This measure ensures the stability of the Cartesian product structure within the feasible set and facilitates the decomposition of the OD pairs, thus enabling efficient solving of the CNaTAP in practical scale networks.

Furthermore, we designed three types of algorithmic frameworks for the CNaTAP and selected the one with the best performance for this study. In addition, we employ the ABB intelligent step-size computation method in the equilibrium subproblems and further accelerate the algorithm’s convergence rate. The numerical experimental results confirm the ability of the proposed algorithm to obtain highly accurate solutions in practical transportation networks.

Numerous research directions merit further exploration. For example, investigate the integration of the CNaTAP with road tolling to develop a rational tolling scheme to alleviate traffic congestion more effectively, extending the CNaTAP to dynamic traffic assignment.

Disclosure

A preprint has previously been published (Hu et al. [40]).

Author Contributions

Wangxin Hu: conceptualization, data curation, formal analysis, writing the original draft, reviewing and editing, visualization, and software. Zhongxiang Huang: funding acquisition, methodology, reviewing and editing, supervision, and resources. Shihao Cao: reviewing and editing, supervision, visualization, and validation.

Funding

This research was jointly supported by the National Natural Science Foundation of China (no. 51978082) and the Postgraduate Scientific Research Innovation Project of Hunan Province (CX20240760).

References

[1] Y. Nie, "A Class of Bush-Based Algorithms for the Traffic Assignment Problem," Transportation Research Part B: Methodological, vol. 44 no. 1, pp. 73-89, DOI: 10.1016/j.trb.2009.06.005, 2010.

[2] O. Perederieieva, M. Ehrgott, A. Raith, J. Y. T. Wang, "A Framework for and Empirical Study of Algorithms for Traffic Assignment," Computers & Operations Research, vol. 54, pp. 90-107, DOI: 10.1016/j.cor.2014.08.024, 2015.

[3] S. A. Gabriel, D. Bernstein, "The Traffic Equilibrium Problem with Nonadditive Path Costs," Transportation Science, vol. 31 no. 4, pp. 337-348, DOI: 10.1287/trsc.31.4.337, 1997.

[4] L. Shen, F. Wang, L. Hu, X. Lyu, H. Shao, "Rescue Vehicle Allocation Problem Based on Optimal Reliable Path under Uncertainty," Journal of Central South University, vol. 29 no. 11, pp. 3779-3792, DOI: 10.1007/s11771-022-5188-1, 2022.

[5] H. K. Lo, A. Chen, "Traffic Equilibrium Problem with Route-specific Costs: Formulation and Algorithms," Transportation Research Part B: Methodological, vol. 34 no. 6, pp. 493-513, DOI: 10.1016/S0191-2615(99)00035-1, 2000.

[6] A. Chen, Z. Zhou, "The α-reliable Mean-Excess Traffic Equilibrium Model with Stochastic Travel Times," Transportation Research Part B: Methodological, vol. 44 no. 4, pp. 493-513, DOI: 10.1016/j.trb.2009.11.003, 2010.

[7] P. Chen, Y. M. Nie, "Bicriterion Shortest Path Problem With a General Nonadditive Cost," Transportation Research Part B: Methodological, vol. 57, pp. 419-435, DOI: 10.1016/j.trb.2013.05.008, 2013.

[8] C. Sun, L. Cheng, S. Zhu, F. Han, Z. Chu, "Multi-criteria User Equilibrium Model Considering Travel Time, Travel Time Reliability and Distance," Transportation Research Part D: Transport and Environment, vol. 66,DOI: 10.1016/j.trd.2017.03.002, 2019.

[9] A. Chen, Z. Zhou, X. Xu, "A Self-Adaptive Gradient Projection Algorithm for the Nonadditive Traffic Equilibrium Problem," Computers & Operations Research, vol. 39 no. 2, pp. 127-138, DOI: 10.1016/j.cor.2011.02.018, 2012.

[10] Z. Liu, X. Chen, Q. Meng, I. Kim, "Remote Park-And-Ride Network Equilibrium Model and its Applications," Transportation Research Part B: Methodological, vol. 117, pp. 37-62, DOI: 10.1016/j.trb.2018.08.004, 2018.

[11] H. Tan, M. Du, A. Chen, "Accelerating the Gradient Projection Algorithm for Solving the Non-additive Traffic Equilibrium Problem with the Barzilai-Borwein Step Size," Computers & Operations Research, vol. 141,DOI: 10.1016/j.cor.2022.105723, 2022.

[12] J. Xie, Y. M. Nie, X. Liu, "A Greedy Path-Based Algorithm for Traffic Assignment," Transportation Research Record: Journal of the Transportation Research Board, vol. 2672 no. 48, pp. 36-44, DOI: 10.1177/0361198118774236, 2018.

[13] C. F. Daganzo, "On the Traffic Assignment Problem with Flow Dependent Costs—I," Transportation Research, vol. 11 no. 6, pp. 433-437, DOI: 10.1016/0041-1647(77)90009-0, 1977.

[14] C. F. Daganzo, "On the Traffic Assignment Problem with Flow Dependent Costs—II," Transportation Research, vol. 11 no. 6, pp. 439-441, DOI: 10.1016/0041-1647(77)90010-7, 1977.

[15] H. Yang, S. Yagar, "Traffic Assignment and Signal Control in Saturated Road Networks," Transportation Research Part A: Policy and Practice, vol. 29 no. 2, pp. 125-139, DOI: 10.1016/0965-8564(94)E0007-V, 1995.

[16] H. Yang, S. Yagar, "Traffic Assignment and Traffic Control in General Freeway-Arterial Corridor Systems," Transportation Research Part B: Methodological, vol. 28 no. 6, pp. 463-486, DOI: 10.1016/0191-2615(94)90015-9, 1994.

[17] J. N. Prashker, T. Toledo, "A Gradient Projection Algorithm for Side-Constrained Traffic Assignment," European Journal of Transport and Infrastructure Research, vol. 4 no. 2, pp. 177-193, DOI: 10.18757/ejtir.2004.4.2.4261, 2004.

[18] A. H. Shahpar, H. Z. Aashtiani, A. Babazadeh, "Dynamic Penalty Function Method for the Side Constrained Traffic Assignment Problem," Applied Mathematics and Computation, vol. 206 no. 1, pp. 332-345, DOI: 10.1016/j.amc.2008.09.014, 2008.

[19] T. Larsson, M. Patriksson, "Side Constrained Traffic Equilibrium Models—Analysis, Computation and Applications," Transportation Research Part B: Methodological, vol. 33 no. 4, pp. 233-264, DOI: 10.1016/S0191-2615(98)00024-1, 1999.

[20] T. Larsson, M. Patriksson, "An Augmented Lagrangean Dual Algorithm for Link Capacity Side Constrained Traffic Assignment Problems," Transportation Research Part B: Methodological, vol. 29 no. 6, pp. 433-455, DOI: 10.1016/0191-2615(95)00016-7, 1995.

[21] S. Lawphongpanich, "Simplicial With Truncated Dantzig–Wolfe Decomposition for Nonlinear Multicommodity Network Flow Problems With Side Constraints," Operations Research Letters, vol. 26 no. 1, pp. 33-41, DOI: 10.1016/S0167-6377(99)00059-0, 2000.

[22] T. Larsson, M. Patriksson, C. Rydergren, "A Column Generation Procedure for the Side Constrained Traffic Equilibrium Problem," Transportation Research Part B: Methodological, vol. 38 no. 1, pp. 17-38, DOI: 10.1016/S0191-2615(02)00092-9, 2004.

[23] Y. Nie, H. M. Zhang, D.-H. Lee, "Models and Algorithms for the Traffic Assignment Problem with Link Capacity Constraints," Transportation Research Part B: Methodological, vol. 38 no. 4, pp. 285-312, DOI: 10.1016/S0191-2615(03)00010-9, 2004.

[24] L. Feng, J. Xie, Y. M. Nie, X. Liu, "Efficient Algorithm for the Traffic Assignment Problem with Side Constraints," Transportation Research Record: Journal of the Transportation Research Board, vol. 2674 no. 4, pp. 129-139, DOI: 10.1177/0361198120912234, 2020.

[25] A. Chen, Z. Zhou, S. Ryu, "Modeling Physical and Environmental Side Constraints in Traffic Equilibrium Problem," International Journal of Sustainable Transportation, vol. 5 no. 3, pp. 172-197, DOI: 10.1080/15568318.2010.488277, 2011.

[26] X. Xu, A. Chen, L. Cheng, "Reformulating Environmentally Constrained Traffic Equilibrium via a Smooth Gap Function," International Journal of Sustainable Transportation, vol. 9 no. 6, pp. 419-430, DOI: 10.1080/15568318.2013.777261, 2015.

[27] D. E. Boyce, B. N. Janson, R. W. Eash, "The Effect on Equilibrium Trip Assignment of Different Link Congestion Functions," Transportation Research Part A: General, vol. 15 no. 3, pp. 223-232, DOI: 10.1016/0191-2607(81)90003-0, 1981.

[28] A. A. Goldstein, "Convex Programming in Hilbert Space," Bulletin of the American Mathematical Society, vol. 70 no. 5, pp. 709-710, DOI: 10.1090/S0002-9904-1964-11178-2, 1964.

[29] E. S. Levitin, B. T. Polyak, "Constrained Minimization Methods," USSR Computational Mathematics and Mathematical Physics, vol. 6 no. 5,DOI: 10.1016/0041-5553(66)90114-5, 1966.

[30] J. Barzilai, J. M. Borwein, "Two-point Step Size Gradient Methods," IMA Journal of Numerical Analysis, vol. 8 no. 1, pp. 141-148, DOI: 10.1093/imanum/8.1.141, 1988.

[31] B. Zhou, L. Gao, Y.-H. Dai, "Gradient Methods with Adaptive Step-Sizes," Computational Optimization and Applications, vol. 35 no. 1, pp. 69-86, DOI: 10.1007/s10589-006-6446-0, 2006.

[32] A. Chen, J. S. Oh, D. Park, W. Recker, "Solving the Bicriteria Traffic Equilibrium Problem with Variable Demand and Nonlinear Path Costs," Applied Mathematics and Computation, vol. 217 no. 7, pp. 3020-3031, DOI: 10.1016/j.amc.2010.08.035, 2010.

[33] O. Perederieieva, A. Raith, M. Schmidt, "Non-additive Shortest Path in the Context of Traffic Assignment," European Journal of Operational Research, vol. 268 no. 1, pp. 325-338, DOI: 10.1016/j.ejor.2018.01.017, 2018.

[34] S. Demeyer, J. Goedgebeur, P. Audenaert, M. Pickavet, P. Demeester, "Speeding up Martins’ Algorithm for Multiple Objective Shortest Path Problems," The Journal, vol. 11 no. 4, pp. 323-348, DOI: 10.1007/s10288-013-0232-5, 2013.

[35] E. Q. V. Martins, "On a Multicriteria Shortest Path Problem," European Journal of Operational Research, vol. 16 no. 2, pp. 236-245, DOI: 10.1016/0377-2217(84)90077-8, 1984.

[36] R. Jayakrishnan, W. T. Tsai, J. N. Prashker, S. Rajadhyaksha, "A Faster Path-Based Algorithm for Traffic Assignment," Transportation Research Record, vol. 1554, pp. 75-83, 1994.

[37] H. Zhang, Z. Liu, J. Wang, Y. Wu, "A Novel Flow Update Policy in Solving Traffic Assignment Problems: Successive over Relaxation Iteration Method," Transportation Research Part E: Logistics and Transportation Review, vol. 174,DOI: 10.1016/j.tre.2023.103111, 2023.

[38] B. Panicucci, M. Pappalardo, M. Passacantando, "A Path-Based Double Projection Method for Solving the Asymmetric Traffic Network Equilibrium Problem," Optimization Letters, vol. 1 no. 2, pp. 171-185, DOI: 10.1007/s11590-006-0002-9, 2007.

[39] R. B. Dial, "A Path-Based User-Equilibrium Traffic Assignment Algorithm that Obviates Path Storage and Enumeration," Transportation Research Part B: Methodological, vol. 40 no. 10, pp. 917-936, DOI: 10.1016/j.trb.2006.02.008, 2006.

[40] W. Hu, Z. Huang, S. Cao, "Efficient Algorithm for the Nonadditive Traffic Assignment Problem with Link Capacity Constraints," , 2024.