# Improving Truncated Newton Method for the Logit-Based Stochastic User Equilibrium Problem

Authors: Min Xu; Bojian Zhou; Jie He
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This study proposes an improved truncated Newton (ITN) method for the logit-based stochastic user equilibrium problem. The ITN method incorporates a preprocessing procedure to the traditional truncated Newton method so that a good initial point is generated, on the basis of which a useful principle is developed for the choice of the basic variables. We discuss the rationale of both improvements from a theoretical point of view and demonstrate that they can enhance the computational efficiency in the early and late iteration stages, respectively, when solving the logit-based stochastic user equilibrium problem. The ITN method is compared with other related methods in the literature. Numerical results show that the ITN method performs favorably over these methods.


---

## Body

## 1. Introduction

The main role of traffic assignment models is to forecast equilibrium link or path flows in a transportation network. These models are widely used in the field of transportation planning and network design. Traditionally, traffic assignment models are formulated as the user equilibrium (UE) or stochastic user equilibrium (SUE) problems, in which no traveler can reduce his/her actual or perceived travel time by unilaterally changing routes at equilibria [1, 2]. Among various types of traffic assignment models in the literature, the logit-based stochastic user equilibrium traffic assignment problem is most widely adopted and extensively studied [3]. This problem incorporates a random error term in the route cost function to simulate travelers’ imperfect perceptions, which follows Gumbel distribution [4]. The logit-based SUE problem can be equivalently formulated as a mathematical programming problem with a unique solution. This feature facilitates its usage in both theoretical and practical studies [5–8].The widespread application of the logit-based SUE model makes its solution approach also receive considerable interest in recent years. In general, there are two classes of solution algorithms for the logit-based SUE problem. The first class is link-based algorithms. This type of algorithms uses link flow as its variable. Since link flow is an aggregate variable of different path flows, link-based algorithms do not require explicit path enumeration. It only assumes an implicit path choice set, such as the use of all efficient paths (Dial [9]; Maher [10]), or all cyclic and acyclic paths (Bell [11]; Akamatsu [12]). The most well-known link-based algorithm is the method of successive averages (MSAs) proposed in Sheffi and Powell’s study [4]. This algorithm uses a stochastic loading procedure to produce an auxiliary link flow pattern, and the search direction equals the difference between the auxiliary link flow and the current link flow. The step size for MSA is a predetermined sequence that is decreasing towards zero, such as 1/k where k is the iteration index. Maher [10] made further modifications to the method of successive averages. In his study, the Davidon–Fletcher–Powell (DFP) method was used to generate a search direction and the cubic (or quadratic) interpolation was applied to estimate the optimal step sizes.The other class is path-based algorithms. This kind of algorithm is built on path flow variables. It requires an explicit choice of a subset of feasible paths prior to or during the assignment. Unlike the link flow variable, path flow variable is a disaggregate variable which cannot be further decomposed. Therefore, different nonlinear programming methods can be utilized in a more flexible way. For example, Damberg et al. [13] extended the disaggregate simplicial decomposition method of Larsson and Patriksson [14] to solve the logit-based SUE problem. This path-based method iteratively solves subproblems that are generated through partial linearization of the objective function. The search direction is obtained by the difference between the solution of the subproblem and the current iteration point. Bekhor and Toledo [15] proposed using the gradient projection (GP) method to solve this problem. In their study, the gradient of the objective function was projected on a linear manifold of the equality constraints, with the scaling matrix being diagonal elements of the Hessian.This study focuses on path-based algorithms for the logit SUE problem. To the best of our knowledge, almost all existing path-based algorithms have a linear or sublinear convergence rate, which is relatively slow when the iteration point is approaching the optimal solution. In order to improve the convergence, it is desirable to develop an algorithm with a superlinear convergence rate. Recently, Zhou et al. [16] proposed a modified truncated Newton (MTN) method to solve the logit-based SUE problem. This method consists of two phases. The major iteration phase is performed in the original space, while the minor iteration phase is performed in the reduced space. At each major iteration, a reduced Newton equation is approximately solved using the preconditioned conjugate gradient (PCG) method. The reduced variables in the reduced Newton equation can be changed dynamically, which facilitates the usage of the PCG method. Zhou et al. proved that the convergence rate of the MTN method is superlinear. It works very fast once the iteration point gets near to the optimal SUE solution. However, there are two important problems that are not resolved in Zhou et al.’s research. First, when the iteration point is far from the optimal SUE solution, the truncated Newton type methods are relatively slow. The reason for this phenomenon is not clear. Second, in Zhou et al.’s research, a dynamic principle on how to choose the basic route is proposed. This is only an intuitive principle. The rationale behind this principle is not explained.With the aim of addressing the above two problems, in this study, we propose an improved truncated Newton (ITN) method for the logit-based SUE problem. The ITN method makes two improvements over the traditional truncated Newton method. First, a preprocessing procedure is introduced. This procedure utilizes the partial linearization method (Patriksson [17]) to generate a good initial point in the original space. It can largely replace the early iteration stage of the traditional truncated Newton method. Second, on the basis of the generated initial point, a static principle on how to partition the coefficient matrix and the variables is developed. With this principle, the computational efficiency of the truncated Newton method in the late iteration stage can be enhanced. Furthermore, the rationale behind these two improvements is analyzed theoretically, which broadens the theoretical significance of this study.The remainder of the paper is organized as follows. Section2 outlines the traditional truncated Newton method for a linear equality constrained optimization problem. Section 3 discusses some implementation issues when applying the traditional truncated Newton method to the logit-based SUE problem. Section 4 proposes a preprocessing procedure to determine a good initial point. Section 5 develops a maximal flow principle for the choice of the basic/nonbasic variable. Numerical experiments are conducted in Section 6. Section 7 wraps up the paper with conclusions and future research directions.


## 2. The Truncated Newton Method for a Linear Equality Constrained Optimization Problem

Consider the following convex problem:(1)P1minimizefxsubject toAx=b,where f:ℜn⟶ℜ is a strictly convex function that is twice continuously differentiable, A is an m×n matrix of full row rank, and b∈ℜm. [P1] can be viewed as a general formulation of the logit-based stochastic user equilibrium problem that will be investigated in this study. We will first show how to solve [P1] by the truncated Newton method.Since [P1] only involves linear equality constraints, it can be transformed into an unconstrained optimization problem using variable reduction technique. Specifically, the matrixA and variable x are partitioned as follows:(2)A=AB,AN,x=xBxN,where AB∈ℜm×m is a nonsingular matrix and AN∈ℜm×n−m, xB∈ℜm, and xN∈ℜn−m. AB is called the basic matrix and its columns correspond to the basic variables xB. AN is called the nonbasic matrix. The columns of AN correspond to the nonbasic variables xN.Therefore, the constraintsAx=b can be rewritten as(3)ABxB+ANxN=b.By rearranging the above equation, the basic variablesxB can be expressed as follows:(4)xB=AB−1Tb−ANxN.Substituting equation (4) into [P1], we obtain the following reduced unconstrained problem:(5)P1−REDminimizexN∈ℜn−mf˜xN=fAB−1Tb−ANxNxN,where f˜xN is referred to as the reduced objective function.LetxNk be a feasible point for [P1-RED]. By approximating f˜ using a second-order Taylor series around xNk, the following subproblem can be obtained:(6)SUB−1minpN∈ℜn−mg˜kTp+12pTH˜kp,where g˜k≜∇f˜xNk and H˜k≜∇2f˜xNk are the reduced gradient and reduced Hessian of f, and p≜xN−xNk is the difference between the nonbasic variable xN and the feasible point xNk.Clearly, [SUB-1] is a quadratic programming problem. A typical method for this problem is the preconditioned conjugate gradient (PCG) method. This method constructs a sequence of conjugate directions using the objective gradient and minimizes the objective function along each of the directions. Interested readers may refer to Chapter 5 in Nocedal and Wright [18] for a detailed description of this method. It is commonly known that for large-scale optimization problems, finding the exact solution of [SUB-1] is computationally intensive. The truncated Newton method is thus designed to alleviate this drawback by solving [SUB-1] approximately if xNk is far from the optimal solution of [P1] and solving [SUB-1] more accurately when the optimal solution is approached.Letp˜k be an approximate solution of [SUB-1] generated by the PCG method. According to Lemma A.2 in Dembo and Steihaug [19], p˜k defines a descent direction with respect to the reduced objective function f˜xN. Hence, by finding an appropriate step size in this direction, the new solution point for the next iteration can be obtained.In what follows, we will give a detailed description of the truncated Newton method for the linear equality constrained optimization problem [P1]. As elaborated above, this method consists of two phases. The major iteration phase transforms the original problem into an unconstrained one and applies the truncated Newton framework to solve it. The minor iteration phase uses the PCG method to solve a quadratic programming subproblem approximately.The detailed steps of major iteration are described in Algorithm1 below. It is performed in the reduced variable space.Algorithm 1:Major iteration.
Step 0: partition the matrix A and vector x as A=AB,AN, x=xBxN, respectively. Transform [P1] into [P1-RED] according to equation (5).
Step 1: let xNk be an initial point in the reduced space. Set k=1.
Step 2: if xNk is a sufficiently accurate approximation to the minimizer of [P1-RED], terminate the algorithm.
Step 3: solve [SUB-1] approximately to generate a search direction p˜k in the reduced space.
Step 4: compute a step size λk along p˜k, for which λk produces a sufficient decrease in the reduced function f˜.
Step 5: set xNk+1=xNk+λkp˜k, k=k+1. Go to Step 2.The minor iteration is elaborated in Algorithm2 below. In each minor iteration, zj is the sequence of iterations, rj is the gradient of the objective of [SUB-1] evaluated at zj, dj is the conjugate search direction, and αj and βj are scalars that are used to determine zj and dj, respectively.Algorithm 2:Minor iteration.
Given a preconditionerMk and a forcing term ηk;
Step 1: set z0=0,r0=g˜k
SolveMkl0=r0 for l0
Setd0=−l0
Setj=0.
Step 2: set αj=rjTlj/djTH˜kdj
Setzj+1=zj+αjdj
Setrj+1=rj+αjH˜kdj
If the termination criterionrj+1/g˜k<ηk is satisfied, return p˜k=zj+1, sk=rj+1 and terminate.
Else continue with Step 3.
Step 3: solve Mklj+1=rj+1 for lj+1
Setβj+1=rj+1Tlj+1/rjTlj
Setdj+1=−lj+1+βj+1dj
Setj=j+1
Go to step 2.From the above description of the truncated Newton method, we have the following two remarks:(1)
The forcing termηk in Algorithm 2 is usually chosen to be (Nocedal and Wright [18] p168)(7)ηk=minρ,g˜k,where ρ is a given positive parameter. Clearly, ηk plays the role of controlling the solution accuracy of [SUB-1] for each major iteration. When the incumbent solution is far from the optimal solution, we have ηk=ρ>0, which means that only a few inner iterations are sufficient to satisfy the termination criterion in Step 2. When the incumbent solution is near to the optimal solution, we have ηk=g˜k⟶0, which implies that more inner iterations should be performed.(2)
In Step 2 of Algorithm2, the reduced Hessian matrix H˜k need not be formed explicitly. Algorithm 2 only requires matrix-vector products, i.e., the value of H˜kd for some vector d. In the next section, we will use this feature to simplify the calculation process of the logit-based SUE problem.Furthermore, we would like to emphasize that the truncated Newton method employed in this study is performed in a fixed reduced space. In other words, once the partitions of the matrixA and vector x are made, they remain unchanged in all major iterations. This is different from the modified truncated Newton method (MTN) proposed in Zhou et al. [16], for which the partitions can be changed from one iteration to another. Using a fixed partition simplifies both theoretical analysis and practical implementation of the algorithm. However, it puts forward higher requirements for the selection of the initial points and the basic/nonbasic variables, which is the reason why the improved truncated Newton (ITN) method is proposed. In what follows, we will first apply the traditional truncated Newton method to the logit-based SUE problem (Section 3) and then propose improvements made by the ITN method (Sections 4 and 5).


## 3. Solving the Stochastic User Equilibrium Problem Using Truncated Newton Method

### 3.1. Stochastic User Equilibrium Problem

As discussed in Section1, stochastic user equilibrium problem is fundamental to the analysis of transportation systems. It concerns the distribution of travel demands to routes in a transportation network under the assumption that travelers have different perception errors when selecting routes. This problem is defined over a transportation network GV,L, where V is the set of nodes and L is the set of directed links in the network. Let W be the set of all origin-destination (OD) pairs in the network, Rw be the set of simple (loop-free) routes between OD pair w∈W, and bw be the travel demand between OD pair w∈W. For a route r∈Rw connecting OD pair w∈W, the route flow is denoted as xrw. Let tava be the travel time on link a∈L, which is assumed to be a continuous and differentiable function of the flow va on that link only. The logit-based SUE problem can be expressed as the following minimization problem (Fisk [3]):(8)SUEminfx=∑a∈A∫0vataτdτ+1θ∑w∈W∑r∈Rwxrwlnxrw,(9)∑r∈Rwxrw=bw,∀w∈W,(10)xrw≥0,∀w∈W,r∈Rw,(11)va=∑w∈W∑r∈Rwxrwδarw,∀a∈L.In the above formulation, equation (8) is the objective function. It consists of an integral term and an entropy term. The parameter θ reflects an aggregate measure of people’s perception of travel costs. Equation (9) defines the demand/route flow conservation conditions Equation (10) indicates the nonnegativity constraints. The link-route flow relationship is characterized by equation (11), in which the indicator δarw=1 if route r between OD pair w uses link a and δarw=0 otherwise.By substituting equation (11) into equation (8), we obtain a minimization problem in terms of route flow variable xrw only. Fisk [3] has proved that the objective function f is strictly convex with respect to x, which ensures the uniqueness of equilibrium route flows. On the other hand, it is well known that the logit model assigns strictly positive flows to all paths in the choice set. Therefore, the nonnegative constraints (10) are not binding at the optimal solution. Consequently, constraints (10) and (11) can be ignored and problem [SUE] is essentially equivalent to the equality constrained minimization problem [P1].The optimal solution to the logit-based SUE problem yields an equilibrium flow distribution, which is fundamental to many transportation planning and network design problems. For practical transportation networks, in order to obtain the equilibrium flow distribution in a reasonably short time, a fast method that can cope with problem size should be employed. As is known, the truncated Newton method is one of the most accurate and fast methods for large-scale problems. In what follows, we will discuss how to use this method to solve [SUE].


### 3.2. Implementing the Truncated Newton Method to Solve the Stochastic User Equilibrium Problem

From a practical point of view, when applying the truncated Newton method to solve the logit-based SUE problem, some implementation issues should be addressed. Next, we will discuss them, respectively.


#### 3.2.1. Application of the Variable Reduction Technique

The coefficient matrixA for equation (9) can be written as follows:(12)A=1,1,...10...001,1,...1...0............00...1,1,...1.This is a block diagonal matrix whose diagonal consists of vectors of all ones (in different length). Clearly, for each row of matrixA, any variable whose coefficient is “1” can be chosen as the basic variable, and the rest variables in this row are nonbasic variables. The basic and nonbasic matrices are then formed by combining columns that correspond to the basic and nonbasic variables. It is easy to see that the basic matrix is an identity matrix, which is obviously invertible.As is known, each variable in the logit-based SUE problem represents a specific route flow. For OD pairw∈W, let xrBw denote the basic route flow variable where rB is the index of the basic route, and xrNw denote the nonbasic route flow variable where rN is the index of the nonbasic route. Therefore, we can express the basic flow variable in terms of the nonbasic flow variables, i.e.,(13)xrBw=bw−∑rN∈RNwxrNw,∀w∈W,where RNw is the set of indices of nonbasic routes between OD pair w∈W.Next, we present an example to show how the coefficient matrixA and the variable x are formed and partitioned for a specific network.Example 1.
Consider the network shown in Figure1, consisting of 4 nodes and 5 links. There are two OD pairs in the network, one is from node 1 to node 4 (indexed by OD pair 1) and the other from node 2 to node 4 (indexed by OD pair 2). The demand between each OD pair is b1=10, b2=5. There are three routes connecting OD pair 1, which are indexed by route (1, 1), route (1, 2), and route (1, 3), respectively. The two routes that connect OD pair 2 are numbered as route (2, 1) and route (2, 2).
The node sequence for each route is presented as follows:route (1, 1): node sequence 1-2-4,route (1, 2): node sequence 1-3-4,route (1, 3): node sequence 1-2-3-4,route (2, 1): node sequence 2-4,route (2, 2): node sequence 2-3-4.
Letx11, x21, x31, x12, and x12 denote the flows through the five routes; then, the demand/route flow conservation conditions for equation (9) is(14)x11+x21+x31=10,(15)x12+x22=5.
The coefficient matrixA for equations (14) and (15) is(16)A=1110000011.
If we choosex21 as the basic flow variable for OD pair 1, and x12 as the basic flow variable for OD pair 2, then the basic flow variables in equations (14) and (15) can be expressed in terms of the nonbasic flow variables, i.e.,(17)x21=10−x11−x31,(18)x12=5−x22.
Equations (17) and (18) can be written in the following matrix form:(19)1001x21x12=105−110001x11x31x22,from which we can observe that the basic and nonbasic for matrix A (c.f. equation (16)) are(20)AB=1001,AN=110001.Figure 1
An example to explain the basic/nonbasic flow variable choice.


#### 3.2.2. Computation of a Search Direction

By substituting equation (13) into equation (8), the reduced objective function can be rewritten as(21)f˜xN=∑a∈A∫0vataτdτ+1θ∑w∈W∑rN∈RNwxrNwlnxrNw+1θ∑w∈Wbw−∑rN∈RNwxrNwlnbw−∑rN∈RNwxrNw.Elements of the reduced gradient are given by(22)∂f˜∂xrNw=∑a∈Atava⋅δarNw−δarBw+1θlnxrNw−lnbw−∑rN∈RNwxrNw.Elements of the reduced Hessian matrix are(23)∂2f˜∂xrNw∂xr˜Nw˜=∑a∈A∂ta∂va⋅δar˜Nw˜−δar˜Bw˜δarNw−δarBw+1θ1xrNw⋅δrNr˜Nww˜+1bw−∑rN∈RNwxrNw⋅δww˜,where indicator δrNr˜Nww˜ is equal to 1 if rN=r˜N, w=w˜, and 0 otherwise, and indicator δww˜ is equal to 1 if w=w˜, and 0 otherwise.As mentioned before, when solving [SUB-1], we only need to compute the following matrix-vector product for some vectord:(24)u=H˜k⋅d.In equation (24), d is any vector whose dimension equals the total number of the nonbasic route flow variables.The coordinates of the vectoru can be calculated by(25)urNw=∑a∈A∂ta∂va∑w˜∈W∑r˜∈RNw˜dr˜Nw˜⋅δar˜Nw˜−∑w˜∈W∑r˜∈RNw˜dr˜Nw˜⋅δarBw˜δarNw−δarBw+1θ⋅1xrNw⋅drNw+1θ1bw−∑rN∈RNwxrNw⋅∑r˜∈RNw˜dr˜Nw,∀w∈W,r∈RNw.The diagonal elements of the reduced Hessian matrix are(26)∂2f˜∂xrNw2=∑a∈A∂ta∂va⋅δarNw−δarBw2+1θ1xrNw−1bw−∑rN∈RNwxrNw.Following the suggestion in Nash [20], the preconditioner matrix Mk is chosen as a diagonal matrix whose elements are equal to the diagonal elements of the reduced Hessian. Obviously, with Mk defined in this way, it can be easily inverted.Up to now, all the ingredients required in the computation of the search direction are available. The search direction can be obtained by iteratively performing Algorithm2.


#### 3.2.3. Determination of the Step Size

For the logit-based SUE problem, although the optimal solution satisfies the nonnegativity constraints (10), some iteration points may violate these constraints during the iterative process. When this happens, the term lnxrNw that appears in equations (21) and (22) becomes undefined. To avoid such a circumstance, at each iteration, a restriction on the step size should be imposed to ensure strictly positive path flows, i.e., the step size λk should satisfy the following two constraints:(27)xrNwk+λkp˜rNwk>0,∀w∈W,rN∈RNw,(28)∑rN∈RNwxrNwk+λkp˜rNwk<bw,∀w∈W.By incorporating the above restriction into the well-known Armijo line search rule, the step size can be determined by(29)λk=βi¯,where i¯ is the smallest nonnegative integer i which satisfies equations (27) and (28) and the following inequality:(30)f˜xNk+βip˜k≤f˜xNk+αβig˜xNkTp˜k,α,β∈0,1.


## 3.1. Stochastic User Equilibrium Problem

As discussed in Section1, stochastic user equilibrium problem is fundamental to the analysis of transportation systems. It concerns the distribution of travel demands to routes in a transportation network under the assumption that travelers have different perception errors when selecting routes. This problem is defined over a transportation network GV,L, where V is the set of nodes and L is the set of directed links in the network. Let W be the set of all origin-destination (OD) pairs in the network, Rw be the set of simple (loop-free) routes between OD pair w∈W, and bw be the travel demand between OD pair w∈W. For a route r∈Rw connecting OD pair w∈W, the route flow is denoted as xrw. Let tava be the travel time on link a∈L, which is assumed to be a continuous and differentiable function of the flow va on that link only. The logit-based SUE problem can be expressed as the following minimization problem (Fisk [3]):(8)SUEminfx=∑a∈A∫0vataτdτ+1θ∑w∈W∑r∈Rwxrwlnxrw,(9)∑r∈Rwxrw=bw,∀w∈W,(10)xrw≥0,∀w∈W,r∈Rw,(11)va=∑w∈W∑r∈Rwxrwδarw,∀a∈L.In the above formulation, equation (8) is the objective function. It consists of an integral term and an entropy term. The parameter θ reflects an aggregate measure of people’s perception of travel costs. Equation (9) defines the demand/route flow conservation conditions Equation (10) indicates the nonnegativity constraints. The link-route flow relationship is characterized by equation (11), in which the indicator δarw=1 if route r between OD pair w uses link a and δarw=0 otherwise.By substituting equation (11) into equation (8), we obtain a minimization problem in terms of route flow variable xrw only. Fisk [3] has proved that the objective function f is strictly convex with respect to x, which ensures the uniqueness of equilibrium route flows. On the other hand, it is well known that the logit model assigns strictly positive flows to all paths in the choice set. Therefore, the nonnegative constraints (10) are not binding at the optimal solution. Consequently, constraints (10) and (11) can be ignored and problem [SUE] is essentially equivalent to the equality constrained minimization problem [P1].The optimal solution to the logit-based SUE problem yields an equilibrium flow distribution, which is fundamental to many transportation planning and network design problems. For practical transportation networks, in order to obtain the equilibrium flow distribution in a reasonably short time, a fast method that can cope with problem size should be employed. As is known, the truncated Newton method is one of the most accurate and fast methods for large-scale problems. In what follows, we will discuss how to use this method to solve [SUE].


## 3.2. Implementing the Truncated Newton Method to Solve the Stochastic User Equilibrium Problem

From a practical point of view, when applying the truncated Newton method to solve the logit-based SUE problem, some implementation issues should be addressed. Next, we will discuss them, respectively.


### 3.2.1. Application of the Variable Reduction Technique

The coefficient matrixA for equation (9) can be written as follows:(12)A=1,1,...10...001,1,...1...0............00...1,1,...1.This is a block diagonal matrix whose diagonal consists of vectors of all ones (in different length). Clearly, for each row of matrixA, any variable whose coefficient is “1” can be chosen as the basic variable, and the rest variables in this row are nonbasic variables. The basic and nonbasic matrices are then formed by combining columns that correspond to the basic and nonbasic variables. It is easy to see that the basic matrix is an identity matrix, which is obviously invertible.As is known, each variable in the logit-based SUE problem represents a specific route flow. For OD pairw∈W, let xrBw denote the basic route flow variable where rB is the index of the basic route, and xrNw denote the nonbasic route flow variable where rN is the index of the nonbasic route. Therefore, we can express the basic flow variable in terms of the nonbasic flow variables, i.e.,(13)xrBw=bw−∑rN∈RNwxrNw,∀w∈W,where RNw is the set of indices of nonbasic routes between OD pair w∈W.Next, we present an example to show how the coefficient matrixA and the variable x are formed and partitioned for a specific network.Example 1.
Consider the network shown in Figure1, consisting of 4 nodes and 5 links. There are two OD pairs in the network, one is from node 1 to node 4 (indexed by OD pair 1) and the other from node 2 to node 4 (indexed by OD pair 2). The demand between each OD pair is b1=10, b2=5. There are three routes connecting OD pair 1, which are indexed by route (1, 1), route (1, 2), and route (1, 3), respectively. The two routes that connect OD pair 2 are numbered as route (2, 1) and route (2, 2).
The node sequence for each route is presented as follows:route (1, 1): node sequence 1-2-4,route (1, 2): node sequence 1-3-4,route (1, 3): node sequence 1-2-3-4,route (2, 1): node sequence 2-4,route (2, 2): node sequence 2-3-4.
Letx11, x21, x31, x12, and x12 denote the flows through the five routes; then, the demand/route flow conservation conditions for equation (9) is(14)x11+x21+x31=10,(15)x12+x22=5.
The coefficient matrixA for equations (14) and (15) is(16)A=1110000011.
If we choosex21 as the basic flow variable for OD pair 1, and x12 as the basic flow variable for OD pair 2, then the basic flow variables in equations (14) and (15) can be expressed in terms of the nonbasic flow variables, i.e.,(17)x21=10−x11−x31,(18)x12=5−x22.
Equations (17) and (18) can be written in the following matrix form:(19)1001x21x12=105−110001x11x31x22,from which we can observe that the basic and nonbasic for matrix A (c.f. equation (16)) are(20)AB=1001,AN=110001.Figure 1
An example to explain the basic/nonbasic flow variable choice.


### 3.2.2. Computation of a Search Direction

By substituting equation (13) into equation (8), the reduced objective function can be rewritten as(21)f˜xN=∑a∈A∫0vataτdτ+1θ∑w∈W∑rN∈RNwxrNwlnxrNw+1θ∑w∈Wbw−∑rN∈RNwxrNwlnbw−∑rN∈RNwxrNw.Elements of the reduced gradient are given by(22)∂f˜∂xrNw=∑a∈Atava⋅δarNw−δarBw+1θlnxrNw−lnbw−∑rN∈RNwxrNw.Elements of the reduced Hessian matrix are(23)∂2f˜∂xrNw∂xr˜Nw˜=∑a∈A∂ta∂va⋅δar˜Nw˜−δar˜Bw˜δarNw−δarBw+1θ1xrNw⋅δrNr˜Nww˜+1bw−∑rN∈RNwxrNw⋅δww˜,where indicator δrNr˜Nww˜ is equal to 1 if rN=r˜N, w=w˜, and 0 otherwise, and indicator δww˜ is equal to 1 if w=w˜, and 0 otherwise.As mentioned before, when solving [SUB-1], we only need to compute the following matrix-vector product for some vectord:(24)u=H˜k⋅d.In equation (24), d is any vector whose dimension equals the total number of the nonbasic route flow variables.The coordinates of the vectoru can be calculated by(25)urNw=∑a∈A∂ta∂va∑w˜∈W∑r˜∈RNw˜dr˜Nw˜⋅δar˜Nw˜−∑w˜∈W∑r˜∈RNw˜dr˜Nw˜⋅δarBw˜δarNw−δarBw+1θ⋅1xrNw⋅drNw+1θ1bw−∑rN∈RNwxrNw⋅∑r˜∈RNw˜dr˜Nw,∀w∈W,r∈RNw.The diagonal elements of the reduced Hessian matrix are(26)∂2f˜∂xrNw2=∑a∈A∂ta∂va⋅δarNw−δarBw2+1θ1xrNw−1bw−∑rN∈RNwxrNw.Following the suggestion in Nash [20], the preconditioner matrix Mk is chosen as a diagonal matrix whose elements are equal to the diagonal elements of the reduced Hessian. Obviously, with Mk defined in this way, it can be easily inverted.Up to now, all the ingredients required in the computation of the search direction are available. The search direction can be obtained by iteratively performing Algorithm2.


### 3.2.3. Determination of the Step Size

For the logit-based SUE problem, although the optimal solution satisfies the nonnegativity constraints (10), some iteration points may violate these constraints during the iterative process. When this happens, the term lnxrNw that appears in equations (21) and (22) becomes undefined. To avoid such a circumstance, at each iteration, a restriction on the step size should be imposed to ensure strictly positive path flows, i.e., the step size λk should satisfy the following two constraints:(27)xrNwk+λkp˜rNwk>0,∀w∈W,rN∈RNw,(28)∑rN∈RNwxrNwk+λkp˜rNwk<bw,∀w∈W.By incorporating the above restriction into the well-known Armijo line search rule, the step size can be determined by(29)λk=βi¯,where i¯ is the smallest nonnegative integer i which satisfies equations (27) and (28) and the following inequality:(30)f˜xNk+βip˜k≤f˜xNk+αβig˜xNkTp˜k,α,β∈0,1.


## 3.2.1. Application of the Variable Reduction Technique

The coefficient matrixA for equation (9) can be written as follows:(12)A=1,1,...10...001,1,...1...0............00...1,1,...1.This is a block diagonal matrix whose diagonal consists of vectors of all ones (in different length). Clearly, for each row of matrixA, any variable whose coefficient is “1” can be chosen as the basic variable, and the rest variables in this row are nonbasic variables. The basic and nonbasic matrices are then formed by combining columns that correspond to the basic and nonbasic variables. It is easy to see that the basic matrix is an identity matrix, which is obviously invertible.As is known, each variable in the logit-based SUE problem represents a specific route flow. For OD pairw∈W, let xrBw denote the basic route flow variable where rB is the index of the basic route, and xrNw denote the nonbasic route flow variable where rN is the index of the nonbasic route. Therefore, we can express the basic flow variable in terms of the nonbasic flow variables, i.e.,(13)xrBw=bw−∑rN∈RNwxrNw,∀w∈W,where RNw is the set of indices of nonbasic routes between OD pair w∈W.Next, we present an example to show how the coefficient matrixA and the variable x are formed and partitioned for a specific network.Example 1.
Consider the network shown in Figure1, consisting of 4 nodes and 5 links. There are two OD pairs in the network, one is from node 1 to node 4 (indexed by OD pair 1) and the other from node 2 to node 4 (indexed by OD pair 2). The demand between each OD pair is b1=10, b2=5. There are three routes connecting OD pair 1, which are indexed by route (1, 1), route (1, 2), and route (1, 3), respectively. The two routes that connect OD pair 2 are numbered as route (2, 1) and route (2, 2).
The node sequence for each route is presented as follows:route (1, 1): node sequence 1-2-4,route (1, 2): node sequence 1-3-4,route (1, 3): node sequence 1-2-3-4,route (2, 1): node sequence 2-4,route (2, 2): node sequence 2-3-4.
Letx11, x21, x31, x12, and x12 denote the flows through the five routes; then, the demand/route flow conservation conditions for equation (9) is(14)x11+x21+x31=10,(15)x12+x22=5.
The coefficient matrixA for equations (14) and (15) is(16)A=1110000011.
If we choosex21 as the basic flow variable for OD pair 1, and x12 as the basic flow variable for OD pair 2, then the basic flow variables in equations (14) and (15) can be expressed in terms of the nonbasic flow variables, i.e.,(17)x21=10−x11−x31,(18)x12=5−x22.
Equations (17) and (18) can be written in the following matrix form:(19)1001x21x12=105−110001x11x31x22,from which we can observe that the basic and nonbasic for matrix A (c.f. equation (16)) are(20)AB=1001,AN=110001.Figure 1
An example to explain the basic/nonbasic flow variable choice.


## 3.2.2. Computation of a Search Direction

By substituting equation (13) into equation (8), the reduced objective function can be rewritten as(21)f˜xN=∑a∈A∫0vataτdτ+1θ∑w∈W∑rN∈RNwxrNwlnxrNw+1θ∑w∈Wbw−∑rN∈RNwxrNwlnbw−∑rN∈RNwxrNw.Elements of the reduced gradient are given by(22)∂f˜∂xrNw=∑a∈Atava⋅δarNw−δarBw+1θlnxrNw−lnbw−∑rN∈RNwxrNw.Elements of the reduced Hessian matrix are(23)∂2f˜∂xrNw∂xr˜Nw˜=∑a∈A∂ta∂va⋅δar˜Nw˜−δar˜Bw˜δarNw−δarBw+1θ1xrNw⋅δrNr˜Nww˜+1bw−∑rN∈RNwxrNw⋅δww˜,where indicator δrNr˜Nww˜ is equal to 1 if rN=r˜N, w=w˜, and 0 otherwise, and indicator δww˜ is equal to 1 if w=w˜, and 0 otherwise.As mentioned before, when solving [SUB-1], we only need to compute the following matrix-vector product for some vectord:(24)u=H˜k⋅d.In equation (24), d is any vector whose dimension equals the total number of the nonbasic route flow variables.The coordinates of the vectoru can be calculated by(25)urNw=∑a∈A∂ta∂va∑w˜∈W∑r˜∈RNw˜dr˜Nw˜⋅δar˜Nw˜−∑w˜∈W∑r˜∈RNw˜dr˜Nw˜⋅δarBw˜δarNw−δarBw+1θ⋅1xrNw⋅drNw+1θ1bw−∑rN∈RNwxrNw⋅∑r˜∈RNw˜dr˜Nw,∀w∈W,r∈RNw.The diagonal elements of the reduced Hessian matrix are(26)∂2f˜∂xrNw2=∑a∈A∂ta∂va⋅δarNw−δarBw2+1θ1xrNw−1bw−∑rN∈RNwxrNw.Following the suggestion in Nash [20], the preconditioner matrix Mk is chosen as a diagonal matrix whose elements are equal to the diagonal elements of the reduced Hessian. Obviously, with Mk defined in this way, it can be easily inverted.Up to now, all the ingredients required in the computation of the search direction are available. The search direction can be obtained by iteratively performing Algorithm2.


## 3.2.3. Determination of the Step Size

For the logit-based SUE problem, although the optimal solution satisfies the nonnegativity constraints (10), some iteration points may violate these constraints during the iterative process. When this happens, the term lnxrNw that appears in equations (21) and (22) becomes undefined. To avoid such a circumstance, at each iteration, a restriction on the step size should be imposed to ensure strictly positive path flows, i.e., the step size λk should satisfy the following two constraints:(27)xrNwk+λkp˜rNwk>0,∀w∈W,rN∈RNw,(28)∑rN∈RNwxrNwk+λkp˜rNwk<bw,∀w∈W.By incorporating the above restriction into the well-known Armijo line search rule, the step size can be determined by(29)λk=βi¯,where i¯ is the smallest nonnegative integer i which satisfies equations (27) and (28) and the following inequality:(30)f˜xNk+βip˜k≤f˜xNk+αβig˜xNkTp˜k,α,β∈0,1.


## 4. A Preprocessing Procedure for the Determination of a Good Initial Point

The first important feature of the improved truncated Newton method is the incorporation of a preprocessing procedure into the traditional truncated Newton method. This procedure overcomes the drawback of the truncated Newton method and generates a good initial point to start with. It in essence replaces the early iteration stage of the truncated Newton method.


### 4.1. A Drawback of the Truncated Newton Method

The restriction in equations (27) and (28) is indispensable when applying the truncated Newton method, since it makes the reduced objective function (21) and the reduced gradient (22) well defined in all iterations. However, when the iteration point is far from the optimal SUE solution (this usually occurs in the early iteration stage), this restriction may deteriorate the performance of the algorithm. The reason is as follows, By the procedure of the Armijo rule (equations (29) and (30)), if anyone of the nonnegative constraints (10) is violated, the step size should be reduced so that each iteration is strictly feasible. Hence, in the early iteration stage, the actual step size is usually much smaller than needed, which will result in a very slow rate of convergence.However, when the iteration point is near the optimal solution (this usually occurs in the late iteration stage), the restriction in equations (27) and (28) does not work. In addition, by Theorem 2.3 in Dembo and Steihaug [19], the truncated Newton method converges superlinearly when the iteration point is within a neighbourhood of the optimal SUE solution. As a result, the performance of the truncated Newton is quite fast at the late iteration stage.The following proposition rigorously validates the above phenomenon.Proposition 1.
Letx˜k be the sequence of iterations generated by Algorithm 1. Suppose that ηk is chosen as equation (7), and the step size is chosen according to equations (27)—(30); then there exists a positive integer k¯ such that the restriction in equations (27) and (28) can be automatically satisfied for all k≥k¯.Proof.
By Theorem 2.1 in Dembo and Steihaug [19], we have(31)limk⟶∞g˜k=0,which, together with step 2 in Algorithm 2 and (7), implies(32)limk⟶∞rj+1g˜k=0.
Sincerj+1 is the gradient of the objective of [SUB-1] evaluated at p˜k, we have(33)rj+1=H˜kp˜k+g˜k.
Substituting equation(33) into equation (29) and noting that H˜k is positive definite, we obtain(34)limk⟶∞p˜k+H˜k−1g˜kg˜k=0.
Define(35)qk=p˜kg˜k,ek=g˜kg˜k.
Then equation (34) can be rewritten as(36)limk⟶0qk+H˜k−1ek=0.
SinceH˜k is positive definite and ek=1, it follows that qk is a bounded sequence. In view of equations (35) and (31), we have(37)limk⟶0p˜k=0.
Hence, for everyε>0, there exist a positive integer k¯, such that for all k≥k¯,(38)p˜k≤ε.
Equation(38) implies that for every component of p˜k,(39)prNwk≤ε,∀w∈W,rN∈RNw.
Assume that at iterationk(k≥k¯), both xrNwk and xrBwk are bounded away from zero. Then there exist δ1>0 and δ2>0 such that(40)xrNwk>δ1,xrBwk>δ2.
Sinceε can be chosen sufficiently small and λk≤1, it follows from Equations (27), (28), and (40) that(41)xrNwk+λkprNwk>δ1−λkε>0,∀w∈W,rN∈RNw,∑rN∈RNwxrNwk+λkprNwk=bw−xrBwk+λk∑rN∈RNwprNwk<bw−δ2+λkεRNw<bw,∀w∈W.
The above two inequalities suggest thatxrNwk+1 and xrBwk+1 are also bounded away from zero. Therefore, the conclusion of Proposition 1 holds true.
Furthermore, we would like to point out that in the early iteration stage of the truncated Newton method, although the subproblem (6) is only solved approximately, usually it still requires more than one inner iteration for each major iteration. This consumes more CPU times when compared with algorithms that do not need any inner iteration.
From all discussions above, if we find an initial point such that it is close to the optimal SUE solution, the drawback of the truncated Newton method can be avoided. This can be achieved by replacing the early iteration stage of the truncated Newton method with a preprocessing procedure, which will be discussed in the next subsection.


### 4.2. Preprocessing Procedure

The preprocessing procedure is proposed to find a good initial point to start with. It can largely replace the early iteration stage of the truncated Newton method. This procedure is based on the partial linearization method in Patriksson [17]. When performing this procedure, all iteration points will strictly satisfy the nonnegativity constraints, such that the step size restriction, (27) and (28), is not needed. Furthermore, at each iteration of the procedure, the generated subproblem has a closed form solution. Therefore, it does not require any inner iteration to solve this subproblem. Next, we elaborate the preprocessing procedure.From equation (8), we know that the objective function for the logit-based SUE problem is composed of two terms. Without loss of generality, we reconsider the problem [P1] in the following equivalent form:(42)P2minfx=f1x+f2xsubject toAx=b.Suppose that at iterationk, a feasible point xk is given. By approximating the first term of the objective function in [P2] with a first-order Taylor series around xk, the following subproblem can be obtained:(43)SUB−2min∇f1xx−xk+f2xsubject toAx=b.Letx¯k be an exact solution to [SUB-2]. By Theorem 2.1 in [17], it is known that if the vector x¯k−xk is nonzero, it is a descent direction with respect to the original objective function fx. The next iteration point is then obtained through a line search along this descent direction.Details of the preprocessing procedure are described in Algorithm 3.Algorithm 3:Preprocessing procedure.
Step 0: let xk be an initial feasible point in the original space. Set k=0.
Step 1: if xk is close to the minimizer of [P1], terminate this procedure.
Step 2: solve [SUB-2] to obtain a solution x¯k, then the search direction in the original space is given by pk=x¯k−xk.
Step 3: compute a step size λk=βi¯ where i¯ is the smallest nonnegative integer i such that fxk+βipk≤fxk+αβigxkTpk,α,β∈0,1.
Step 4: set xk+1=xk+λkpk, k=k+1. Go to Step 1.Clearly, by lettingf1x=∑a∈A∫0vataτdτ, f2x=1/θ∑w∈W∑r∈Rwxrwlnxrw in [P2], we obtain the logit-based SUE problem. The resulting subproblem that corresponds to [SUB-2] is(44)SUB−2−SUEmin∑w∈W∑r∈Rwcrwkxrw+1θ∑w∈W∑r∈Rwxrwlnxrw,(45)∑r∈Rwxrw=bw,∀w∈W.where crwk is the travel cost on path r∈Rw,w∈W based on the vector of path flows at iteration k.Unlike the complicated solution approach to [SUB-1], the solution to [SUB-2-SUE] in the preprocessing process can be given in closed form. The next proposition presents a detailed derivation of this solution.Proposition 2.
The subproblem [SUB-2-SUE] has a closed form solution, which can be explicitly expressed as(46)x¯rw=bwexp−θcrwk∑l∈Rwexp−θclwk,∀w∈W,r∈Rw.Proof.
Consider the following Lagrange function:(47)Lx,μ=∑w∈W∑r∈Rwcrwkxrw+1θ∑w∈W∑r∈Rwxrwlnxrw+∑w∈Wμw∑r∈Rwxrw−bw,where μw is the Lagrange multiplier associated with equation (45). Then, [SUB-2-SUE] can be transformed into the following minimization problem:(48)minLx,μ.
The first-order conditions for the above problem state that(49)∂Lx,μ∂xrw=0,∂Lx,μ∂μw=0,∀w∈W.
Hence, it follows from equation (49) that(50)∂Lx,μ∂xrw=∑w∈W∑r∈Rwcrwk+1θlnxrw+1+μw=0,∀w∈W,r∈Rw.
Solving the above equation yields(51)xrw=exp−θcrwk−1−θμw.
Insertingxrw into equation (45), we have(52)x¯rw=bwexp−θcrwk∑l∈Rwexp−θclwk,∀w∈W,r∈Rw.
This completes the proof.
In conclusion, the preprocessing procedure is an application of the partial linearization method to the logit-based SUE problem. Compared with the truncated Newton method, the preprocessing procedure has two advantages. On one hand, it is strictly feasible for all iteration points. Therefore, the step size restriction, (27) and (28), is not necessary so that the step size will not be forced to be reduced. On the other hand, it utilizes the special structure of the logit-based SUE model so that the resulting subproblem has a closed form solution. As a result, inner iterations are not needed, which can save a lot of CPU times. However, since the preprocessing procedure only uses a first-order approximation of the objective function, its convergence rate is sublinear. It may become very slow during the late iteration stage. This is why this procedure is only used to the take place of the truncated Newton method in the early iteration stage.
Convergence of the partial linearization method is established in Patriksson [17], which ensures that it is easy for the preprocessing procedure to find a good initial point. Since the preprocessing procedure is performed in the original space, it does not need to determine the basic and nonbasic variables. However, the outcome of this procedure lays a good foundation for what variable should be chosen as the basic/nonbasic variables.


## 4.1. A Drawback of the Truncated Newton Method

The restriction in equations (27) and (28) is indispensable when applying the truncated Newton method, since it makes the reduced objective function (21) and the reduced gradient (22) well defined in all iterations. However, when the iteration point is far from the optimal SUE solution (this usually occurs in the early iteration stage), this restriction may deteriorate the performance of the algorithm. The reason is as follows, By the procedure of the Armijo rule (equations (29) and (30)), if anyone of the nonnegative constraints (10) is violated, the step size should be reduced so that each iteration is strictly feasible. Hence, in the early iteration stage, the actual step size is usually much smaller than needed, which will result in a very slow rate of convergence.However, when the iteration point is near the optimal solution (this usually occurs in the late iteration stage), the restriction in equations (27) and (28) does not work. In addition, by Theorem 2.3 in Dembo and Steihaug [19], the truncated Newton method converges superlinearly when the iteration point is within a neighbourhood of the optimal SUE solution. As a result, the performance of the truncated Newton is quite fast at the late iteration stage.The following proposition rigorously validates the above phenomenon.Proposition 1.
Letx˜k be the sequence of iterations generated by Algorithm 1. Suppose that ηk is chosen as equation (7), and the step size is chosen according to equations (27)—(30); then there exists a positive integer k¯ such that the restriction in equations (27) and (28) can be automatically satisfied for all k≥k¯.Proof.
By Theorem 2.1 in Dembo and Steihaug [19], we have(31)limk⟶∞g˜k=0,which, together with step 2 in Algorithm 2 and (7), implies(32)limk⟶∞rj+1g˜k=0.
Sincerj+1 is the gradient of the objective of [SUB-1] evaluated at p˜k, we have(33)rj+1=H˜kp˜k+g˜k.
Substituting equation(33) into equation (29) and noting that H˜k is positive definite, we obtain(34)limk⟶∞p˜k+H˜k−1g˜kg˜k=0.
Define(35)qk=p˜kg˜k,ek=g˜kg˜k.
Then equation (34) can be rewritten as(36)limk⟶0qk+H˜k−1ek=0.
SinceH˜k is positive definite and ek=1, it follows that qk is a bounded sequence. In view of equations (35) and (31), we have(37)limk⟶0p˜k=0.
Hence, for everyε>0, there exist a positive integer k¯, such that for all k≥k¯,(38)p˜k≤ε.
Equation(38) implies that for every component of p˜k,(39)prNwk≤ε,∀w∈W,rN∈RNw.
Assume that at iterationk(k≥k¯), both xrNwk and xrBwk are bounded away from zero. Then there exist δ1>0 and δ2>0 such that(40)xrNwk>δ1,xrBwk>δ2.
Sinceε can be chosen sufficiently small and λk≤1, it follows from Equations (27), (28), and (40) that(41)xrNwk+λkprNwk>δ1−λkε>0,∀w∈W,rN∈RNw,∑rN∈RNwxrNwk+λkprNwk=bw−xrBwk+λk∑rN∈RNwprNwk<bw−δ2+λkεRNw<bw,∀w∈W.
The above two inequalities suggest thatxrNwk+1 and xrBwk+1 are also bounded away from zero. Therefore, the conclusion of Proposition 1 holds true.
Furthermore, we would like to point out that in the early iteration stage of the truncated Newton method, although the subproblem (6) is only solved approximately, usually it still requires more than one inner iteration for each major iteration. This consumes more CPU times when compared with algorithms that do not need any inner iteration.
From all discussions above, if we find an initial point such that it is close to the optimal SUE solution, the drawback of the truncated Newton method can be avoided. This can be achieved by replacing the early iteration stage of the truncated Newton method with a preprocessing procedure, which will be discussed in the next subsection.


## 4.2. Preprocessing Procedure

The preprocessing procedure is proposed to find a good initial point to start with. It can largely replace the early iteration stage of the truncated Newton method. This procedure is based on the partial linearization method in Patriksson [17]. When performing this procedure, all iteration points will strictly satisfy the nonnegativity constraints, such that the step size restriction, (27) and (28), is not needed. Furthermore, at each iteration of the procedure, the generated subproblem has a closed form solution. Therefore, it does not require any inner iteration to solve this subproblem. Next, we elaborate the preprocessing procedure.From equation (8), we know that the objective function for the logit-based SUE problem is composed of two terms. Without loss of generality, we reconsider the problem [P1] in the following equivalent form:(42)P2minfx=f1x+f2xsubject toAx=b.Suppose that at iterationk, a feasible point xk is given. By approximating the first term of the objective function in [P2] with a first-order Taylor series around xk, the following subproblem can be obtained:(43)SUB−2min∇f1xx−xk+f2xsubject toAx=b.Letx¯k be an exact solution to [SUB-2]. By Theorem 2.1 in [17], it is known that if the vector x¯k−xk is nonzero, it is a descent direction with respect to the original objective function fx. The next iteration point is then obtained through a line search along this descent direction.Details of the preprocessing procedure are described in Algorithm 3.Algorithm 3:Preprocessing procedure.
Step 0: let xk be an initial feasible point in the original space. Set k=0.
Step 1: if xk is close to the minimizer of [P1], terminate this procedure.
Step 2: solve [SUB-2] to obtain a solution x¯k, then the search direction in the original space is given by pk=x¯k−xk.
Step 3: compute a step size λk=βi¯ where i¯ is the smallest nonnegative integer i such that fxk+βipk≤fxk+αβigxkTpk,α,β∈0,1.
Step 4: set xk+1=xk+λkpk, k=k+1. Go to Step 1.Clearly, by lettingf1x=∑a∈A∫0vataτdτ, f2x=1/θ∑w∈W∑r∈Rwxrwlnxrw in [P2], we obtain the logit-based SUE problem. The resulting subproblem that corresponds to [SUB-2] is(44)SUB−2−SUEmin∑w∈W∑r∈Rwcrwkxrw+1θ∑w∈W∑r∈Rwxrwlnxrw,(45)∑r∈Rwxrw=bw,∀w∈W.where crwk is the travel cost on path r∈Rw,w∈W based on the vector of path flows at iteration k.Unlike the complicated solution approach to [SUB-1], the solution to [SUB-2-SUE] in the preprocessing process can be given in closed form. The next proposition presents a detailed derivation of this solution.Proposition 2.
The subproblem [SUB-2-SUE] has a closed form solution, which can be explicitly expressed as(46)x¯rw=bwexp−θcrwk∑l∈Rwexp−θclwk,∀w∈W,r∈Rw.Proof.
Consider the following Lagrange function:(47)Lx,μ=∑w∈W∑r∈Rwcrwkxrw+1θ∑w∈W∑r∈Rwxrwlnxrw+∑w∈Wμw∑r∈Rwxrw−bw,where μw is the Lagrange multiplier associated with equation (45). Then, [SUB-2-SUE] can be transformed into the following minimization problem:(48)minLx,μ.
The first-order conditions for the above problem state that(49)∂Lx,μ∂xrw=0,∂Lx,μ∂μw=0,∀w∈W.
Hence, it follows from equation (49) that(50)∂Lx,μ∂xrw=∑w∈W∑r∈Rwcrwk+1θlnxrw+1+μw=0,∀w∈W,r∈Rw.
Solving the above equation yields(51)xrw=exp−θcrwk−1−θμw.
Insertingxrw into equation (45), we have(52)x¯rw=bwexp−θcrwk∑l∈Rwexp−θclwk,∀w∈W,r∈Rw.
This completes the proof.
In conclusion, the preprocessing procedure is an application of the partial linearization method to the logit-based SUE problem. Compared with the truncated Newton method, the preprocessing procedure has two advantages. On one hand, it is strictly feasible for all iteration points. Therefore, the step size restriction, (27) and (28), is not necessary so that the step size will not be forced to be reduced. On the other hand, it utilizes the special structure of the logit-based SUE model so that the resulting subproblem has a closed form solution. As a result, inner iterations are not needed, which can save a lot of CPU times. However, since the preprocessing procedure only uses a first-order approximation of the objective function, its convergence rate is sublinear. It may become very slow during the late iteration stage. This is why this procedure is only used to the take place of the truncated Newton method in the early iteration stage.
Convergence of the partial linearization method is established in Patriksson [17], which ensures that it is easy for the preprocessing procedure to find a good initial point. Since the preprocessing procedure is performed in the original space, it does not need to determine the basic and nonbasic variables. However, the outcome of this procedure lays a good foundation for what variable should be chosen as the basic/nonbasic variables.


## 5. A Maximal Flow Principle for the Choice of the Basic/Nonbasic Variables

The second feature of the improved truncated Newton method is the development of a practical principle for choosing the basic variables in the logit-based SUE problem. This principle is applied after the preprocessing procedure. It relies on the information contained in the initial point that belongs to the original space. Such a principle accelerates the PCG method used in the minor iteration, and it in essence improves the computation efficiency of the truncated Newton method in the late iteration stage.


### 5.1. The Principle and Its Rationale

The preprocessing procedure generates an initial point that is near to the optimal SUE solution. Since this point is in the original space, it can provide us valuable information on how to optimally partition the variables and the coefficient matrix. In view of the special structure of matrix A, it is known that any variable whose coefficient is nonzero can be potentially used as the basic route flow variable. However, in practice, different choice of the basic route flow variables affects the performance of the PCG method differently and thus yields different convergence rates of the inner iteration.As discussed in Nocedal and Wright [18], the convergence behavior of the PCG method is strongly dependent on the condition number of the quadratic optimization problem (6), which is defined as(53)condH˜k=H˜kH˜k−1=σnσ1,where σn and σ1 are the largest and smallest eigenvalues of the matrix H˜k. Nocedal and Wright [18] showed that the larger the condition number, the slower the likely convergence of the PCG method. For the logit-based SUE problem, it is obvious that a different choice of the basic route flow variables yields a different H˜k. Therefore, among all possible choices, the most suitable way is to select a group of basic route variables such that the condition number condH˜k is as small as possible. However, in practice, the value of condH˜k is difficult to evaluate, which makes it hard to sort different condition numbers. Fortunately, based on the information contained in the initial point, we can at least avoid the case in which the condition number is very large. The next principle presents a strategy to avoid such a case.


#### 5.1.1. The Maximum Flow Principle

Letxrww∈W,r∈RwInt. be the initial point generated by the preprocessing procedure. The basic route can be set as the one that corresponds to the maximum flow route of each OD pair at the initial point. In other words, the basic route index rB for OD pair w should satisfy(54)rB=argmaxr∈RwxrwInt.,∀w∈W.It is worth to note that the maximum flow principle proposed in this study is static. Once the basic routes are selected, they are unchanged in all remaining iterations. In [16], a dynamic principle with a similar idea is also discussed. However, that is only an intuitive principle inspired by a simple example. The rationale behind that principle is not clear. In this study, we will fill this gap by rigorously explaining the rationale behind the maximum flow principle. Details of the explanation are given in the following two propositions.Proposition 3.
At a certain initial point, if the basic route for each OD pair is chosen as the one whose flow is nearly zero, then the condition number of the reduced Hessian matrix defined in equation (53) will be very large.Proof.
Elements of the reduced Hessian matrix are given in equation (23). This equation can be decomposed into three parts.
The first part is(55)∑a∈A∂ta∂va⋅δar˜Nw˜−δar˜Bw˜δarNw−δarBw.
The second part is(56)1θ⋅1xrNw⋅δrNr˜Nww˜.
The third part is(57)1θ⋅1bw−∑rN∈RNwxrNw⋅δww˜.
Assume that for each OD pair, we choose the basic route as the one whose flow is nearly zero. Under this assumption, the values of equations(55) and (56) are limited, whereas the value of equation (57) will be quite large. The reason is that bw−∑rN∈RNwxrNw is the flow on the basic route. By the assumption that bw−∑rN∈RNwxrNw is nearly zero, 1/bw−∑rN∈RNwxrNw will tend to infinity. Therefore, we can omit equations (55) and (56) and only consider equation (57). As a result, equation (23) approximately equals(58)∂2f˜∂yrw∂yr˜w˜≈1θ⋅1bw−∑rN∈RNwxrNw⋅δww˜.
From equation (58), we can observe that if the flow on the basic route is nearly zero, the reduced Hessian matrix H˜k can be approximately viewed as a block diagonal matrix. Elements in each block are nearly equal. Hence, there exist nearly linearly dependent columns in H˜k. This means H˜k is almost singular, and it has eigenvalues that are nearly zero. Since the reduced Hessian matrix H˜k is positive definite, all its eigenvalues are positive. As a result, the smallest eigenvalue σ1 presented in equation (53) is very small, which means that condH˜k is quite large.
From Proposition 3, we know that at the initial point, if the basic route for each OD pair is chosen as the one whose flow is away from zero, we can avoid the case of very large condition numbers. At the remaining iteration points, this conclusion still holds. Proposition 4 below establishes it theoretically.Proposition 4.
Letx0 be an initial point that is close to the optimal SUE solution. At x0, if the basic route is chosen according to the maximum flow principle, then for all remaining iteration points xk (k>0), the flow on the basic route is always bounded away from zero.Proof.
Atx0, let xrBw0 be the flow on the basic route. By the maximum flow principle, xrBw0 is the maximum variable for OD pair w∈W. Therefore, xrBw0 is bounded away from zero, and there exists δ>0 such that(59)xrBw0≥δ>0,∀w∈W.
Sincex0 is close to the optimal SUE solution x∗, we can assume that(60)x0−x∗≤ε,where ε a very small positive number.
The above inequality implies that(61)xrBw0−xrBw∗≤ε,∀w∈W.
By Theorem 2.1 in Dembo and Steihaug [19], xk converges to the SUE solution x∗ as k⟶∞, and accordingly, the basic route flow variable xrBwk also converges to xrBw∗. i.e.,(62)limk⟶∞xrBwk=xrBw∗,∀w∈W.
In view of equation (61) and the definition of limit, we have that for all k≥0,(63)xrBwk−xrBw∗≤ε.
Adding equation (61) in equation (63) yields(64)xrBw0−xrBw∗+xrBwk−xrBw∗≤2ε.
Using the triangle inequality, we have(65)xrBw0−xrBwk≤2ε.
Sinceε is very small, as a result of equations (59) and (65), we obtain that for all k>0,(66)xrBwk≥xrBw0−2ε≥δ−2ε>0.which implies that the flow on the basic route is away from zero for all k>0.
The above two propositions rigorously explain the rationality of the maximum flow principle. In Proposition 3, the special structure ofH˜k is crucial to prove the resultant assertion. In order to illuminate H˜k more clearly, we introduce an example below.


### 5.2. An Illustrative Example

Example 2.
Consider the grid network in Figure2. There are 9 nodes and 12 links in the network with only one OD pair between node 1 and node 9. The OD demand is 100. There are six different routes connecting the origin and destination.
The link travel time function is defined as follows:(67)tava=ta01+βvaCan,where ta0, Ca, tava, respectively, are link a’s free-flow travel time, capacity, and travel time with flow va, and β, n are deterministic parameters. In this example, we set β=0.6, n=4. We assume that Ca=100 for all links; for links (4, 5), (5, 6), and (7, 8), ta0=1; and for the remaining links ta0=2.
In this example, there are 6 route flow variables. The initial point is given byx0=13.93,0.12,15.09,1.46,6.57,62.83. At this point, if route i (i=1,…,6) is chosen as the basic route, the reduced Hessian matrix corresponding to this route can be formed, and its condition number can be calculated. For abbreviation, we use “the reduced Hessian matrix for route i” and “the condition number for route i” to express the above meanings.
Tables1–6 present the reduced Hessian matrices for routes 1–6. Substituting these matrices into equation (53), we can calculate the condition number for each route, which is shown in the lower left corner below each table.
By comparing these tables above, we can observe that the column vectors in Table2 are nearly linearly dependent. This coincides with the fact that route 2 is the route whose flow is nearest to 0 (f2=0.12). Clearly, the condition number for route 2 is 397.28, which is much larger than the condition number for any other routes.Figure 2
A simple grid network.Table 1
The reduced Hessian matrix for route 1.
1.399761
0.022992
0.022817
0.022992
0.012136
0.022992
0.085817
0.045225
0.074526
0.041262
0.022817
0.045225
0.159536
0.04525
0.011961
0.022992
0.074526
0.04525
0.107955
0.049282
0.012136
0.041262
0.011961
0.049282
0.051935
Condition number: 65.88.Table 2
The reduced Hessian matrix for route 2.
1.399761
1.376769
1.376944
1.376769
1.387625
1.376769
1.439594
1.399177
1.428303
1.405895
1.376944
1.399177
1.513663
1.399202
1.376769
1.376769
1.428303
1.399202
1.461731
1.413915
1.387625
1.405895
1.376769
1.413915
1.427424
Condition number: 397.28.Table 3
The reduced Hessian matrix for route 3.
0.085817
0.062825
0.040592
0.011292
0.044556
0.062825
1.439594
0.040417
0.011292
0.0337
0.040592
0.040417
0.154903
0.011316
0.011292
0.011292
0.011292
0.011316
0.04472
0.019312
0.044556
0.0337
0.011292
0.019312
0.055229
Condition number: 81.69.Table 4
The reduced Hessian matrix for route 4.
0.159536
0.136719
0.114311
0.114286
0.147575
0.136719
1.513663
0.114486
0.114461
0.136894
0.114311
0.114486
0.154903
0.143587
0.143612
0.114286
0.114461
0.143587
0.176991
0.151607
0.147575
0.136894
0.143612
0.151607
0.187549
Condition number: 88.17.Table 5
The reduced Hessian matrix for route 5.
0.107955
0.084963
0.033429
0.062705
0.058672
0.084963
1.461731
0.033429
0.06253
0.047816
0.033429
0.033429
0.04472
0.033404
0.025408
0.062705
0.06253
0.033404
0.176991
0.025383
0.058672
0.047816
0.025408
0.025383
0.061325
Condition number: 78.84.Table 6
The reduced Hessian matrix for route 6.
0.051935
0.039799
0.010673
0.039974
0.002653
0.039799
1.427424
0.021529
0.050655
0.013509
0.010673
0.021529
0.055229
0.043937
0.035917
0.039974
0.050655
0.043937
0.187549
0.035942
0.002653
0.013509
0.035917
0.035942
0.061325
Condition number: 68.24.


## 5.1. The Principle and Its Rationale

The preprocessing procedure generates an initial point that is near to the optimal SUE solution. Since this point is in the original space, it can provide us valuable information on how to optimally partition the variables and the coefficient matrix. In view of the special structure of matrix A, it is known that any variable whose coefficient is nonzero can be potentially used as the basic route flow variable. However, in practice, different choice of the basic route flow variables affects the performance of the PCG method differently and thus yields different convergence rates of the inner iteration.As discussed in Nocedal and Wright [18], the convergence behavior of the PCG method is strongly dependent on the condition number of the quadratic optimization problem (6), which is defined as(53)condH˜k=H˜kH˜k−1=σnσ1,where σn and σ1 are the largest and smallest eigenvalues of the matrix H˜k. Nocedal and Wright [18] showed that the larger the condition number, the slower the likely convergence of the PCG method. For the logit-based SUE problem, it is obvious that a different choice of the basic route flow variables yields a different H˜k. Therefore, among all possible choices, the most suitable way is to select a group of basic route variables such that the condition number condH˜k is as small as possible. However, in practice, the value of condH˜k is difficult to evaluate, which makes it hard to sort different condition numbers. Fortunately, based on the information contained in the initial point, we can at least avoid the case in which the condition number is very large. The next principle presents a strategy to avoid such a case.


### 5.1.1. The Maximum Flow Principle

Letxrww∈W,r∈RwInt. be the initial point generated by the preprocessing procedure. The basic route can be set as the one that corresponds to the maximum flow route of each OD pair at the initial point. In other words, the basic route index rB for OD pair w should satisfy(54)rB=argmaxr∈RwxrwInt.,∀w∈W.It is worth to note that the maximum flow principle proposed in this study is static. Once the basic routes are selected, they are unchanged in all remaining iterations. In [16], a dynamic principle with a similar idea is also discussed. However, that is only an intuitive principle inspired by a simple example. The rationale behind that principle is not clear. In this study, we will fill this gap by rigorously explaining the rationale behind the maximum flow principle. Details of the explanation are given in the following two propositions.Proposition 3.
At a certain initial point, if the basic route for each OD pair is chosen as the one whose flow is nearly zero, then the condition number of the reduced Hessian matrix defined in equation (53) will be very large.Proof.
Elements of the reduced Hessian matrix are given in equation (23). This equation can be decomposed into three parts.
The first part is(55)∑a∈A∂ta∂va⋅δar˜Nw˜−δar˜Bw˜δarNw−δarBw.
The second part is(56)1θ⋅1xrNw⋅δrNr˜Nww˜.
The third part is(57)1θ⋅1bw−∑rN∈RNwxrNw⋅δww˜.
Assume that for each OD pair, we choose the basic route as the one whose flow is nearly zero. Under this assumption, the values of equations(55) and (56) are limited, whereas the value of equation (57) will be quite large. The reason is that bw−∑rN∈RNwxrNw is the flow on the basic route. By the assumption that bw−∑rN∈RNwxrNw is nearly zero, 1/bw−∑rN∈RNwxrNw will tend to infinity. Therefore, we can omit equations (55) and (56) and only consider equation (57). As a result, equation (23) approximately equals(58)∂2f˜∂yrw∂yr˜w˜≈1θ⋅1bw−∑rN∈RNwxrNw⋅δww˜.
From equation (58), we can observe that if the flow on the basic route is nearly zero, the reduced Hessian matrix H˜k can be approximately viewed as a block diagonal matrix. Elements in each block are nearly equal. Hence, there exist nearly linearly dependent columns in H˜k. This means H˜k is almost singular, and it has eigenvalues that are nearly zero. Since the reduced Hessian matrix H˜k is positive definite, all its eigenvalues are positive. As a result, the smallest eigenvalue σ1 presented in equation (53) is very small, which means that condH˜k is quite large.
From Proposition 3, we know that at the initial point, if the basic route for each OD pair is chosen as the one whose flow is away from zero, we can avoid the case of very large condition numbers. At the remaining iteration points, this conclusion still holds. Proposition 4 below establishes it theoretically.Proposition 4.
Letx0 be an initial point that is close to the optimal SUE solution. At x0, if the basic route is chosen according to the maximum flow principle, then for all remaining iteration points xk (k>0), the flow on the basic route is always bounded away from zero.Proof.
Atx0, let xrBw0 be the flow on the basic route. By the maximum flow principle, xrBw0 is the maximum variable for OD pair w∈W. Therefore, xrBw0 is bounded away from zero, and there exists δ>0 such that(59)xrBw0≥δ>0,∀w∈W.
Sincex0 is close to the optimal SUE solution x∗, we can assume that(60)x0−x∗≤ε,where ε a very small positive number.
The above inequality implies that(61)xrBw0−xrBw∗≤ε,∀w∈W.
By Theorem 2.1 in Dembo and Steihaug [19], xk converges to the SUE solution x∗ as k⟶∞, and accordingly, the basic route flow variable xrBwk also converges to xrBw∗. i.e.,(62)limk⟶∞xrBwk=xrBw∗,∀w∈W.
In view of equation (61) and the definition of limit, we have that for all k≥0,(63)xrBwk−xrBw∗≤ε.
Adding equation (61) in equation (63) yields(64)xrBw0−xrBw∗+xrBwk−xrBw∗≤2ε.
Using the triangle inequality, we have(65)xrBw0−xrBwk≤2ε.
Sinceε is very small, as a result of equations (59) and (65), we obtain that for all k>0,(66)xrBwk≥xrBw0−2ε≥δ−2ε>0.which implies that the flow on the basic route is away from zero for all k>0.
The above two propositions rigorously explain the rationality of the maximum flow principle. In Proposition 3, the special structure ofH˜k is crucial to prove the resultant assertion. In order to illuminate H˜k more clearly, we introduce an example below.


## 5.1.1. The Maximum Flow Principle

Letxrww∈W,r∈RwInt. be the initial point generated by the preprocessing procedure. The basic route can be set as the one that corresponds to the maximum flow route of each OD pair at the initial point. In other words, the basic route index rB for OD pair w should satisfy(54)rB=argmaxr∈RwxrwInt.,∀w∈W.It is worth to note that the maximum flow principle proposed in this study is static. Once the basic routes are selected, they are unchanged in all remaining iterations. In [16], a dynamic principle with a similar idea is also discussed. However, that is only an intuitive principle inspired by a simple example. The rationale behind that principle is not clear. In this study, we will fill this gap by rigorously explaining the rationale behind the maximum flow principle. Details of the explanation are given in the following two propositions.Proposition 3.
At a certain initial point, if the basic route for each OD pair is chosen as the one whose flow is nearly zero, then the condition number of the reduced Hessian matrix defined in equation (53) will be very large.Proof.
Elements of the reduced Hessian matrix are given in equation (23). This equation can be decomposed into three parts.
The first part is(55)∑a∈A∂ta∂va⋅δar˜Nw˜−δar˜Bw˜δarNw−δarBw.
The second part is(56)1θ⋅1xrNw⋅δrNr˜Nww˜.
The third part is(57)1θ⋅1bw−∑rN∈RNwxrNw⋅δww˜.
Assume that for each OD pair, we choose the basic route as the one whose flow is nearly zero. Under this assumption, the values of equations(55) and (56) are limited, whereas the value of equation (57) will be quite large. The reason is that bw−∑rN∈RNwxrNw is the flow on the basic route. By the assumption that bw−∑rN∈RNwxrNw is nearly zero, 1/bw−∑rN∈RNwxrNw will tend to infinity. Therefore, we can omit equations (55) and (56) and only consider equation (57). As a result, equation (23) approximately equals(58)∂2f˜∂yrw∂yr˜w˜≈1θ⋅1bw−∑rN∈RNwxrNw⋅δww˜.
From equation (58), we can observe that if the flow on the basic route is nearly zero, the reduced Hessian matrix H˜k can be approximately viewed as a block diagonal matrix. Elements in each block are nearly equal. Hence, there exist nearly linearly dependent columns in H˜k. This means H˜k is almost singular, and it has eigenvalues that are nearly zero. Since the reduced Hessian matrix H˜k is positive definite, all its eigenvalues are positive. As a result, the smallest eigenvalue σ1 presented in equation (53) is very small, which means that condH˜k is quite large.
From Proposition 3, we know that at the initial point, if the basic route for each OD pair is chosen as the one whose flow is away from zero, we can avoid the case of very large condition numbers. At the remaining iteration points, this conclusion still holds. Proposition 4 below establishes it theoretically.Proposition 4.
Letx0 be an initial point that is close to the optimal SUE solution. At x0, if the basic route is chosen according to the maximum flow principle, then for all remaining iteration points xk (k>0), the flow on the basic route is always bounded away from zero.Proof.
Atx0, let xrBw0 be the flow on the basic route. By the maximum flow principle, xrBw0 is the maximum variable for OD pair w∈W. Therefore, xrBw0 is bounded away from zero, and there exists δ>0 such that(59)xrBw0≥δ>0,∀w∈W.
Sincex0 is close to the optimal SUE solution x∗, we can assume that(60)x0−x∗≤ε,where ε a very small positive number.
The above inequality implies that(61)xrBw0−xrBw∗≤ε,∀w∈W.
By Theorem 2.1 in Dembo and Steihaug [19], xk converges to the SUE solution x∗ as k⟶∞, and accordingly, the basic route flow variable xrBwk also converges to xrBw∗. i.e.,(62)limk⟶∞xrBwk=xrBw∗,∀w∈W.
In view of equation (61) and the definition of limit, we have that for all k≥0,(63)xrBwk−xrBw∗≤ε.
Adding equation (61) in equation (63) yields(64)xrBw0−xrBw∗+xrBwk−xrBw∗≤2ε.
Using the triangle inequality, we have(65)xrBw0−xrBwk≤2ε.
Sinceε is very small, as a result of equations (59) and (65), we obtain that for all k>0,(66)xrBwk≥xrBw0−2ε≥δ−2ε>0.which implies that the flow on the basic route is away from zero for all k>0.
The above two propositions rigorously explain the rationality of the maximum flow principle. In Proposition 3, the special structure ofH˜k is crucial to prove the resultant assertion. In order to illuminate H˜k more clearly, we introduce an example below.


## 5.2. An Illustrative Example

Example 2.
Consider the grid network in Figure2. There are 9 nodes and 12 links in the network with only one OD pair between node 1 and node 9. The OD demand is 100. There are six different routes connecting the origin and destination.
The link travel time function is defined as follows:(67)tava=ta01+βvaCan,where ta0, Ca, tava, respectively, are link a’s free-flow travel time, capacity, and travel time with flow va, and β, n are deterministic parameters. In this example, we set β=0.6, n=4. We assume that Ca=100 for all links; for links (4, 5), (5, 6), and (7, 8), ta0=1; and for the remaining links ta0=2.
In this example, there are 6 route flow variables. The initial point is given byx0=13.93,0.12,15.09,1.46,6.57,62.83. At this point, if route i (i=1,…,6) is chosen as the basic route, the reduced Hessian matrix corresponding to this route can be formed, and its condition number can be calculated. For abbreviation, we use “the reduced Hessian matrix for route i” and “the condition number for route i” to express the above meanings.
Tables1–6 present the reduced Hessian matrices for routes 1–6. Substituting these matrices into equation (53), we can calculate the condition number for each route, which is shown in the lower left corner below each table.
By comparing these tables above, we can observe that the column vectors in Table2 are nearly linearly dependent. This coincides with the fact that route 2 is the route whose flow is nearest to 0 (f2=0.12). Clearly, the condition number for route 2 is 397.28, which is much larger than the condition number for any other routes.Figure 2
A simple grid network.Table 1
The reduced Hessian matrix for route 1.
1.399761
0.022992
0.022817
0.022992
0.012136
0.022992
0.085817
0.045225
0.074526
0.041262
0.022817
0.045225
0.159536
0.04525
0.011961
0.022992
0.074526
0.04525
0.107955
0.049282
0.012136
0.041262
0.011961
0.049282
0.051935
Condition number: 65.88.Table 2
The reduced Hessian matrix for route 2.
1.399761
1.376769
1.376944
1.376769
1.387625
1.376769
1.439594
1.399177
1.428303
1.405895
1.376944
1.399177
1.513663
1.399202
1.376769
1.376769
1.428303
1.399202
1.461731
1.413915
1.387625
1.405895
1.376769
1.413915
1.427424
Condition number: 397.28.Table 3
The reduced Hessian matrix for route 3.
0.085817
0.062825
0.040592
0.011292
0.044556
0.062825
1.439594
0.040417
0.011292
0.0337
0.040592
0.040417
0.154903
0.011316
0.011292
0.011292
0.011292
0.011316
0.04472
0.019312
0.044556
0.0337
0.011292
0.019312
0.055229
Condition number: 81.69.Table 4
The reduced Hessian matrix for route 4.
0.159536
0.136719
0.114311
0.114286
0.147575
0.136719
1.513663
0.114486
0.114461
0.136894
0.114311
0.114486
0.154903
0.143587
0.143612
0.114286
0.114461
0.143587
0.176991
0.151607
0.147575
0.136894
0.143612
0.151607
0.187549
Condition number: 88.17.Table 5
The reduced Hessian matrix for route 5.
0.107955
0.084963
0.033429
0.062705
0.058672
0.084963
1.461731
0.033429
0.06253
0.047816
0.033429
0.033429
0.04472
0.033404
0.025408
0.062705
0.06253
0.033404
0.176991
0.025383
0.058672
0.047816
0.025408
0.025383
0.061325
Condition number: 78.84.Table 6
The reduced Hessian matrix for route 6.
0.051935
0.039799
0.010673
0.039974
0.002653
0.039799
1.427424
0.021529
0.050655
0.013509
0.010673
0.021529
0.055229
0.043937
0.035917
0.039974
0.050655
0.043937
0.187549
0.035942
0.002653
0.013509
0.035917
0.035942
0.061325
Condition number: 68.24.


## 6. Numerical Results

In order to numerically justify the theoretical analysis conducted in this research, this section presents some performance comparisons between the ITN method and other typical methods in the literature. These methods include the MSA method (Sheffi and Powell [4]), the MTN method (Zhou et al. [16]), and the GP method (Bekhor and Toledo [15]). Among the above four methods, the MSA method is an implicit enumeration method. It uses Dial’s STOCH procedure for network loading, which implicitly considers all efficient paths of the network. The other three methods are explicit enumeration methods, which require an explicit choice of a subset of feasible paths prior to the traffic assignment. These methods are tested on the Sioux Falls network and Winnipeg network. The Sioux Falls network is a medium-size network. It is composed of 76 links, 24 nodes, and 528 OD pairs. The Winnipeg network is a real-size network. It consists of 2836 links, 1052 nodes, and 4344 OD pairs. Both networks are taken from Bar-Gera [21].To provide a common basis for the comparison of the above three explicit enumeration methods, each of them is performed on the same working route set. We use a combination of the link elimination method (Azevedo et al. [22]) and link penalty method (De La Barra et al. [23]) to generate this working route set. For the Sioux Falls network, the average number of generated routes is 7.3 per OD pair, and the maximum number of generated routes is 11 for any OD pair. For the Winnipeg network, the average and maximum number of generated routes per OD pair are 20.3 and 29, respectively.The performances of the three explicit enumeration methods are compared in terms of iteration number and CPU times under different circumstances. For the MTN method, we only concern its major iterations, because the number of minor iterations can be inherently reflected by its CPU times. Similarly, for the ITN method, its iteration number is calculated by the sum of the iteration number of the preprocessing procedure (Algorithm3) and the major iteration number of the truncated Newton method (Algorithm 1). The convergence criterion for the ITN method is based on the root mean square error (RMSE) of the reduced gradient (Zhou et al. [24]), i.e.,(68)RMSE=g˜kH≤ε,where g˜k is the reduced gradient at the kth iteration and H is the total number of routes.In Step 2 of Algorithm1, we use ε=10−4 to terminate the ITN method. In Step 1 of Algorithm 3, we terminate the preprocessing procedure when RMSE is smaller than one-tenth of its initial value for the Sioux Falls network, and one-half of its initial value for the Winnipeg network. The parameter ρ in the forcing term (7) in Algorithm 2 is set to be 0.5.For all of the computation instances in this section, the start point is obtained by averagely assigning the demand to each route in the working route set between each OD pair. Our computer programs are coded in MATLAB and executed on a notebook computer.


### 6.1. Algorithm Performance

The algorithm performance test is carried out to compare the convergence progress of different methods at different iteration stages. In this test case, we useln(RMSE) as the convergence progress measure. The dispersion parameter θ is assumed to be 0.5.Figure3 displays the convergence performance of the explicit and implicit methods in terms of iteration numbers on the Sioux Falls and Winnipeg network. This figure in essence illustrates the convergence rate of the four methods. We find that in the early iteration stage, the convergences of the MTN and GP method are relatively slow in comparison with the ITN and MSA method on both networks. Since the MTN method belongs to the category of the truncated Newton method, this phenomenon in essence confirms the analysis developed in Section 4.1, which indicates that the step size restriction may deteriorate the convergence rate of the truncated Newton type methods in the early iteration stage. As for the GP method, it also requires a maximum step size restriction to ensure strictly positive path flows in each iteration (c.f. Bekhor and Toledo [15]). With a similar reasoning, it is also slower than the ITN method in the early iteration stage.Figure 3
Convergence performance in terms of iteration numbers. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)We can also observe from Figure3 that in the late iteration stage, the slope of the curve for the MSA method is much more gentle than the other three methods. The reason is that the predetermined step size employed in the MSA method diminishes to zero as k approaches infinity, which deteriorates its rate of convergence. Furthermore, if we compare the three explicit methods, we find that the slopes of the curve for the ITN method and the MTN method are similar, both of which are steeper than the GP method. This observation is consistent with the fact that the ITN and MTN methods are locally superlinearly convergent, while the convergence rate for GP method is only linear. Furthermore, the fast convergence rate of the ITN method in the late iteration stage also validates Proposition 1, which suggests that the restriction on the step size does not work when the iteration point is close to the optimal SUE solution.Figure4 shows the convergence performance of the explicit and implicit methods in terms of CPU times on the Sioux and Winnipeg network. This figure in essence illustrates the computational efficiency of the four methods. From this figure, we can see that the computational efficiency of the ITN method is very high, which validates the effectiveness of the maximum flow principle proposed in Section 5.1. Compared with the MTN and GP methods, the ITN method can reduce the CPU times by roughly 50%–80%. The reason for the greater efficiency of the ITN method is two-fold. First, the number of iterations required by the ITN method is smaller than that of the MTN and GP method. Secondly, the preprocessing procedure of the ITN method does not entail any inner iteration, thus saving a lot of CPU times in the early iteration stage.Figure 4
Convergence performance in terms of CPU times. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)


### 6.2. Sensitivity Analysis

By varying the value ofθ or multiplying the model demand with different factors, we can examine the impact of the dispersion parameter or demand factor on the performance of different methods. As we learned from Section 6.1, the MSA method (i.e., the implicit method) is much inferior to the other three methods, especially when computing exact solutions. In order for the curves of different methods contained the relevant figures to be comparable, in the sensitivity analysis, we will omit the MSA method and only examine the performance of the ITN method with the MTN and GP method under various conditions.


#### 6.2.1. Sensitivity Tests for Different Dispersion Parameters

The iteration numbers and CPU times by differentθ for the three methods on the Sioux Falls and Winnipeg network are illustrated in Figures 5 and 6. As can be observed from Figure 5, the ITN method is superior to the MTN and GP method in terms of iteration number for all values of θ on both networks. This suggests that the average convergence rate of the ITN method is higher than the other two methods. As for the required CPU times shown in Figure 6, we can observe the following fact. For the Sioux Falls network, the ITN method consumes less CPU times than the MTN and GP method in all θ. For the Winnipeg network, the CPU times for the three methods are similar when θ is very small (i.e., θ = 0.1), but the ITN method consumes much less CPU times than the other two methods when θ is relatively large. This indicates that the ITN method is more efficient than the MTN and GP method for most practical cases.Figure 5
Effect of the dispersion parameter on iteration numbers. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Figure 6
Effect of the dispersion parameter on CPU times. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Furthermore, we can find that the larger the dispersion parameterθ is, the more advantages the ITN method gains. The reason for this phenomenon is that if θ is large, the start point is much farther from the optimal solution than the case when the θ is small. Since the ITN method contains a preprocessing procedure which is tailored to avoid the step size restriction in the early iteration stage, the effects of this procedure will become more prominent if the start point is farther from the optimal solution.


#### 6.2.2. Sensitivity Tests for Different Levels of Demand

Figures7 and 8 present the iteration numbers and CPU times required by different demand levels for the three methods on the Sioux Falls and Winnipeg network. As illustrated in Figure 7, for the Sioux Falls network, the number of iterations of the ITN method is smaller than the MTN and GP method for all demand factors. For the Winnipeg network, when the demand level is relatively small, the ITN method requires fewer iterations in comparison with the other two methods. When the demand level is high, the performance of the ITN and that of MTN method are comparable, both of which are superior to the GP method. From Figure 8, we can observe similar results for the curves of the three methods in terms of CPU times.Figure 7
Effect of the demand factor on iteration numbers. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Figure 8
Effect of the demand factor on CPU times. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)In view of all discussions above, we conclude that the ITN method performs better than the MTN and GP method for most practical cases.


## 6.1. Algorithm Performance

The algorithm performance test is carried out to compare the convergence progress of different methods at different iteration stages. In this test case, we useln(RMSE) as the convergence progress measure. The dispersion parameter θ is assumed to be 0.5.Figure3 displays the convergence performance of the explicit and implicit methods in terms of iteration numbers on the Sioux Falls and Winnipeg network. This figure in essence illustrates the convergence rate of the four methods. We find that in the early iteration stage, the convergences of the MTN and GP method are relatively slow in comparison with the ITN and MSA method on both networks. Since the MTN method belongs to the category of the truncated Newton method, this phenomenon in essence confirms the analysis developed in Section 4.1, which indicates that the step size restriction may deteriorate the convergence rate of the truncated Newton type methods in the early iteration stage. As for the GP method, it also requires a maximum step size restriction to ensure strictly positive path flows in each iteration (c.f. Bekhor and Toledo [15]). With a similar reasoning, it is also slower than the ITN method in the early iteration stage.Figure 3
Convergence performance in terms of iteration numbers. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)We can also observe from Figure3 that in the late iteration stage, the slope of the curve for the MSA method is much more gentle than the other three methods. The reason is that the predetermined step size employed in the MSA method diminishes to zero as k approaches infinity, which deteriorates its rate of convergence. Furthermore, if we compare the three explicit methods, we find that the slopes of the curve for the ITN method and the MTN method are similar, both of which are steeper than the GP method. This observation is consistent with the fact that the ITN and MTN methods are locally superlinearly convergent, while the convergence rate for GP method is only linear. Furthermore, the fast convergence rate of the ITN method in the late iteration stage also validates Proposition 1, which suggests that the restriction on the step size does not work when the iteration point is close to the optimal SUE solution.Figure4 shows the convergence performance of the explicit and implicit methods in terms of CPU times on the Sioux and Winnipeg network. This figure in essence illustrates the computational efficiency of the four methods. From this figure, we can see that the computational efficiency of the ITN method is very high, which validates the effectiveness of the maximum flow principle proposed in Section 5.1. Compared with the MTN and GP methods, the ITN method can reduce the CPU times by roughly 50%–80%. The reason for the greater efficiency of the ITN method is two-fold. First, the number of iterations required by the ITN method is smaller than that of the MTN and GP method. Secondly, the preprocessing procedure of the ITN method does not entail any inner iteration, thus saving a lot of CPU times in the early iteration stage.Figure 4
Convergence performance in terms of CPU times. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)


## 6.2. Sensitivity Analysis

By varying the value ofθ or multiplying the model demand with different factors, we can examine the impact of the dispersion parameter or demand factor on the performance of different methods. As we learned from Section 6.1, the MSA method (i.e., the implicit method) is much inferior to the other three methods, especially when computing exact solutions. In order for the curves of different methods contained the relevant figures to be comparable, in the sensitivity analysis, we will omit the MSA method and only examine the performance of the ITN method with the MTN and GP method under various conditions.


### 6.2.1. Sensitivity Tests for Different Dispersion Parameters

The iteration numbers and CPU times by differentθ for the three methods on the Sioux Falls and Winnipeg network are illustrated in Figures 5 and 6. As can be observed from Figure 5, the ITN method is superior to the MTN and GP method in terms of iteration number for all values of θ on both networks. This suggests that the average convergence rate of the ITN method is higher than the other two methods. As for the required CPU times shown in Figure 6, we can observe the following fact. For the Sioux Falls network, the ITN method consumes less CPU times than the MTN and GP method in all θ. For the Winnipeg network, the CPU times for the three methods are similar when θ is very small (i.e., θ = 0.1), but the ITN method consumes much less CPU times than the other two methods when θ is relatively large. This indicates that the ITN method is more efficient than the MTN and GP method for most practical cases.Figure 5
Effect of the dispersion parameter on iteration numbers. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Figure 6
Effect of the dispersion parameter on CPU times. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Furthermore, we can find that the larger the dispersion parameterθ is, the more advantages the ITN method gains. The reason for this phenomenon is that if θ is large, the start point is much farther from the optimal solution than the case when the θ is small. Since the ITN method contains a preprocessing procedure which is tailored to avoid the step size restriction in the early iteration stage, the effects of this procedure will become more prominent if the start point is farther from the optimal solution.


### 6.2.2. Sensitivity Tests for Different Levels of Demand

Figures7 and 8 present the iteration numbers and CPU times required by different demand levels for the three methods on the Sioux Falls and Winnipeg network. As illustrated in Figure 7, for the Sioux Falls network, the number of iterations of the ITN method is smaller than the MTN and GP method for all demand factors. For the Winnipeg network, when the demand level is relatively small, the ITN method requires fewer iterations in comparison with the other two methods. When the demand level is high, the performance of the ITN and that of MTN method are comparable, both of which are superior to the GP method. From Figure 8, we can observe similar results for the curves of the three methods in terms of CPU times.Figure 7
Effect of the demand factor on iteration numbers. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Figure 8
Effect of the demand factor on CPU times. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)In view of all discussions above, we conclude that the ITN method performs better than the MTN and GP method for most practical cases.


## 6.2.1. Sensitivity Tests for Different Dispersion Parameters

The iteration numbers and CPU times by differentθ for the three methods on the Sioux Falls and Winnipeg network are illustrated in Figures 5 and 6. As can be observed from Figure 5, the ITN method is superior to the MTN and GP method in terms of iteration number for all values of θ on both networks. This suggests that the average convergence rate of the ITN method is higher than the other two methods. As for the required CPU times shown in Figure 6, we can observe the following fact. For the Sioux Falls network, the ITN method consumes less CPU times than the MTN and GP method in all θ. For the Winnipeg network, the CPU times for the three methods are similar when θ is very small (i.e., θ = 0.1), but the ITN method consumes much less CPU times than the other two methods when θ is relatively large. This indicates that the ITN method is more efficient than the MTN and GP method for most practical cases.Figure 5
Effect of the dispersion parameter on iteration numbers. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Figure 6
Effect of the dispersion parameter on CPU times. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Furthermore, we can find that the larger the dispersion parameterθ is, the more advantages the ITN method gains. The reason for this phenomenon is that if θ is large, the start point is much farther from the optimal solution than the case when the θ is small. Since the ITN method contains a preprocessing procedure which is tailored to avoid the step size restriction in the early iteration stage, the effects of this procedure will become more prominent if the start point is farther from the optimal solution.


## 6.2.2. Sensitivity Tests for Different Levels of Demand

Figures7 and 8 present the iteration numbers and CPU times required by different demand levels for the three methods on the Sioux Falls and Winnipeg network. As illustrated in Figure 7, for the Sioux Falls network, the number of iterations of the ITN method is smaller than the MTN and GP method for all demand factors. For the Winnipeg network, when the demand level is relatively small, the ITN method requires fewer iterations in comparison with the other two methods. When the demand level is high, the performance of the ITN and that of MTN method are comparable, both of which are superior to the GP method. From Figure 8, we can observe similar results for the curves of the three methods in terms of CPU times.Figure 7
Effect of the demand factor on iteration numbers. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)Figure 8
Effect of the demand factor on CPU times. (a) Sioux falls network. (b) Winnipeg network.
(a)
(b)In view of all discussions above, we conclude that the ITN method performs better than the MTN and GP method for most practical cases.


## 7. Conclusions and Future Research

This study investigated truncated Newton type algorithms for the logit-based SUE problem. We showed that using the traditional truncated Newton method to solve this problem may be very slow in the early iteration stage, and we thus proposed an improved truncated Newton method to overcome this drawback. We compared the ITN method with GP and MTN method on the Sioux Falls network. Numerical results validate that the ITN method can indeed improve the performance of these methods.Further research work can be undertaken in the following aspects. First, the current research concentrates on solving the logit-based SUE model. Extending the ITN method to the more general traffic assignment models, such as the combined distribution and assignment model [25], the cross-nested logit model [26], and the generalized extreme value (GEV) family of model [27] is highly anticipated. Second, when implementing the ITN method, routes are generated prior to the traffic assignment. It is desirable to incorporate a column generation scheme to the ITN method, which can generate the path set dynamically during the assignment.



---

*Source: 7313808-2019-10-09.xml*