Parameter optimization methods for estimating dynamic origin-destination trip-tables

Transpn Res.-B, Vol. 31, No. 2, pp. 141-157, 1997 0 1997 Elsevier Science Ltd All rights reserved. Printed in Great Britain 0191-2615/97 $17.00+0.00

Pergamon

PII: SO191-261S@5)OOO12-4

PARAMETER OPTIMIZATION METHODS FOR ESTIMATING DYNAMIC ORIGIN-DESTINATION TRIP-TABLES HANIF D. SHERALI, NAMITA ARORA Department of Industrial and Systems Engineering, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061-0118, U.S.A.

and ANTOINE G. HOBEIKA Department of Civil Engineering, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061-0105, U.S.A. (Received 7 September 1994; in revised form 6 February 1996)

Abstract-Dynamic origin-destination tables help in on-line control of traffic facilities and, consequently, are of significant use in alleviating traffic congestion. Such tables find useful applications in the contexts of Advanced Traffic Management Systems and Advanced Traveler Information Systems. This paper considers theestimation of split parameters that prescribe an origindestination trip-table based on dynamic information regarding entering and exiting traffic volumes through an intersection or a small freeway segment. Two models are developed and motivated for this problem, one based on a least-squares estimation approach and the other based on a least absolute norm approach. Both models enhance existing dynamic origindestination trip-table estimation models in that they also consider freeway segments having differing time-dependent transfer lags between different pairs of entrances and exits. A projected conjugate gradient scheme is employed for solving the constrained least-squares problem and is compared against a standard commercial software. The least absolute norm estimation problem is posed as a linear programming problem and is also solved using a commercial software for the sake of comparison. Computational results are presented on a set of test problems using synthetic as well as realistic simulated data, involving the determination of origindestination trip tables for both intersection and freeway scenarios, in order to demonstrate the viability of the proposed methods. These results exhibit that, unlike as reported in the literature based on previous efforts, properly designed parameter optimization methods can indeed provide accurate estimates in a real-time implementation framework. Hence, these methods provide competitive aitematives to the iterative statistical techniques that have been heretofore used because of their real-time processing capabilities, despite their inherent inaccuracies. 0 1997 Elsevier Science Ltd. All rights reserved

I. INTRODUCTION

Dynamic origindestination (O-D) trip-tables provide key information for effective on-line route guidance in the contexts of Advanced Traffic Management Systems (ATMS) and Advanced Traveler Information Systems (ATIS). Whenever events occur that cause congestion such as accidents, road maintenance work, or chemical spills, fast and efficient traffic control and diversion strategies are required. A crucial information in developing successful control and diversion strategies in such situations is the knowledge of the destinations of travelers. O-D trip-tables are one way of representing flows through a traffic facility. An GD trip-table is a two-dimensional matrix of elements whose cell values represent the number of trips made between various O-D zone pairs in a given region. For implementing reliable on-line traffic control strategies, it is advisable that the O-D trip-table for the traffic facility under consideration be updated frequently, in keeping with changes in the traffic pattern. Dynamic O-D estimation procedures address the problem of estimating such trip-tables in real time. The goal here is to estimate the proportion of traffic that flows between each I41

142

Hanif D. Sherali et al.

pair of entrances and exits of a traffic facility such as an intersection or a freeway section. These proportions are referred to by several names such as turning proportions, split parameters, or split ratios. The matrices formed by these split parameters are called O-D trip-tables or split parameter matrices. Several approaches can be found in the literature for estimating split parameter matrices for on-line traffic control. Among the competing methods are the recursive methods given by Cremer and Keller (1983) and Nihan and Davis (1989), the Kalman filtering approach introduced by Cremer and Keller (198 1) and later modified by Nihan and Davis (1987), and the cross-correlation method of Cremer and Keller (1983). A major drawback of these methods is that they do not explicitly incorporate the flow conservation restrictions on the split parameters in the form of constraints. Rather, these constraints have to be enforced later by additional normalization steps, hence inducing inaccuracies in the solution obtained. The concept of travel time delays in this context was introduced by Bell (1991), where it was assumed that the time taken by a vehicle to traverse a junction follows a certain distribution that might span a number of time intervals. The method suggested used a constrained recursive least-squares algorithm that satisfies the conservation constraints on the split parameters only partially. Also, the computational time required by this method is very high, thereby making this approach unsuitable for real-time implementation. Cremer and Keller (1987) also present an optimization approach for computing the split parameters. Such techniques are usually referred to as parameter optimization methods. In comparison with other statistical estimation techniques, they point out that their parameter optimization method gave the most accurate results, but required the maximum computational effort. It was therefore considered less appropriate for on-line control. The present research effort focuses on developing enhanced models for dynamic O-D estimation as well as for designing more efficient parameter optimization strategies that would be capable of producing accurate estimates in real time. Two models and related algorithms are suggested for estimating the split parameter matrix in real time. In the first approach, the model is posed as a constrained least-squares estimation problem which is solved by using a projected conjugate gradient scheme. In the second approach, the model is stated as an L, estimation problem which is transformed to an equivalent linear programming problem. Such an L1 estimation strategy is more appropriate than a least-squares approach whenever there exist outliers in the data. Both models include an enhancement that permits one to consider not only an intersection, but also a limited freeway section for which there might be varying, time-dependent transfer lags in traversing between each pair of entrances and exits. These model formulations are presented in Section 2, and suitable solution procedures are developed in Section 3. Section 4 presents some test results that were obtained using synthetic as well as realistic simulated traffic data. A comparison of our approach with the Kalman filter method discussed by Cremer and Keller (1987) is also provided in this section to demonstrate that the proposed projected conjugate direction procedure used to solve the constrained least-squares problem offers a very competitive alternative. Not only does this approach produce more accurate results because of its explicit treatment of the problem constraints, but also a proper algorithmic design permits one to obtain optimal solutions with comparable efficiency, easily admitting a real-time implementation capability. The linear programming Li estimation model is also evaluated with respect to efficiency and relative accuracy in providing split parameter estimates. A computational study on the effect of outliers is also included. Section 5 concludes the paper. 2. PROBLEM FORMULATION

We consider the flow of traffic through a limited facility such as an intersection or a small linear freeway section. It is assumed that traffic volume counts taken at constant intervals of time are available for each entrance and exit of the facility. These time intervals are referred to as sampling intervals or time periods. The task is to use the traITlc count

Parameter optimization methods for O-D trip-tables

143

data that are collected over a horizon of several time periods in order to estimate a split parameter matrix [&I, where /$j is the proportion of traffic entering entrance i that is estimated to leave the facility via exitj. It is important to consider small sampling intervals to capture the variations in traffic with time so as to predict accurate split parameter matrices. Accordingly, if the time required to estimate a split parameter matrix is less than the duration of one sampling interval, then the estimation procedure is considered suitable for on-line implementation. Furthermore, since we are also considering the case of a (small) freeway section, we assume the availability of data on the average number of time periods tii taken for a vehicle to traverse from entrance i to exit j. This duration might be time-dependent and pre-estimated based on historical data, and the values of rij might therefore vary from one run of the model to the next when executed in a rolling horizon fashion. More specifically, suppose that we consider a horizon of K time periods, indexed sequentially as 1,. . .,K, each period being of duration T time units, and that we have traffic counts y,(k) of vehicles that leave the facility via each exit j during the kth time interval [kT,(k + l)Z’l for k = 1; . .,K. If such an exit j is accessible from some entrance i, we would then be interested in the number of vehicles q&k-t& entering the facility at entrance i during the corresponding shifted time intervals indexed by 1--rii; . ., K-t,, where the indexing of time periods is carried consecutively over nonpositive integers as well. The split parameters or proportions &>O would then be determined so as to minimize a measure of discrepancy between xi qi(k - tii)& and y,(k) over all exits j and time periods k = 1; . .,K, where xi /?v = 1Ventrances i. Below, we adopt two such measures of discrepancy. One is the least-squares approach that minimizes the sum of squared errors, and the other is an L1 norm approach that minimizes the sum of absolute deviations. This latter measure is more suitable in cases where there might exist outliers in the data because the resulting estimates using this approach are affected to a lesser degree by the presence of such outliers (see some related test results in Section 4). The O-D trip matrix [&] obtained using either approach can then be updated in a rolling horizon framework in which, as time progresses, the horizon used for estimating the split parameters is rolled forward by at least one additional period, perhaps maintaining a fixed horizon length, although not necessarily so. Accordingly, depending on a historical profile of the traffic flow pattern (or the use of an accompanying dynamic traffic flow model for the freeway section), the values of tii can be modified, and the revised problem can then be re-solved to update the O-D trip matrix. The linear section of freeway is assumed to be sufficiently small so that its flow dynamics can be adequately modeled to yield reasonable estimates for the transfer lags tij An advanced start solution based on the solution of the previous run can be used to accelerate this re-optimization phase. Remark 1: Note that, in the above model, we have assumed that vehicles traverse from entrance i to exit j in some fixed number of time periods tii, where this transfer time might vary from one run of the model to the next in a rolling horizon framework. However, as in Bell (1991), we can consider a variability in tii within each run itself by assuming, for example, that vehicles traversing from entrance i to exit j might require rij + s time periods, where s = - 1, 0, or 1. Accordingly, we can introduce corresponding split parameters /?&O, and hence attempt to minimize the discrepancy between Ci C, q@ - rq - s)& and y,(k) over all exits j and time periods k, where cj C, piis = 1V entrances i. This would yield a model having a similar structure as that considered herein, and an identical solution methodology can be adopted, although the resulting increase in problem size might hamper its real-time implementability. As an alternative, we can simplify this enhanced model by assuming that certain externally determined parameters CYST have been estimated such that fiiis can be approximated as aij,&Vij,s. Using this identity and suitably scaling the resulting model, produces a model having the same size and structure as that considered herein. Hence, for the sake of simplicity, we will suppress the explicit consideration of variability in tii as above, and note that the models analyzed in this paper can incorporate such a variability in the foregoing described manner. 0 To formally present a mathematical formulation of a given horizon’s problem, let us define the following attributes of the traffic facility under consideration.

144

i= l,...,rn: j= 1,. . .,n: k= l,...,K:

Ji

qXk)

Hanif D. Shed

et al.

Indices for the entrances of the traffic facility. Indices for the exits of the traflic facility. Indices for the time periods in the horizon corresponding to exiting volume data that is being considered for estimation purposes. = (i: exit j is permissible for users entering at entrance i}, V entrances i= 1,. . .,m. {i: entrance i admits users that can take exit j}, V exits j= 1,. . .,n. = Duration of the discretized sampling interval or time period. =Average number of time intervals that a vehicle needs to traverse from entrance i to exit j, V&h, j = 1,. f .,n. = Traffic volume entering the Ph entrance during the kth time interval [kT,(k+ l)T), for each i= I,... ,m. (Note that we assume the availability of data for periods k = 1 - maximumjEJirij,---, K- minimumjEJirij, for each i= 1;. .,m.) -qz{k-tii)Vich,j= l;..,n, k= l,...,K. =Traffic volume leaving via the jth exit during the kth time interval [kT,(k+ 1)7’j, for each k= l,...,K,j= l;..,n. =To be determined split parameters that estimate the portion of the volume entering at entrance i that will ultimately leave through exit j, for each ich, j= 1;. .,n.

Based on the foregoing discussion, we can now model the problem of determining the split parameters /Iii so as to minimize the sum of squared deviations between the estimated total proportion of entering volumes at the various entrances that ought to leave through exit j, and the actual exiting volume, over all the exits j= 1,. . .,n and time periods k= 1,. . .,K, as the following constrained least-squares problem LSQ. LSQ: Minimize

subject to

c Bij

=

1 Vi=

l,...,m

(lb)

i= l,...,m.

(14

jeJi

/3~~0

Vjc Ji,

In matrix form, the objective function f(p) given by eqn (la) can be written as

where, for each exit j= 1,. . .,n, yj is a column vector having components y,(k), k = 1,. . .,K, bj is a column vector having components &, ich, and Qj is a KX IIjl matrix whose (k,z’)* element is q#. Note that f is strictly convex (i.e. LSQ is ‘nondegenerate’) if and only if Q,!Qj is positive definite for each j= 1,. .-, n , and this holds true if and only if for each exit j= 1,. . .,n, the columns of Qj are linearly independent. Note that each column of Qj corresponds to the flow values qiik=qr(k-tii) over the different periods k = 1,. . .,K for each entrance $4, and hence this condition in essence requires that for each exit j= 1,. . .,n, the columns of entry measurements over the different periods for the corresponding entrances should be linearly independent. For instance, these entry flows should not be proportional to a given set of flows. Under this condition, we can be guaranteed of a unique optimum to problem LSQ. (However, we do not assume this to be the case.)

145

Parameter optimization methods for C&D trip-tables

An alternative LI estimation model that minimizes the sum of absolute deviations between observed and estimated flows can be formulated in a similar fashion as stated below. As mentioned earlier, the use of such a model is more attractive than a leastsquares approach in the presence of outliers in the input data that are significantly nonconformant with respect to the flow model. In such cases, an attempt to minimize squared deviations as in eqn (la) might cause these outliers to greatly influence the solution, whereas minimizing absolute deviations might have a lesser effect on the solution. (We illustrate this phenomenon computationally in Section 4.) Moreover, the L, estimation model is readily able to accommodate additional side-constraints that further control the determined split parameter estimates. (However, side-constraints of this type might render the problem infeasible in the presence of significant outliers, and hence should be avoided in such contexts.) More specifically, let us define yj as the average volume leaving through exit j over all the periods k = 1; . .,K, and let qii represent the average volume that enters the traffic facility at entrance i during the corresponding shifted time periods (k-tti), for k= l,..., K. Hence, we have 1 K r/=zCyj(k)

and

~~=~~qi(k-q)--J-~ K k=l

k=l

K

(2) k=l

qok*

We can now require that in addition to eqn (1b,c), while there might exist a variability between the observed exiting volumes and the flow model estimates over the different periods, the split parameters should satisfy a set of restrictions which enforce that at least on average, these estimates are exact. Such linear constraints are readily accommodated within a linear programming approach that will be used to solve the LI estimation problem. This model, with the additional restrictions stated as eqn (3~) is specified below. Ll: Minimize K

n

k=l

j=l

cc

(34

subject to

&9(jP4=.?j I fiq?O

VjcJi,

Vi=

l,...,m,

(3b)

Vj=

l,...,n,

(34

i=l,...,m.

(34

3. SOLUTION PROCEDURES

In this section, we propose suitable solution procedures for solving problems LSQ and L, given by systems (1) and (3), respectively. These are treated in turn below. 3.1. Algorithm for problem LSQ Note that this problem seeks to minimize a convex quadratic objective function f subject to simple linear proportionality or generalized upper bounding types of constraints. There exist a host of nonlinear programming algorithms that one can possibly adopt to solve this problem (see Bazaraa et al., 1993, for example). One approach might be to use unconstrained optimization techniques in concert with the popular augmented Lagrangian penalty function method. However, this method would not adequately exploit the special, simple, linear structure of the constraints eqn (1b,c). On the other hand, this problem can be effectively solved by a (generalized) reduced gradient algorithm that employs secondorder approximations to the objective function in the direction-finding subproblem via the

Hanif D. Sherali et al.

146

celebrated BFGS quasi-Newton scheme. Commercial packages such as MINOS 5.4 and GINO adopt such a strategy, and include sophisticated techniques to ensure numerical stability and to safeguard against loss of positive definiteness in the Hessian approximations. Later we will provide results using such a package that is linked to the mathematical programming language GAMS (see Brooke et al., 1992). However, if the user does not have access to such a package or, more relevantly, if such a procedure is to be implemented as a subroutine within a larger traffic management system, it would be desirable to have an efficient algorithm that is both robust and easy to implement. A third approach might be to use an active set approach in which the problem is solved via a sequence of subproblems obtained by projecting the original problem onto suitable linear subspaces defined by the constraints. Motivated by the interest in designing an effective, yet easy-toimplement procedure, we develop below a specialized projected conjugate gradient algorithm of this type, and demonstrate computationally that it provides a competitive alternative to the statistical estimation techniques discussed in Cremer and Keller (1987). This is in contrast with the parameter optimization technique implemented by Cremer and Keller (1987) based on the so-called complex approach due to Box (1965). Like quasi-Newton methods, conjugate gradient algorithms are also conjugate direction procedures that converge in at most N iterations to an optimum for an N-variable unconstrained quadratic optimization problem (see Bazaraa et al., 1993). For our case, since we have a constrained problem, we need to implement such an approach over projected subspaces of the feasible region, as mentioned above. This is readily accomplished for problem LSQ due to the simple nature of the constraints. Figure 1 illustrates the main loop of the proposed algorithm. Below, we present details for certain specific steps required to operate this algorithm. 3.1.1. Starting solution. Typically, in a rolling horizon framework, the algorithm can use the split parameter estimates that were obtained using the previous period’s data as a starting solution, and then could update these estimates given the new period’s data. However, when starting from scratch based on data pertaining to a given set of recent periods, an initial solution could be obtained by inexactly solving the unconstrained problem via an approximate solution to the linear, first-order optimality conditions. More specifically, setting the gradient equal to zero [see eqn (Id)], we obtain the system of equations

Note that for each exit s, the foregoing system yields a separable set of equations in Bis, &Z,. If the corresponding ] Z, 1 x ] Z, ] symmetric matrix of coefficients QsfQs whose (r,i)th element is given by [C, qtikqrsk] is nonsingular, i.e. the columns of Qs are linearly independent as discussed in Section 2, then this system would yield a unique solution. However, in order to conserve effort, we solve this system approximately by aggregating the equations for each exit s, and then letting each fiis, icZ,, equal the average value obtained if these were solved for one at a time using this aggregate equation. This yields the solution

C Bis=

5 qrsmW

~1, kc1

ViEIs,

s= l,...,n.

k&k%sk)

1 r, 1 rs s

The initial feasible solution is then obtained by projecting the above solution onto the feasible region given by eqn (lb,c), by using the simple algorithm given in Bitran and Hax (1976), and Sherali and Shetty (1980). Remark 2. We also attempted an initial projected Newton step to enhance the starting solution obtained above. Restricting any split parameter that is less than 0.01 to zero, we computed the reduced gradient over the equality subspace dellned by eqn (1 b). Setting this reduced gradient equal to zero and solving the resulting linear system of equations, we

Parameter optimization methods for O-D trip-tables

tiak&al:

Setf x 1 ad compute aninitial Slution

0’.

147

lxt do i 0.

Determine the &ecrion d mtion d’ = g~j + ytd I-’ using Quati~n (6)

1 Compute the optimal step length Ir using Equarim

(7) and

I

Reset by letting d’- ’ I 0 if the l&going step length k equals )cmrr in Equation (7d). Fig. 1. The main loop of the proposed algorithm.

derived a revised starting solution. If this solution violated eqn (lc), we examined the maximum step-length that could be taken from the initial starting solution to this revised solution while maintaining feasibility, and used the corresponding solution to commence further iterations of the algorithm. Some test results using this modified strategy are presented in Section 4.0 3.1.2. Computation of the gradient g’ at se, for iteration .Ql. Suppose that at any iteration L> 1, we have a current feasible solution @ (see Fig. 1). We compute the gradient Of (j?e)&=(&t) off at p = Be as follows [see eqn (Id)]:

W where ~(k)=~&&&-y#) id,

Vs= I,...,n,and

k= l,...,K.

(4b)

148

Hanif D. Sherali et al.

3.1.3. Computation of the projected anti-gradient g&,j. Recall that in unconstrained convex optimization, a solution is optimal if and only tf the gradient is zero, failing which the anti-gradient provides a direction of descent. Since we have a constrained problem, we need to compute the projection $,. of the anti-gradient -g” onto the set of feasible directions D, and then Be is optima P,C 1 and only if $rroj = 0; else, $rroj yields an improving, feasible direction. (In Section 3.1.3, we further deflect this direction in the latter case.) The projected anti-gradient $rroj is given by the solution to the problem minimize

{11d - (-ge) 11: d E D}

@a)

where, denoting f* = {j E Ji : & = O}Vi = 1, . . . , m, the set of feasible directions D is given by d:xdg=O

Vi=l,...,m

dg?O

and

Vjc$,,

(W

jdi

The computation of $rroj is easily accomplished by the following procedure due to Bitran and Hax (1976) and Sherali and Shetty (1980). Note from eqns (5) that this projection problem is separable in i= 1,. . .,m. Hence, for each entrance i= 1; . .,m, the (g&,j) ,, for jEJi of the vector $rrrojare computed as follows. (The complexity B of this scheme for determining $rroj is 0(mn2).) Initialization: Set the iteration counter t = 0. Define y” = (-dj, jEJi), and MO= Ji. Step 1: Given y’ and M,, compute

components

PI =

-

[

1

CI$ / IMA W,

and

V’ =

yr+ ptur

where ut is a vector of liWrlones. Hence, cjcJi fj = 0. If pj 2 OVjE t,, II M,, then terminate with the projected direction components given by, (5c) Else, go to Step 2. Step 2: Set Mt+r = Mt - j E Jn : fj < O}, y’+’ = { vf : j E Mt+i}, increment t by 1 { and return to Step 1. Having thus determined Groj, if 11$proj 11for some small tolerance so (we used so = 10m4), then /I” is a (near) optimal solution to LSQ, and so the algorithm can be terminated. 3.1.4. Determination of the direction of motion 6e. Continuing to work in the projected subspace defined via the set of active constraints as in eqn (5b), we apply the conjugate gradient strategy to compute the direction of motion by deflecting the projected antigradient direction using the previous direction of motion dl-’ as follows: d =

&j

+

qed-*,

#a)

where $e is a deflection parameter computed as in Fletcher and Reeves (1964) by

II&j II2 ee=llg’,;d,*

(6’9

Note that de-’ is taken as zero if either e = 1, or if the algorithm is reset as noted next, as well as noted in Section 3.1.6 below. If either 11de 115 EOor if d is not sufficiently a descent direction as assessed by the condition de&-e0 (or alternatively, if de +G,,,j 5 0.8 llg’proj112), we reset the algorithm by taking de = $PrOjitself. 3.1.5. Determination of optimal step length le. Having determined a direction of motion de, we next find the optimal step length Ai to take along this direction, exploiting

149

Parameter optimization methods for WD triptables

the special structure of our problem. Note that the maximum possible step length imX that can be taken along dl from the point /I’ and yet maintain feasibility is given by )imaX= minimum~,~,...,m min2 i

Furthermore, equation

the unconstrained

F.

jeJj-4

: 4j .c0 .

(74

1

optimal step length can be computed

by solving the

Now, defining

(7b) and recalling the definition of q(k) from eqn (4b), the foregoing equation becomes

Solving for 1, we get,

(7c)

where the denominator unconstrained

is non-zero

since de is a descent direction and LSQ has an

minimum. Also, note that d2f/d)c2 = C,“,

&

[$(k)]

*> 0. Hence, an

unconstrained global minimum to the function f along the direction dl occurs at 1= 1*. The constrained optimal step length Aceis therefore given by )CC= minimum{k*, Amax).

(7d)

3.1.6. New iterate, termination check, and reset criteria. By using eqns (6) and (7), the revised iterate @+I is then obtained as #I”’ l =Be+&de, and the iteration counter L is incremented by one. To safeguard against the event of a slow asymptotic convergence, in addition to the theoretical termination criterion based on g&,j as described in Section 3.1.3, we adopt a practical criterion by terminating the algorithm in case the objective function value fails to improve by at least 1% over each of four consecutive iterations. Furthermore, the procedure is reset by putting 8-l in eqn (6a) if a new subspace of active constraints is encountered as detected by the foregoing step length being equal to k max< )L*. The algorithm then loops back and repeats the main step, this process continuing until termination occurs. (In real-time implementation, alternative termination criteria based on elapsed solution time or a maximum number of iterations could also be employed. This is also computationally tested in Section 4.) 3.2. Algorithm for problem L1 The L1 estimation problem (3) can be equivalently formulated as the following linear programming problem LP, by denoting each term within the absolute value symbol 1.1as

150

Hanif D. Sherali et al.

the difference of two nonnegative variables S$ - Sk;.,so that the absolute value of this term itself equals the sum of these two nonnegative variables. LP: Minimize K

n

k=l

j=l

cc

@a)

subject to ~q$~ij-(~&-8~)=yj(k)

Vk= l,...,K,

j=

1,

, n,

@b)

kl,

w

l,...,m,

1 Vi=

&=

c

icX

/Iv>0

Vjc Ji,

?i&,S~~O

(84

i= l,*..,m,

Vk=l,...,K,

j=l,...,n.

The dual to LP can be written as follows: DP: Maximize

k=l

j=l

i=l

j=l

subject to K c

q+Wji)

+ WY) + 4ijWy’

+

Wf’

=

0

vj

E Ji,

i=

1, . . . , m

W

k=l

-l
-

Vjzl,...,

n,

unrestricted Vj= 1,. wt2) I unrestricted Vi = 1,. . . ,m ’ wi3) J w!? > 0 Y -

Vj E J.I,

k=l

....K

(9c)

. .,n,

i= l,...,m.

(94

Note that the dual DP has typically far fewer structural constraints in eqn (9b) than the number of structural constraints present in eqns (8b-d) for LP. The bounds on the variables in DP can be readily handled by either the simplex method or by the interior point affine scaling algorithm (see the algorithmic description and test results presented in Sherali et al., 1988). Moreover, modifications in the number of periods K are readily handled since they result in a variation in the columns of DP. We hence consider the dual problem DP for solving the Li estimation problem. In our experiments, we used the affine scaling interior point algorithm to solve the problem DP and also solved the same problem using the commercial software MINOS linked to the modelling language GAMS (see Brooke et al., 1992). GAMS-MINOS performed extremely well for this particular problem, given its moderate size, requiring an order of magnitude less computational time to solve the problem than did the interior point algorithm. We hence advocate this approach for solving LP via DP. Furthermore, in order to compare the least-squares and the LI estimation solutions, we ran problem LP (problem Ll) with constraints eqn (8d) [constraints eqn (3c)] relaxed. Let us refer to this relaxed linear program as LPr, and to its dual as DPr. Test results on solving both DP and DPr are presented in the following section.

Parameter optimization methods for O-D trip-tables 4. COMPUTATIONAL

151

RESULTS

In this section, we present computational results on several experiments that were conducted to test the proposed models and algorithms. Unless otherwise noted, all runs were made on an IBM 3090-300 E/3 computer. These experiments and their results are presented below. 4.1. Evaluation of the proposed algorithm for problem LSQ and comparison with MINOS For the purpose of this study, we generated three categories of realistic, simulated test problems according to the following scenarios, for various entrances (m), exits (n), and time periods (K). Type I: Type II: Type III:

Travel time between each entrance-exit pair is zero (case of an intersection). Travel time between each entrance-exit pair is non-zero and all exits are accessible from each entrance (case of a small freeway section). Travel time between each entrance-exit pair is non-zero and not all exits are necessarily accessible from each entrance (case of unidirectional flows or closed ramps on the freeway or a ‘no left turn’ scenario for an intersection).

The data were generated using SIMUL, a MATLAB-based simulation program written by Van der Zijpp (1992). In this program, for each entrance i, a total volume parameter I’i (over the horizon) is generated uniformly on [2000,6000], and this is divided by the number of periods K. The volume entering each entrance i in each time period k is then generated based on a Poisson distribution having a mean of Vi/K. Next, to determine exit volumes, an a priori split matrix of parameters &~[0,1] is first randomly generated on a uniform distribution, and the rows are scaled so that cj & = 1Vi. In addition, for the test cases 13-18 of Type III, some 1, 1, 1, 2, 2, and 1 entrance-exit pairs, respectively, were randomly prohibited. Now, for each entrance i, using & as the probability that an entering vehicle will take exit j, for each time period, the entering vehicles are each assigned an exit over the horizon. The travel time tii between each entrance i and exit i is taken as zero for problems of Type I, and are computed using unit link travel times for problems of Types II and III, based on some corridor network structure that is assumed within SIMUL. The program SIMUL also permits these link travel times to vary depending on the departure time, based on the previous period’s travel time and some trafIlc dispersion factors. For each entering vehicle, at each departure period, knowing the intended exit and the estimated travel time, we hence compute the exit volume data over the problem horizon. Table 1 presents the test results obtained for the projected conjugate gradient algorithm of Section 3.1, versus the commercial package MINOS 5.3 used in conjunction with GAMS. (MINOS uses a second-order reduced gradient method with BFGS quasi-Newton Hessian approximations.) An examination of these results reveals that with the termination criterion employed, the proposed algorithm produces quite accurate solutions using far fewer iierations than does MINOS. (The ‘mean absolute deviation’ records the value Ci cj 1#?u- &I / mn, where /I is the solution produced by the proposed algorithm and /I* is the solution produced by MINOS.) The execution time is also seen to be very competitive in comparison with MINOS, and is always well under one second. This is considerably less than a typical sampling interval of about a minute (Cremer and Keller, 1987), hence suggesting that the procedure is easily real-time implementable. As a point of interest, we also computed the mean absolute deviation between the optimal solution p* to problem LSQ and the randomly generated split parameters p used within SIMUL to generate the problem data. These mean absolute deviation values,

152

Hanif D. Sherali et

al.

Table 1. Comparison of MINOS and the proposed algorithm for problem LSQ Type

I

II

III

Problem no.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 I5 16 17 18

(m,n,K)

(3,3,60) (4,4,40) (4,4,60) (5,5,20) (5,5,40) (5,560) (3,360) (4,440) (4,4,60) (5,5,20) (5,5,40) (5,560) (3,360) (4,440) (4,460) (595720) (5,540) (5,560)

given by CjCiIS~-S;I/

MINOS

Proposed algorithm

Obj. value

No. iterations

CPU time (s)

5010.767 4408.419 5645.724 2466.596 4476.194 9136.497 3022.268 4529.627 3087.029 1975.759 5748.598 9487.264 5402.935 2116.822 3711.435 560.023 2185.037 5319.210

:;

0.150 0.140 0.150 0.150 0.160 0.170 0.140 0.140 0.150 0.140 0.150 0.170 0.140 0.140 0.150 0.140 0.150 0.160

30 48 49 43 :: 31 43 48 44 :: 5 41 41 41

Cpu time Obj.value No. iterations (s) Qpt.vaIue 3 6 5 6 6 5 8 8 6 8 II 8 5 5 8 7 10 12

0.049 174 0.065649 0.086566 0.057188 0.098971 0.130071 0.096298 0.079072 0.131744 0.069080 0.140545 0.167801 0.059264 0.089197 0.180895 0.127238 0.277812 0.396871

Mean absolute deviation

1.OOOOOO 0.000297 1.OOOOOO0.000225 1.OOOOOO0.000211 1.007539 0.024331 1.000017 0.001128 1.000011 0.001108 1.007440 0.040649 I .000065 0.002948 1.027884 0.045008 1.003034 0.020497 1.000093 0.002376 1.000319 0.005124 1.000000 0.000244 1.OOOOOO0.000226 1.000000 0.000212 1.000000 0.000212 1.006527 0.023965 1.006514 0.021130

mn, turned out to be as follows, respectively, for the 18 test

problems of Table 1: 0.074667, 0.024556, 0.026333,

0.042375, 0.032437, 0.070250,

0.069187, 0.064437, 0.075125,

0.082320, 0.062 160, 0.073880,

0.089240, 0.054400, 0.077800,

0.098480, 0.034840, 0.052840.

Next, we tested the proposed algorithm by restricting the maximum number of iterations to 7, and we also tested a second strategy in which the enhanced starting solution of Remark 2 was used along with a maximum iteration limit of 5. Table 2 presents the results obtained. With a maximum iteration limit of 7, the results for problems of Type I were unaffected, while for Types II and III the final solution value remained within about 2% of optimality, except for problem 18 where the objective value obtained was 5.47% greater than the optimal value. Examining the results obtained for the second experiment reported in Table 2, we notice that the initial solution value produced by the method of Section 3.1.1 is on average five times greater than the optimal evidenced by Table 1 and the results of the first experiment in Table 2; this solution is quickly improved by the proposed algorithm. On the other hand, the enhanced starting solution step suggested in Remark 2 improves this solution to within l-5% of optimality. This enables the algorithm to determine solutions within a 5 iteration limit that are comparable to those produced within 7 iterations of the original algorithm, with the effort also being comparable because of the added effort required by the enhanced starting procedure. Hence, while both procedures are competitive, we recommend the simpler-to-implement initialization method that is described in Section 3.1.1. As an incidental remark, we also attempted applying the projected Newton step of Remark 2 at other than just the initial iteration, i.e. whenever the proposed algorithm yielded less than a 1% improvement in objective value. This improved the accuracy of the solutions produced but required somewhat more computational effort, and hence, in the interest of efficiency and simplicity, this was not pursued further. 4.2. Evaluation of the L1 estimation models In Section 3.2, we proposed two linear programming models LP and LPr for generating least absolute deviation estimates for the split parameters. Table 3 presents results obtained by using GAMS-MINOS to solve the linear programming problems for the

I 2 3 4 5 6 7 8 9 10 II 12 13 14 15 16 17 18

Problem

no.

3 6 5 6 6 5 7 7 6 7 7 7 5 5 7 7 7 7

No. iterations

0.049174 0.065649 0.086566 0.057188 0.098971 0.130071 0.067588 0.069306 0.131744 0.061514 0.103127 0.146962 0.059264 0.089197 0.095678 0.059895 0.102343 0.143790

Cpu time (s)

l.OOOOOO l.OOOoOO l.OOOOOO 1.007539 1.000017 1.000011 1.010672 1.001067 1.027884 1.012897 1.015325 1.003811 1.OOOOOO 1.OOOOOO l.OOOOOl I.013364 1.020614 1.054780

Opt.value

= 7

o.ooo297 0.000225 0.000211 0.02433 1 0.001128 0.001108 0.050368 0.010174 0.045008 0.030103 0.027259 0.016283 o.OOG244 0.000226 0.000390 0.024033 0.040403 0.049714

Mean absolute deviation

no. iterations

Obj.value

Maximum

iteration

5 1 1 1 1 1 1 5 5 5

1

1 5 5 5 5

1

1

No. iterations

Table 2. Results using a maximum

0.047767 0.054314 0.072181 0.074907 0.119648 0.166100 0.068897 0.056251 0.105519 0.059128 0.089665 0.121229 0.048298 0.054375 0.073705 0.076238 0.121611 0.165468

Cpu time (s)

1.033182 l.OOOOOO 1.OOOOOO l.OOOOOO l.OOOOOO 1.OOOOOO l.OOOOOO I .029327 1.059539 1.076922

I .oooooO

1.018141

1.000727

l.OoOOOO l.OOOOOO l.OOOOOO 1It09647 1.001007

Optvalue

Obj.value

0.000277 0.000225 0.000211 0.027071 0.009866 0.008926 0.060179 0.000260 0.050516 0.000245 0.000224 0.000256 0.000244 0.000226 0.000212 0.027789 0.068797 0.051825

1.43 9.53 16.61 4.15 3.84 4.05 11.14 2.95 7.73 3.81 2.11 3.97 1.29 5.59 1.26 4.59 3.12 2.67

Opt. value

= 5 initial value

8 f

1.009 1.009 1.998 1.009 1.046 1.031 1.027

B

5 R v)

B 1.010 1.009

30 3 2

[

4r.

$ p 3 B 3

1.077 1.OlO

1.041 1.009 1.038 1.041 1.111

1.009 1.009

Opt value

Enhanced

no. iterations

Initial value

and maximum

Mean absolute deviation

solution

solution

starting

starting

Enhanced

limit and an enhanced

Hanif D. Sherah et al.

154

Table 3. Results for the Li estimation problems LP and LPr of Section 3.2 Problem no.

Model LP No. iterations

1

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

Model LPr No.

Cpu time Value E LSQ Mean (s) Opt. LSQ value absolute deviation

202 204 351 144 316 514 163 179 299 183 270 521 142 204 309 154 304 446

0.250 0.240 0.290 0.200 0.270 0.340 0.250 0.240 0.290 0.200 0.260 0.340 0.250 0.240 0.290 0.190 0.260 0.340

1.015858 1.048549 1.014790 1.064669 1.016569

1.044566 1.005233 1.034438 1.026456 1.098874 1.031260

1.022869 1.010161 1.021882 1.108323 1.127897 1.020083 1.018566

iterations

0.051111 0.071938 0.037937 0.087440 0.032040 0.070880 0.019556 0.047313 0.061437 0.065320 0.031480 0.036600 0.027 111 0.038312 0.059500 0.067400 0.030120 0.042240

262 260 355 206 353 532 215 245 398 219 325 471 226 208 322 162 271 562

Cpu time Value E LSQ Mean (s) Opt. LSQ value absolute deviation 1.019536 1.065351 1.017447 1.079950 1.022865 1Xi49856

0.250 0.240 0.290 0.200 0.270 0.340 0.260 0.240 0.290 0.200 0.260 0.340 0.240 0.240 0.300 0.190 0.270 0.340

1.025935 1.048432 1.042057 1.062851 1.026676 1.033538 1.005098 1.054438 1.033594 1.134172

1.020834 1.025095

0.056889 0.073187 0.040034 0.082840 0.028960 0.066600 0.041111 0.035625 0.032000 0.055120 0.025209 0.34000 0.019889 0.051563 0.045188 0.038200 0.026120 0.042640

same set of 18 test problems. The table also presents an evaluation of the linear programming solutions when used as estimates for solutions to the corresponding problems LSQ. Specifically, the column ‘value in LSQ/opt. LSQ value’ gives the ratio of the objective function value of problem LSQ corresponding to the solution produced by the linear programming model, to its (LSQ’s) optimal objective value. The mean absolute deviations in these two solutions for each problem are also reported. The results indicate that, on average, LP produces solutions that are within about 4% of optimality for problem LSQ. The mean absolute deviation in the solutions is 0.045 at an average. The cpu time required by MINOS to solve LP is somewhat greater than that required to solve problem LSQ, although still under 0.35 s, thereby easily admitting real-time implementability. (Note that a more specialized linear programming package such as CPLEX might further enhance the solution efficiency. We used MINOS for the sake of uniformity in comparisons.) Similar results are obtained for the relaxed model LPr, although the more restricted model LP produces somewhat more accurate solutions. Hence, if a standard linear (but not nonlinear) programming package is available, the use of an Li estimation model is a viable alternative. (Also, see the following section.) 4.3. Tests using outliers To test the relative effect of outliers on the least-squares model versus the L1 estimation model, we ran the following experiment. For some five problems from Table 1, as indicated in Table 4, we perturbed the data using two levels of perturbations to simulate outliers in data. In ‘perturbation l’, some five exit volumes were randomly selected and a value of 10 was added to each of these volumes. In ‘perturbation 2’, a value of 20 was added to each of the five randomly selected exit volumes. Table 4 gives the effect of this Table 4. Mean absolute deviations under outliers/perturbations Problem type I II II II III

Problem no. 1 3 9 11 13

(m,n,K)

(3,360) (4,460) (4,460) (5,540) (3,360)

Model LPr

Model LSQ Perturbation 1

Perturbation 2

Perturbation 1

Perturbation 2

0.022 0.010 0.023 0.011 0.005

0.044 0.019 0.045 0.022 0.009

0.004 0.000 0.000 0.000 0.000

0.004 0.000 0.000 0.000 0.009

15s

Parameter optimization methods for (TD trip-tables

perturbation on the split parameter values produced for problem LSQ versus those produced for the similarly constrained Li estimation model LPr. For each case, we present the mean absolute deviation (to three significant digits) between the solution produced for the perturbed problem, and that for the corresponding unperturbed problem. Hence, this records the effect of the perturbation on the optimal solution for each problem, and for each model. As expected, the effect of outliers (perturbations) in data is far more pronounced on the least-squares model than it is on the L, estimation model. 4.4. Tests using synthetic data and comparisons with a Kalman jilter approach Cremer and Keller (1987) proposed and tested a Kalman filter approach for the O-D split parameter estimation problem that is applicable to the case of an intersection (‘Type I’ problems). This approach is shown to be very competitive in comparison with other statistical estimation techniques (correlation method and recursive technique), and is quite popular in practice. To test this against the proposed parameter optimization approach, as well as to assess the relative accuracy of solutions produced, we generated five problems of Type I that are ‘synthetic’ in nature, in that the exit volumes were determined exactly using the generated a priori split parameter matrix, and then rounded to integer values. The proposed algorithm as well as the Kalman filter approach were run on these five problems. Table 5 gives the results obtained. The proposed algorithm produced the exact solution with an average effort of 0.093 cpu s. The Kalman filter approach, however, produced solutions that were on average 22.5% greater in objective function value than the optimum (where, as in Table 5, these optimum values range from 227.80 to 778.90) having a mean absolute deviation from the optimal split parameter solution of about 0.0156 on average. The computational effort required was only slightly less, being 0.074 s CPU.Hence, unlike the experience reported in Cremer and Keller (1987) the proposed algorithm demonstrates that it is possible to design a simple, accurate, yet efficient parameter optimization approach that is competitive with statistical estimation procedures. 4.5. Experiment

using a rolling horizon framework

To further test the proposed algorithm, two rolling horizon tests were conducted using a problem of Type I and size (m,n) = (4,4). The first case considered 10 initial periods, and applied the proposed algorithm to determine the associated split parameter matrix. Thereafter, for the next 10 periods, new data were introduced for one additional period at a time, and that for the oldest existing period was dropped, and the split parameter estimates were accordingly updated, using the previous run’s solution as an advanced starting solution for the current period run. Each period’s problem was also solved using GAMSMINOS for testing the accuracy of the solutions obtained. The test for the second case is similar, except that 20 initial periods were used to obtain the first split parameter estimates. These runs were made on a Power Macintosh 6100/66 machine using Symantec C + + 8.0 to demonstrate some implementation results on a personal computer. Table 6 gives the results obtained. Note that for both cases and in each period, the proposed algorithm found the exact solution taking 1-4 iterations to perform an update of the split parameter estimates, including the first update. On average, each step Table 5. Test runs using synthetic data and comparison with a Kalman filter approach

WbK)

Proposed algorithm

Opt.

Kalman filter approach

value

(3,360) (4,440) (4,4,&I)

227.80 343.93 446.32

No. iterations

Cpu time (s)

3 6 8

0.049781 0.065569

0.108571 0.106550 0.136947

Obj. value Opt. value

Mean absolute deviation

O.OOOOOO O.OOOOOO 1.OOOOOO O.OOOOOO 1.OOOOOO 0.000000 1.OOOOOO O.OOObO9 1.OOOOOO 1.oooooO

Cpu time (s)

Obj. value Opt. value

Mean absolute deviation

0.043000 0.047100 0.070100 0.083900 0.124200

1.241884 1.421902 I .212827 1.069472 I.121544

0.010146 0.021645 0.010923 0.022680 0.015907

156

Hanif D. Sherali er al. Table 6. Rolling

Opt. obj value (MINOS)

Proposed No. iterations

Obj.value Opt.value

Case (i): 10 1.0375 1.1177 0.8766 1.1260 1.1361 1.1357 1.1914 1.5069 1.4181 1.4684

horizon

tests

algorithm

Kalman

Mean absolute deviation

*CPU time

Obj.value

(s)

Opt.value

Mean absolute deviation

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.050000 0.033333 0.033333 0.016667 0.016667 0.033333 0.033333 0.016667 0.033333 0.016667

48.047 8.308 16.276 21.922 26.749 36.483 60.962 27.366 14.664 41.936

0.0253 0.0224 0.0260 0.0249 0.0255 0.0276 0.0259 0.0221 0.0219 0.0188

initial periods 4 1 4 2 1 3 1 3 1 1

1.000 1.000 1.000 1.oOO 1.000 1.000 1.000 1.000 1.000 1BOO

*Total cpu time = 0.2833 s Case (i) 2.6320 2.5601 2.6347 2.6820 2.5650 2.7740 2.8221 3.0133 3.0931 3.1266

filter approach

*Total cpu time = 0.1833 s

: 20 initial periodv 4 1 3 3 2 3 1 3

1 1

1.000 1.000 1.ooo 1.000 1.000 1.000 1.OOo 1BOO 1.000 1.000

0.0000

o.ocmo 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0090

0.050000 0.033333 0.050000 0.033333 0.033333 0.033333 0.033333 0.033333 0.033333 0.033333

0.0176 0.0145 0.0142 0.0151 0.0155 0.0166 0.0158 0.0140 0.0153 0.0134

*Total cpu time = 0.2000 s

*Total cpu time = 0.3667 s *Note: Cpu time is based on runs made on a Power Macintosh

44.365 8.411 18.434 17.365 9.365 19.616 34.354 17.506 18.895 11.831

6100/66 using Symantec

C+ + 8.0.

consumed about 0.02833 s cpu time for the first case, and 0.03667 s cpu time for the second case, a very nominal effort. In comparison, although the Kalman filter estimates were determined more rapidly (at an average cpu time of 0.01833 s for the first case and 0.02 s for the second case), they were far less accurate. At an average, for the first case where 10 initial periods were used, the objective function value of the solutions obtained was about 30 times the optimal value (where from Table 6, the optimum value was 1.2014 on average), with a mean absolute deviation of the split parameters from the optimal solution being 0.024. For the case of 20 initial periods, because of the added smoothing effect, these figures were, respectively, 20 (the optimum value being 2.790 at an average), and 0.0152. Note also that, unlike the proposed approach, the Kalman filter procedure is unable to drop older data except by way of diminishing its effect as the periods progress, unless if it is re-initialized from scratch. Being wasteful, the latter alternative is not implemented in practice. In fact, we can observe from Table 6 that the accuracy of the solution found for the first initial set of periods was no better than that of the remaining updated solutions that have accumulated the information for one additional period at a time without dropping older data. 5. CONCLUSIONS

AND COMPARISONS

In this paper, two parameter optimization models and algorithms for estimating split parameter matrices have been presented. The first is a constrained least-squares approach, and the second is an L1 estimation approach, posed as a linear programming problem. A projected conjugate gradient method is used for solving the constrained least-squares problem and the solution to the linear programming problem is obtained using a commercial simplex-based algorithmic package. The results exhibit that properly designed parameter optimization methods can indeed be operated in real-time implementations, providing solutions that are more accurate than those obtained via statistical techniques,

Parameter optimization methods for O-D trip-tables

157

but with a very nominal effort that is competitive with that required for such methods. Moreover, we have demonstrated that if access to a linear programming package is available, then an Li estimation problem can be solved quite efficiently to compute O-D trip-tables having a reasonable proximity to least-squares solutions. In fact, this might be a recommended alternative in cases of inaccuracies in data. Furthermore, it is important to mention that our formulation also incorporates travel time delays that have been largely ignored in the other research efforts in this area (see Remark 1). This is an important factor which makes our methods applicable to (limited) linear freeway sections. Cremer and Keller (1987) have stated that their formulations can incorporate travel time delays at the cost of increased computational effort, although no details or computational results have been presented on this account. Our research has demonstrated that travel time delays can be incorporated with a modest increase in computational effort, still well within the computational limit imposed by real-time implementation considerations. Acknorvledgemen~s-This work has been supported by General Motors Inc., the Federal Highway Administration (RCE Grant), and the US and Virginia Departments of Transportation. The authors are grateful to three anonymous referees whose detailed comments served to significantly improve the contents and presentation in this paper. The authors also acknowledge Mr Taehyung Park for his assistance with the computational experiments.

REFERENCES Bazaraa, M. S., Sherali, H. D. and Shetty, C. M. (1993) Nonlinear Programming: Theory and Algorithms, second edition. John Wiley and Sons, New York. Bell, M. G. H. (1991) The real-time estimation of origindestination flows in the presence of platoon dispersion. Transpn Res. 2SB(2/3), 115125. Bitran, G. and Hax, A. (1976) On the solution of convex knapsack problems with bounded variables. Proceedings of IXth International Symposium on Mathematical Programming, Budapest, pp. 357-367. Box, M. J. (1965) A new method of constrained optimization and a comparison with other methods. Comput. I. 8,42-52.

Brooke, A., Kendrick, D. and Meeraus, A. (1992) GAMS: A Users Guide, Release 2.25. The Scientific Press, San Francisco, CA. Cremer, M. and Keller, H. (1981) Dynamic identification of flows from traffic counts at complex intersections. Eighth Internatiorud Symposium on Transportation and Trajlc Theory, pp. 121-142. Cremer, M. and Keller, H. (1983) A systems dynamic approach to the estimation of entry and exit O-D flows. Ninth International Symposium on Transportation and Trafic Theory, pp. 431450. Cremer, M. and Keller, H. (1987) A new class of dynamic methods for the identification of origin-destination flows. Trwupn Res. 25B(2/3), 117-132. Fletcher, R. and Reeves, C. (1964) Function minimization by conjugate gradients. Comput. I. 7, 149-154. Sherali, H. D. and Shetty, C. M. (1980) On the generation of deep disjunctive cutting planes. Naval Res. Logistics Quart. 27(3), 453-475. Sherali, H. D., Skarpness, B. 0. and Kim, B. (1988) An assumption-free convergence analysis for a perturbation of the scaling algorithm for linear programs, with application to the Ll estimation problem. Naval Res. Lagistics Quart. 35(S), 473492. Nihan, N. L. and Davis, G. A. (1987) Recursive estimation of C&D matrices from inuut/outnut _. _ counts. Transtm Res. 21B(2), 149-163. Nihan, N. L. and Davis, G. A. (1989) Application of prediction-error minimization and maximum likelihood to estimate intersection O-D matrices from trafBc counts. Transpn Sci. 23(2), 77-W. Van der Zijpp (1992) Personal communication, University of Delft, The Netherlands.