As automated vehicles gradually become available to travelers, many cities will experience a mixed trac ow consisting of both legacy and automated vehicles. Al- though the overall market penetration of automated vehicles may be known, the pro- portion of automated vehicles may vary in space and time due to spatial and temporal variations in automated vehicle demand. Since automated vehicles are expected to behave di erently than legacy vehicles, this results in a ow-density relationship that varies in both time and space with the local proportion of automated vehicles. We model this scenario using a multiclass kinematic wave theory. Assuming a triangular ow-density relationship (with shape parameters varying with the automated vehicle proportion), we develop a multiclass Newell's method for nding exact solutions to the multiclass kinematic wave theory. The solution method takes the form of a lin- ear program with postprocessing. We then extend this method to a multiclass link transmission model. We develop a faster solution method for the receiving ow which consists of iteratively solving a system of linear equations. Numerical results from dy- namic trac assignment on the downtown Austin city network demonstrate the com- putational tractability of this method and explore the e ects of automated vehicles on trac congestion. 1 Introduction Many previous studies have predicted that automated vehicles will behave signi cantly dif- ferently in trac than legacy (human-driven) vehicles. Of course, for platooning (cooperative adaptive cruise control), automated vehicles (AVs) have a signi cantly di erent ow-density relationship (Shladover et al., 2012; Melson et al., 2018; Ye and Yamamoto, 2018). Even when automated and legacy vehicles are mixed together, previous studies have suggested that di ering reaction times and car-following behaviors could lead to greater capacity (Levin and 1 Boyles, 2016; Ghiasi et al., 2017) and stability (Schakel et al., 2010; Talebpour and Mah- massani, 2016) of trac ow. Since AVs will likely comprise only a fraction of all vehicles in the network for many years, predicting the future impacts of a mixed trac ow is a relevant problem for practitioners and policymakers. This necessitates a model in which the trac ow varies based on the dynamic proportion of AVs. Although the overall market penetration of AVs in a city may be known, spatial and temporal variations in demand and turning movements may result in space- and time-varying AV proportions that di er signi cantly from the network-wide market penetration. For instance, AVs may be substantially more expensive than legacy vehicles. Higher-income residences in the network may generate a signi cantly higher proportion of AV trips than lower-income residences, resulting in spatial variations across the network. Perhaps the most obvious modeling approach is to construct microsimulation models in which individual vehicles use di erent car-following models based on their driver type (legacy or automated) (e.g. Schakel et al., 2010; Zhu and Zhang, 2018). Unfortunately, the detailed dynamics within microsimulation require high computation times for large city networks. Solving dynamic trac assignment (Chiu et al., 2011) to predict the e ects of route choice on trac congestion further increases the computation time. Solving trac assignment is important because Braess (1968) and Daganzo (1998) showed that improvements to the network capacity could result in greater congestion, which Melson et al. (2018) demonstrated for AVs speci cally. Although they are not as detailed as microsimulation, dynamic network loading models based on the kinematic wave theory of trac ow (Lighthill and Whitham, 1955; Richards, 1956) are used in trac ow modeling for their ability to model time-dependent trac ow with manageable computation times. For instance, the dynamic trac assignment model VISTA (Ziliaskopoulos and Waller, 2000) uses the cell transmission model (Daganzo, 1994, 1995) to propagate ow. We propose a multiclass kinematic wave theory in which the ow-density relationship varies in space and time with respect to the local proportion of AVs. Levin and Boyles (2016) developed a multiclass cell transmission model to solve the multiclass kinematic wave theory and used it for dynamic network loading. However, the cell transmission model is only an approximate solution to the kinematic wave theory and contains numerical errors, particularly in the propagation of congested states. Newell (1993) observed that the kinematic wave theory is easily solved along certain characteristic curves, which Yperman et al. (2005) converted into the link transmission model (LTM). We develop a multiclass LTM for the multiclass kinematic wave theory, which is made challenging by the fact that characteristic curves may have to pass through space- and time-dependent variations in the ow-density relationship. The kinematic wave theory can be alternatively solved by transforming it to Lagrangian coordinates. For instance, van Wageningen-Kessels et al. (2009, 2010) used Lagrangian co- ordinates to solve a multiclass kinematic wave theory, although their work did not consider a mix of human and automated vehicles changing the ow-density relationship. Lagrangian coordinates are de nitely worth consideration for modeling a mix of human and AV ows. However, our goal is to construct a model suitable for dynamic network loading in a network 2 of links and nodes. To the best of the authors' knowledge, all kinematic wave theory-based models for dynamic network loading on realistic city networks are solved using Eulerian co- ordinates. This is because intersections would be \moving" in the Lagrangian perspective while in Eulerian coordinates they remain stationary. Therefore, we have chosen to address this problem with Eulerian coordinates, and we construct a link ow model which is compat- ible with the large body of existing work on dynamic network loading and dynamic trac assignment based on the kinematic wave theory. The contributions of this paper are as follows. We formulate a multiclass kinematic wave theory and derive some general results that assist with the structure of the solution method. We then assume that the ow-density relationship is triangular, with the congested wave speed varying with the proportion of automated vehicles. This assumption is used to derive a multiclass version of Newell (1993)'s method to nd exact solutions to the multi- class kinematic wave theory. The solution method involves solving a linear program then post-processing the optimal solution with other constraints to nd an exact value for the cumulative count at an arbitrary ( t;x) point. The primary motivation of this paper is to extend the multiclass Newell's method to a multiclass LTM. We present a faster solution algorithm for the multiclass LTM which involves iteratively solving a system of linear equa- tions. Numerical results on a calibrated large city network demonstrate solving a dynamic trac assignment model using this multiclass LTM, and explore the impacts of AV market penetration on city-wide congestion. The remainder of this paper is organized as follows. Section 2 reviews literature rele- vant to multiclass modeling of trac ow, in particular automated vehicles. In Section 3, we formally de ne the multiclass kinematic wave theory. Section 4 derives the multiclass Newell's method, rst presenting a simpli ed example, then presenting the general solution. Section 5 adapts the multiclass Newell's method to the speci c characteristics of LTM. Nu- merical results from dynamic trac assignment are presented in Section 6, and we conclude in Section 7. 2 Background For decades, trac assignment models have modeled travel times as a continuous function of trac volumes for their favorable analytical properties (Beckmann et al., 1956). Such models have recently been extended to mixed legacy and AV ow. Levin and Boyles (2015) modi ed the standard Bureau of Public Roads travel time function to include a capacity that varied with the proportion of AVs. However, static trac assignment models cannot capture how time-varying proportions of automated vehicles a ect trac congestion. Although the over- all proportion of AVs in a network may be known, local variations in the proportion a ect the local capacity and propagation of congested waves. The e ects of time-varying conges- tion (separate from AV ow di erences) have long been recognized and modeled through microsimulation or mesoscopic ow models such as the kinematic wave theory (Lighthill and Whitham, 1955; Richards, 1956). Research on ecient solutions to the kinematic wave the- ory has been used for dynamic network loading, yielding dynamic trac assignment models 3 that can predict time-dependent congestion across city networks (Chiu et al., 2011). The multiclass ow of a mix of legacy and automated vehicles creates additional challenges for these network models. In addition to time-varying densities, the ow-density relationship itself can vary in space and time with the proportion of AVs present. These variations in the ow-density relationship move with the trac ow, and are therefore interdependent on the solution to the kinematic wave theory. Modeling multiclass ow is important because of the demonstrated behavioral di er- ences of AVs. Results from microsimulation (Kesting et al., 2007; Davis, 2007; Shladover et al., 2012; Zhu and Zhang, 2018; Liu et al., 2018) and live demonstrations (Ploeg et al., 2011; Milan es et al., 2013) on homogeneous segments of road have shown increases in capac- ity from cooperative adaptive cruise control (Shladover et al., 2015) that reduces following headways. It is not clear whether these capacity improvements apply to networks as well, but these studies nevertheless demonstrate that multiclass ow may behave di erently from legacy vehicle ow. Connected and cooperative vehicle control has also resulted in improve- ments in string stability (Schakel et al., 2010; Ploeg et al., 2011; Talebpour and Mahmassani, 2016), although these results may not be true of non-cooperative AVs. One major limita- tion demonstrated trac ow improvements is driver comfort (Nowakowski et al., 2010). Improvements in capacity and stability from automated driving largely result from lower re- action times admitting shorter following headways (Kesting and Treiber, 2008; Chen et al., 2017; Makridis et al., 2019). Until vehicle manufacturers and drivers are comfortable with short headways, capacity increases may be limited (Vander Werf et al., 2002; Shladover et al., 2012; Makridis et al., 2019), or capacity may even decrease (Calvert et al., 2017; James et al., 2019). Nevertheless, most studies suggest signi cant (positive or negative) changes in trac ow resulting from AVs. Some of these trac ow di erences are assumed to be limited to AV-exclusive lanes (Chen et al., 2016; Vander Laan and Sadabadi, 2017; Liu and Song, 2019; Melson et al., 2018) where connectivity and automation are minimally interrupted by hu- man driving. Other studies have explored technologies for ecient automated driving within mixed trac ow (Kesting et al., 2007; Gong and Du, 2018) hat vary with the AV market penetration (Talebpour and Mahmassani, 2016; Zhu and Zhang, 2018; Ye et al., 2018). This paper focuses on the latter category. Since the exact characteristics of future automated driving are dicult to determine, the focus of this study is on solving a multiclass kinematic wave theory where the ow-density relationship can be adjusted to model AV di erences in trac ow. Although many previous studies have used car-following models of automated driving to model trac ow di erences, the computation times required for these models often discour- age their use for solving dynamic trac assignment on large city networks. Static multiclass trac assignment models (Levin and Boyles, 2016; Wang et al., 2019) have been studied but fail to capture the time-varying changes in congestion, especially due to the varying ow-density relationship. The kinematic wave theory is often an useful middle ground with dynamic trac ow but more computational eciency than microsimulation. Thus far, multiclass variants of the kinematic wave theory have only been solved by approximation methods like the cell transmission model. This is partially because solving the standard 4 kinematic wave theory is challenging by itself. Although the kinematic wave theory was pro- posed decades ago by Lighthill and Whitham (1955) and Richards (1956), ecient solutions were discovered much later. The cell transmission model (Daganzo, 1994, 1995), a Godunov (1959) approximation, is still widely used for dynamic network loading. Newell (1993)'s method for nding exact solutions with a triangular ow-density relationship was used in the relatively recent development of LTM (Yperman et al., 2005). A Lax-Hopf method for exact solutions to the general kinematic wave theory was developed by Claudel and Bayen (2010a,b). Given the diculty in solving the standard kinematic wave theory, it is perhaps not sur- prising that most previous work on multiclass kinematic wave theory has relied on numerical approximation methods. Levin and Boyles (2016) developed a multiclass cell transmission model which was used by Patel et al. (2016) to analyze city networks. Tiaprasert et al. (2017) developed an extended multiclass cell transmission model with better modeling of rst-in- rst-out and overtaking behaviors. Zhu and Ukkusuri (2018) and Chen et al. (2020) used multiclass cell transmission models to study how variable speed limit control for AVs would a ect surrounding legacy vehicle trac. The multiclass kinematic wave theory has previously been studied outside of the con- text of AVs (Logghe and Immers, 2008; Jin, 2012; van Wageningen-Kessels et al., 2014). The multiple classes could be motivated by heterogeneity in physical vehicle characteris- tics (Van Lint et al., 2008; Liu et al., 2015) or driver behavior (Van Aerde and Rakha, 1995). van Wageningen-Kessels et al. (2009, 2010) used a Lagrangian transformation to solve the problem. To the best of our knowledge, a LTM for eciently solving other multiclass kine- matic wave theory problems has yet to be proposed. Although the motivation and context for this paper is focused on mixed legacy and automated vehicle trac, the ideas may be relevant to the multiclass kinematic wave theory in other contexts. 3 Multiclass kinematic wave theory Consider a length of road from x= 0 tox=L. LetNm(t;x) be the cumulative count of vehicles of class m, i.e the number of vehicles that have passed point xat or before time t. The total cumulative counts are given by N(t;x) =P m2MNm(t;x) whereMis the set of all vehicle classes. Let qm(t;x) =@Nm(t;x) @tandkm(t;x) =@Nm(t;x) @xbe the ow and density of class m, respectively. For completeness we de ne the total ow and density as q(t;x) =P m2Mqm(t;x) andk(t;x) =P m2Mkm(t;x), respectively. Assume that the cumulative counts for each class satisfy a conservation law, i.e. @qm(t;x) @x+@km(t;x) @t= 0 (1) which results in a conservation law for the total cumulative counts also. Let q(t;x) =f(k(t;x)) (2) 5 be the vector-valued ow-density relationship which speci es the ow of each class as a function of the class-proportions of densitykm(t;x) kand total density. For the results in Section 3, f(k) may be a general ow-density relationship with any concave shape. Assume thatf(k) is not separable by class, so that the partial di erential equations cannot be solved separately. For instance, ow may be 0 when the total densityP m2Mkm(t;x) is equal to some maximum density value, meaning the ow of one class depends on the density of other classes. Assume rst-in- rst-out (FIFO) behavior within the link. The problem is to nd Nm(t;x) for any (t;x) point given initial and boundary condi- tions. Initial condition c(0;x) speci es the initial cumulative counts at time t= 0. Up- stream boundary conditions c(t;0) and downstream boundary conditions c(t;L) specify re- strictions on entering and exiting ow, respectively. These conditions are generally relevant. LTM (Yperman et al., 2005) uses only upstream and downstream boundary conditions, and the generalized Lax-Hopf formula (Claudel and Bayen, 2010a,b) is constructed for similar initial and boundary conditions, albeit for a single class. For formality, we observe that c(t ;x)can be de ned everywhere, except with c(t;x) =1where (t;x) is not an initial point (t= 0), upstream point ( x= 0), or downstream point ( x=L). The class-speci c values are necessary to propagate class proportions according to the speci ed ow-density relationship f(k). For instance, Levin and Boyles (2016) proposed a ow-density relation- ship for mixed human and automated-vehicle trac ow that varies with the proportion of automated vehicles. Claudel and Bayen (2010a) observed that the single-class kinematic wave theory can be written as a Hamilton-Jacobi partial di erential equation: @N(t;x) @t  @N(t;x) @x = 0 (3) which leads to an exact formula for solving N(t;x) (for a single class): N(t;x) = inf u;T0fc(tT;x+Tu) +T'(u)g (4) wherec(t;x) de nes the conditions imposed on the solution, and '(u) is the Legendre- Fenchel transform of the ow-density relationship. For a triangular ow-density relationship, this formula simpli es to Newell (1993)'s method (Jin, 2015). Unfortunately, neither the Hamilton-Jacobi theory nor the method of characteristics are designed to propagate class proportions in time and space, in addition to ows and densities. Hence, Levin and Boyles (2016) proposed a cell transmission model to nd an approximate solution. However, due to numerical errors in such approximations, an exact solution is preferable. 3.1 Class proportion regions Assume that the order of vehicles does not change within x2[0;L], i.e. each vehicle can be labeled with a cumulative count number that remains constant within x2[0;L]. This is equivalent to assuming rst-in- rst-out behavior. In reality, this assumption is broken by 6 lane-changing. However, within the kinematic wave theory where vehicle speed is a function of density, this assumption is reasonable. This assumption ensures that class proportions move with vehicles. Shockwaves occur in the standard (single class) kinematic wave theory as the separation of regions of di erent density. Some shockwaves move with vehicles, and others move backwards in space, meaning that vehicles can cross a shockwave boundary. In the multiclass kinematic wave theory, shockwaves separate regions of di erent total density as well as regions of di erent class proportions. Since class proportions move with vehicles, the class proportions can be described with respect to vehicle labels. Consider Idi erent regions of class proportions. Region iis de ned as the vehicles between vehicle bi1andbiwhich are the boundary vehicles. For instance, a region might be de ned as the vehicles that enter the upstream boundary between a speci ed time, i.e.bi1=N(ti1;0) andbi=N(ti;0) whereti1<ti. Within a region, vehicle class proportions are constant. However, note that class proportions are not all-or-nothing within a region. It is fully possible to have k1(t;x)>0 andk2(t;x)>0 with a corresponding ow-density relationship f(k). The key de ning property is that the proportionkm(t;x)P m0km0(t;x) is constant throughout each region. For example, this could be used to model a region of density with 50% legacy vehicles and 50% AVs, uniformly distributed throughout the region. Proposition 1. When trac ow obeys rst-in- rst-out behavior, class proportion regions move with vehicles. Proof. Suppose that at x= 0, the class proportions between b1andb2areNm(t2;0)Nm(t1;0) N(t2;0)N(t1;0) whereN(t1;0) =b1andN(t2;0) =b2. At any points ( t0 1;x0 1) and (t0 2;x0 2) such that x0 1L andx0 2LwithN(t0 1;x0 1) =b1andN(t0 2;x0 2) =b2, Nm(t0 2;x0 2)Nm(t0 1;x0 1) N(t0 2;x0 2)N(t0 1;x0 1)=Nm(t2;0)Nm(t1;0) N(t2;0)N(t1;0)(5) due to FIFO. Each region has an unique ow-density relationship q(t;x) =fi(k(t;x)) that is speci c to the class proportions of region i. When class proportions are constant within regions, the class-dependent ow-density relationship f(k) can be simpli ed to a scalar function q(t;x) = fi(k(t;x)) that describes the ow within that region. Similarly, the initial and boundary conditions can be written as a scalar c(t;x), with the bilabels indicating the boundary vehicles between di erent class proportion regions. Figure 1 shows space-time trajectories of vehicles of two classes as they enter a queue. The speed of the shockwave separating the queue from the uncongested region behind it varies based on the ow-density relationship of the vehicle class. However, the class regions move with the trac. Consequently, the multiclass cumulative counts Nm(t;x) can be simpli ed to tracking a single cumulative count N(t;x) while also tracking the size and location of class proportion regions as function of time. 7 𝑥 𝑡Vehicle class 1 Vehicle class 2𝐿Figure 1: Illustration of space-time diagram for vehicles of two di erent classes, with di erent ow-density relationships. 3.2 Space- and time-varying ow-density relationship The standard kinematic wave theory is typically used with a space-varying ow-density rela- tionship (varying across links due to capacity and other factors). A space- and time-varying ow-density relationship is admissible also. The ow-density relationship fiis included via the Legendre-Fenchel transform ' i(ui) = sup k2Dom(fi)fkui+fi(k)g (6) corresponding to region i(Claudel and Bayen, 2010a). Notice that the boundaries of these space- and time-varying regions depend on vehicle trajectories. This is not limiting because the trajectory of the lower bound of region i(here lower bound refers to the rst vehicle of regioniwhich has the lowest cumulative count label among all vehicles in i) depends on the movements of vehicles ahead of it, i.e. vehicles in region i1. Hence the location of the boundary between regions i1 andican be known before computing the cumulative counts within region i. Figure 2 illustrates how the class regions might a ect the search along characteristics. Passing the boundary between regions 1 and 2, a di erent congested wave speed is used for the characteristic speed. This relationship is formalized in Proposition 2. Proposition 2. Letc(t;x)describe the initial or boundary conditions. The Lax-Hopf formula for a ow-density relationship that varies for Idi erent space-time regions is N(t;x) = inf ui2[vi;wi];Ti2R+( c tIX i=1Ti;x+IX i=1Tiui! +IX i=1Ti' i(ui)) (7) 8 Region 2𝑡1,𝑥1 −𝑤2 𝑣−𝑤1 𝑡,𝑥Region 1𝑥 𝑡𝐿Figure 2: Illustration of searching along characteristics to evaluate N(t;x). such that the point ti+1P j=ITj;x+i+1P j=ITjuj! corresponds to the trajectory of the boundary vehiclebi. Proof. By induction on I.Basis: By Theorem 3.1 of Claudel and Bayen (2010a), the Lax-Hopf formula for a single region is N(t;x) = inf u12[v1;w1];T12R+fc(tT1;x+T1u1) +T1' 1(u1)g (8) = inf u12[v1;w1];T12R+( c t1X i=1Ti;x+1X i=1Tiui! +1X i=1Ti' i(ui)) (9) Inductive step: Given the space-time trajectory of vehicle bI, the Lax-Hopf formula for regionI+ 1 is N(t;x) = inf uI+12[vI+1;wI+1];TI+12R+ c(tTI+1;x+TI+1uI+1) +TI+1' I+1(uI+1) (10) = inf uI+12[vI+1;wI+1];TI+12R+ bI+TI+1' I+1(uI+1) (11) because the lower boundary of region I+ 1 has a cumulative count number of bI. The space- time trajectory of bIis determined by the Lax-Hopf formula for Iregions, i.e. equation (7), which then yields N(t;x) = inf ui2[vi;wi];Ti2R+( c t0IX i=1Ti;x0+IX i=1Tiui! +IX i=1Ti' i(ui) +TI+1' I+1(uI+1)) (12) 9 with (t0;x0) on the space-time trajectory of bI. Equation (12) can be simpli ed to N(t;x) = inf ui2[vi;wi];Ti2R+( c tI+1X i=1Ti;x+I+1X i=1Tiui! +I+1X i=1Ti' i(ui)) (13) forI+ 1 regions. Notice that Proposition 2 does not depend on the shape of the ow-density relationship (other than requiring concavity). As discussed in Claudel and Bayen (2010a), Jensen's inequality simpli es the search over any curve connecting ( t;x) to a boundary point ( t0;x0) such thatc(t0;x0) is de ned. Given an arbitrary ( t;x) point, the Lax-Hopf formula (7) can detect when ( t;x) lies outside of region i. Proposition 3 shows that if ( t;x) is outside of region i, then calculating N(t;x) through equation (7) using the rst iregions will return a value greater than bi, indicating that ( t;x) lies beyond region i. It may not return the correct value of N(t;x), but the relationship N(t;x)>biis sucient to determine that regions i+ 1,i+ 2,i+ 3, . . . need to be included to obtain the correct N(t;x) value. Proposition 3. Suppose that (t;x)lies outside of region i, i.e.N(t;x)> bi. Then the Lax-Hopf formula (7)using regions 1:::i will return a value N0(t;x)>bi. Proof. Since vehicle biis within region i, the trajectory of bican be determined by the Lax- Hopf formula (7) using the rst iregions. Therefore, if N(t;x)>bi, the Lax-Hopf formula (7) must return a (possibly incorrect) value N0(t;x)> bisince the vehicle label at ( t;x) is greater than bi. 4 Multiclass Newell's method For the remainder of this paper, like Newell (1993) we assume that the ow-density relation- shipfi(k) for each region ihas the triangular shape fi(k) = minfvk;wi(kK)g (14) whereKis the jam density, wiis the congested wave speed, and vis the free ow speed. Sections 4{4.2 and 5 apply for any triangular ow-density relationship. Only the numerical results in Sections 4.3 and 6 use a speci c shape of the ow-density relationship from Levin and Boyles (2016). To further simplify the problem, we make assumptions based on the expected nature of our intended context, i.e. a mixture of human and automated vehicles. Since Kdepends on the physical length of vehicles, we assume that Kis independent of the class proportions. We also assume that vis independent of the class proportions because it is based on the speed limit for the road. A single free ow speed value of vsimpli es the search among forward-moving characteristics, but the method could easily be generalized to a region- speci c free ow speed. However, wiis expected to vary based on class proportions. Levin 10 and Boyles (2016) predicted that higher wivalues were likely for automated vehicles, resulting in correspondingly higher capacities. Let kc ibe the critical density of region icorresponding to capacityQi=vkc i. The critical densities and capacities also depend on the class proportions because they vary with the congested wave speed. The Legendre-Fenchel transform of fi(k) is ' i(u) = sup k2Dom(fi)fku+fi(k)g=vkc i+ukc i (15) de ned for wave speeds u2[v;wi]. As observed by Newell (1993) and Himpe et al. (2016), since' i(u) is linear, the in mum in equation (7) occurs at a characteristic wave speed of u=voru=wi. 4.1 Two-region multiclass Newell's method To develop intuition about the general solution method, we rst present the multiclass Newell's method for two regions. The upper bound of region 1 is b1, so region 1 is active for vehicles [1;b1], and region 2 is active for vehicles ( b1;1]. As illustrated in Figure 2 and proven in Proposition 2, the search along backwards-moving characteristics may involve two di erent congested wave speeds w1andw2corresponding to regions 1 and 2 of the ow- density relationship. Since the free ow speed is uniformly v, forward-moving characteristics always move at speed v, which avoids the need to adjust the speed of the forward-moving characteristics when passing region boundaries for this problem. The objective is to nd N(t;x) for an arbitrary point ( t;x)2R+[0;L]. With two class- proportion regions, ( t;x) could fall into either of the two regions. It is not known a priori which region ( t;x) belongs to, as that depends on the trajectory of the boundary vehicle b1. By Proposition 3, the rst step is to calculate N(t;x) using only region 1: N(t;x) = inf u12[v;w1];T12R+fc(tT1;x+T1u1) +T1' 1(u1)g (16) This can be evaluated using Newell's method, tracing the forward- and backward-moving characteristics to the known condition c(t;x). If the returned value of N(t;x) satis es N(t;x)b1, then (t;x) is in region 1, and the returned N(t;x) value is correct. Other- wise,N(t;x) must be calculated using region 2. Because the triangular ow-density relationship is used, there are two possible charac- teristic speeds to check: vandwi(Newell, 1993). Figure 2 illustrates that the problem is to nd where the trajectory of vehicle b1intersects the line of slope w2at point (t1;x1) satisfying x1x=w2(t1t) (17) The point ( t;x) is where we want to nd the cumulative count N(t;x). (t1;x1) is the point on the trajectory of the boundary vehicle b1intersecting a line of slope w2from (t;x). However, the trajectory of b1may not be known exactly. That trajectory is determined as the set of (t;x) points such that N(t;x) =b1| which can be evaluated using the Lax-Hopf formula (7). Figure 3 illustrates this search for ( t1;x1) along the two characteristic lines of 11 Region 1 Region 2𝑡1,𝑥1 −𝑤2 𝑣𝑡,𝑥𝑥 𝑡𝐿(a) Uncongested characteristic Region 2𝑡1,𝑥1 −𝑤2−𝑤1 𝑡,𝑥Region 1𝑥 𝑡𝐿 (b) Congested characteristic Figure 3: Illustration of searching along characteristics to evaluate ( t1;x1). slopevandw1that determine N(t1;x1). These result in two systems of equations that can be evaluated separately. 4.1.1 Uncongested characteristic for (t1;x1) The problem is to nd ( t1;x1) such that N(t1;x1) =b1, (t1;x1) lies on the characteristic of slopew2passing through ( t;x), and that ( t1;x1) lies on a characteristic of slope vsuch thatc(t0;x0) =b1and x0x1=v(t0t1) (18) Equations (17) and (18) together have variables ( t0;x0) and (t1;x1). However, the point (t0;x0) can be simpli ed by splitting it into two cases: either ( t0;x0) is part of an initial condition with t0= 0, or (t0;x0) is part of an upstream condition with x0= 0. SinceN(t;x) is constant along a characteristic of slope v(Newell, 1993), in either case the remaining variable of x0ort0can be determined by evaluating where c(t0;x0) =b1. That results in two equations (17) and (18) with two variables ( t1;x1), which can be solved to obtain solution (tu 1;xu 1) where the superscript u is used to denote the solution from the uncongested characteristic. 4.1.2 Congested characteristic for (t1;x1) When a second congested characteristic of slope w1determines ( t1;x1), the second equation to determine point ( t0;x0) on the initial or boundary condition becomes x0x1=w1(t0t1) (19) SinceN(t;x) increases at the rate of Kxalong backwards-moving characteristics (Newell, 1993),N(t1;x1)N(t0;x0) =K(x0x1). The point ( t1;x1) satis esN(t1;x1) =b1, but 12 sincex1is a variable, the point ( t0;x0) cannot be determined by b1alone. Instead, a third equation can be written: c(t0;x0) +K(x0x1) =b1 (20) The system of equations (17), (19), and (20) has 3 variables | ( t1;x1) and (t0;x0) with t0= 0 (initial condition) or x0=L(downstream boundary condition). Solving this system yields (tc 1;xc 1) where the superscript c is used to denote the solution from the congested characteristic. If multiple solutions exist, the correct solution minimizes b1+K(x0x1) due to the in mum in the Lax-Hopf formula (7). 4.1.3 Combining (tu 1;xu 1)and (tc 1;xc 1) Due to the in mum, the correct solution is the characteristic with the lowest N(t;x) value. Sinceb1is a constant vehicle label, we have two possible values for the point where the trajectory of b1intersects the line with slope w2through (t0;x0): (tu 1;xu 1) and (tc 1;xc 1). Since @N @t0 and@N @x0, and the congested characteristic through ( t0;x0) has a negative slope ofw2,N(t;x) is nondecreasing along that characteristic. Therefore, the correct value of (t1;x1) for the trajectory of vehicle 1 is t1= maxftu 1;tc 1gsince that corresponds to a later arrival time of vehicle b1, or equivalently, a lower value of N(t;x) along the characteristic of slopew2. There are two particular corner cases associated with the value ( t1;x1). First, observe thatx1> x(andt1< t). Otherwise, point ( t;x) is in region 1, and that case was handled already. Second, it is possible that x1> L, i.e. the intersection of the backwards-moving characteristic of slope w2with the trajectory of vehicle b1is outside of the range x2[0;L]. In that case, c(t;x) is the limiting condition, and the standard Lax-Hopf formula for region 2 can be used: N(t;x) = inf u22[v;w2];T22R+fc(tT2;x+T2u2) +T2' 2(u2)g (21) Once the correct value of ( t1;x1) is determined, with x < x 1L, then the maximum value ofN(t;x) based on the congested characteristic is b1+K(xx1) (Newell, 1993) since N(t1;x1) =b1and the cumulative count increases by Kxalong a characteristic moving at speedw2in region 2. The second bound on N(t;x) is based on the forward moving characteristic: nding point ( tu 0;xu 0) such that c(tu 0;xu 0) is de ned (or less than 1) and (tu 0;xu 0) is on a line of slope vpassing through ( t;x), i.e.xu 0x=v(tu 0t). Since the cumulative count does not change along this characteristic, N(t;x) can be calculated as N(t;x) = minfN(t1;x1) +K(xx1);N(tu 0;xu 0)g (22) 4.2 Multi-region multiclass Newell's method The problem is to nd points ( ti;xi) such that N(ti;xi) =bi. As suggested in Section 4.1, these points form a system of linear equations as they are connected via characteristic lines of slopewi. With multiple regions,the calculation of the backwards-moving characteristic is 13 more complicated because there are multiple characteristics. The objective is to nd a point (ti;xi) on the characteristic of slope wi+1through point ( ti+1;xi+1) such that N(ti;xi) =bi. We rst discuss how to handle a downstream boundary condition (de ned for x=L), then discuss the changes needed for an initial condition (de ned for t= 0) in Section 4.2.3. 4.2.1 Downstream boundary condition For a downstream boundary condition, the objective is to nd Ti0 for alli2[1;I] such that t0=tIX i=1Ti (23) and x0=x+IX i=1Tiwi=L (24) wherex0=Lbecause backwards moving waves intersect the downstream boundary of L. The intermediate points where region iends are labeled ( ti;xi) and can be found by ti=tIX j=i+1Tj (25) and xi=x+IX j=i+1Tjwj (26) Proposition 4 formally proves that de ning ( ti;xi) throughTias in equations (25) and (26) achieves a point ( ti;xi) on the line of slope wi+1through point ( ti+1;xi+1). Proposition 4. Fori2[2;I1], point (ti1;xi1)is on the line with slope withrough point (ti;xi). Point (t;x)is on the line with slope wIthrough point (tI1;xI1). Proof. Point (tI1;xI1) is de ned by equations (25) and (26) as tI1=tIX i=ITi=tTI (27) and xI1=x+IX i=ITiwi=x+TIwI (28) ThereforettI1=TIandxxI1=TIwI, so (tI1;xI1) is on a line of slope wI through (t;x). Similarly, as de ned by equations (25) and (26), titi1=tIX j=i+1Tjt+IX j=iTj=Ti (29) 14 and xixi1=x+IX j=i+1TjwjxIX j=iTjwj=Tiwi (30) so (ti1;xi1) is on a line of slope withrough (ti;xi). Newell (1993) showed that along characteristics moving at speed wi, the cumulative count increases at the rate of Kx. The objective is for intermediate points ( ti;xi) to satisfy N(ti;xi) =bi (31) To achieve this, we impose the constraint c(t0;x0) +K(x0xi)bi (32) Theterm in constraint (32) comes from the fact that N(ti;xi) is the minimum of uncon- gested (upstream) conditions and congested (downstream) conditions. However if c(t0;x0) + K(x0xi)< bithenN(ti;xi)< bialso, which violates condition (31). Notice that if c(t0;x0)bi, then constraint (32) admits xi= 0 (meaningiP j=1Tj= 0). The intuitive in- terpretation is that if c(t0;x0)bi, then regions j= 1:::i are not part of the congested characteristic as the region boundary vehicle bipassedx=Lbefore the congested charac- teristic reaches x=L. To achieve a linear program, we assume that c(t0;L) is linear, i.e. a constant ow constraint, such as a capacity. A piecewise linear de nition of c(t0;L) can be handled by treating each piece separately via the inf-morphism property. When (ti;xi) is determined by upstream (uncongested) conditions, constraint (32) may not hold with equality. Since N(t;x) is constant along uncongested characteristics of speed v(Newell, 1993), we can nd the smallest value ^tisuch thatc^ti;0 =bi, then trace a line of slopevto (ti;xi): N(ti;xi)c^ti;0 =bi (33) with ^ti=tixi v(34) As discussed in Section 4.1.3, the correct value of tiis the larger one. Combining equations (33) and (34) yields ti^ti+xi v(35) meaning that vehicle bican intersect the line of slope wi+1at timetior greater (if congested conditions determine bi's trajectory). If bidoes not cross the line x= 0 during t2[0;1), there should be a maximum value ^ xi2[0;L] such that c(0;^xi) =bi. (If ^xi<0, then ^ti>0. If ^xi>L, then region iis never active.) These two possibilities are illustrated in Figure 4. If we have an ^ xi2[0;L] and not a ^ti2[0;1), then equation (33) can be rewritten as N(ti;xi)c(0;^xi) =bi (36) 15 𝑥 𝑡𝑏2 𝑣 𝑣𝑏1𝐿 Ƹ𝑡2,0(0,ො𝑥1)𝑡1,𝑥1 𝑡2,𝑥2Figure 4: Illustration of boundary vehicles crossing the initial or upstream conditions. Equation (36) is equivalent to tixi^xi v(37) becausevtiis the minimum time for bito traverse the distance xi^xiat free ow speed v. SinceN(ti;xi) is de ned by the minimum value of Nfrom tracing the forward and back- ward characteristics, the objective is to nd the minimum value of tisubject to constraints (32) and (35) or (37) (whichever is appropriate). This leads to the following linear program: minIX i=1ti (38a) s:t: t i=tIX j=i+1Tj 8i2[0;I] (38b) xi=x+IX j=i+1Tjwj 8i2[0;I] (38c) x0=L (38d) c(t0;x0) +K(x0xi)bi8i2[1;I1] (38e)( ti^ti+xi vif^ti2[0;1) tixi^xi vif ^xi2[0;L]8i2[1;I1] (38f) 16 Ti0 8i2[1;I] (38g) Constraint (38f) is noteworthy as it encapsulates two constraints. As discussed in the deriva- tion of equations (35) and (37), either a ^ti2[0;1) or a ^xi2[0;L] de nes the point at which vehiclebipasses the initial or upstream boundaries. Whichever is appropriate should be used. Once intermediate points ( ti;xi) are known, then N(t;x) can be determined as the min- imum of several constraints. First, N(t;x) is upper-bounded by the forward-moving charac- teristic of speed vpassing through point ( t;x). N(t;x)c tx v;0 (39) The calculation of the backwards moving characteristic is more complicated because ( ti;xi) has two possible values: N(ti;xi) =bi, orxi=LwithN(ti;xi)> bi. In the latter case, Ti= 0 (region iis not active). Their e ect on N(t;x) can be summarized as follows. If xI1<L, then N(t;x)bI1+K(xI1x) (40) because along characteristics moving at speed wI1in regionI1, the cumulative count increases at the rate of Kx(Newell, 1993). Regardless of which regions are active, N(t;x)c(t0;x0) +K (xI1x0) +IX i=0K(xixi+1)! =c(t0;x0) +K(x0x) (41) Unfortunately, if Iregions are considered but N(t;x)< bI, then constraint (38e) will cause linear program (38) to be infeasible. One possible algorithm for avoiding the infeasi- bility is to solve J?= maxfJ2[1;I] : linear program (38) is feasible g. This requires solving a series of linear programs for each point ( t;x). Instead, we suggest an alternative by allow- ing constraints (38e) and (38f) to be exceeded, but penalizing the excess. Rewrite constraint (38e) as c(t0;x0) +K(x0xi)bii (42) wherei0 is the amount by which constraint (38e) is exceeded. Similarly, rewrite con- straints (38f) as( ti^tixi viif^ti2[0;1) tixi^xi viif ^xi2[0;L](43) Then objective (38a) is modi ed to minIX i=1ti+ IX i=1i (44) If is suciently large, then i= 0 is optimal if linear program (38) is feasible. If not, then the revised linear program will be feasible but the optimal solution will have i>0, which 17 will indicate that ( ti;xi) was not found because ( t;x) occurs before region i. Linear program (38) can thus be written as minIX i=1ti+ IX i=1i (45a) s:t: t i=tIX j=i+1Tj 8i2[0;I] (45b) xi=x+IX j=i+1Tjwj 8i2[0;I] (45c) x0=L (45d) c(t0;x0) +K(x0xi)bii8i2[1;I1] (45e)( ti^tixi viif^ti2[0;1) tixi^xi viif ^xi2[0;L]8i2[1;I1] (45f) Ti0 8i2[1;I] (45g) i0 8i2[1;I] (45h) The number of variables in linear program (45) scales linearly with the number of class proportion regions considered. After solving linear program (45), N(t;x) is the maximum value satisfying several constraints, including (39) and (41). Constraint (40) is modi ed as follows: nd the largest active region ~Iwith~I= 0, i.e. ~I= maxfi2[1;I] :i= 0g (46) Ifi>0, then constraints (38e) or (38f) do not hold, meaning that region ioccurs after point (t;x). Then if ~I >0, N(t;x)b~I1+K x~I1x (47) 4.2.2 Piecewise-linear downstream boundary conditions Linear program (45) can be modi ed for a piecewise-linear downstream boundary condition by solving with each piece separately with the additional constraint t[t0t](48) wheret[andt]are the lower and upper bounds of tfor which the condition applies, re- spectively. Then, using the inf-morphism property (Claudel and Bayen, 2010a), the correct solution is the one that obtains the lower N(t;x) value. However, two additional issues arise due to the addition of constraint (48). First, the value of ( ti;xi) in the optimal solution may not have either (45e) or (45f) as an active constraint, i.e. N(ti;xi)>bi, to achieve feasibility with constraint (48). This issue is illustrated in Figure 5(a). If w2> w 1, the correct value 18 Region 2𝑡1,𝑥1 −𝑤2−𝑤1 𝑡,𝑥 Region 1𝑥 𝑡𝐿𝑡♭𝑡♯(a)N(ti;xi)>bi Region 2𝑡1,𝑥1−𝑤2−𝑤1 𝑡,𝑥 Region 1𝑡♭𝑡♯ 𝑥 𝑡𝐿 (b)Ti>0 withi>0 Figure 5: Illustration of issues caused by the addition of constraint (48). Dotted lines indicate the optimal characteristics without constraint t02 t[;t] . of (t1;x1) may be too high for a line of slope w1to reach the boundaries of the constraint t0t]. In this case, constraint (40) should only be applied when constraints (45e) or (45f) are active. Second, due to the shallowness of a slope wi, it is possible to have Ti+1>0 while i+1>0, which means that N(ti;xi)<bi. Figure 5(b) illustrates this as N(t1;x1)<b 1yet becausew1< w 2,T2>0 is needed to achieve t0t[. In this case, ( t0;x0) should not be achievable with constraints (38e) and (38f) in place. This results in a incorrect (too small) value ofc(t0;x0) being applied to constraint (41). In this case, constraint (41) on N(t;x) should not be used. 4.2.3 Initial condition Linear program (45) can be used with a slight modi cation to handle an initial condition (de ned for t= 0). Replace constraint (45d) with the following: t0= 0 (49) The remainder of the constraints are suciently general to hold for both a downstream and initial condition. 4.3 Numerical demonstration of multiclass Newell's method We demonstrate the multiclass calculation of N(t;x) by considering a 2-class scenario of legacy and automated vehicles. Using the triangular ow-density relationship derived by Levin and Boyles (2016), the congested wave speed for region iwith an automated vehicle proportion of pAV iis wi=` pAV i AV+ (1pAV i) HV(50a) 19 with capacity of Q=v v(pAV i AV+ (1ipAV i) HV) +`(50b) where HVand AVare the reaction times of legacy and automated vehicles, respectively, and`=1 Kis the vehicle length. For this demonstration, we used speci c values of HV= 1:5s, AV= 0:25s,v= 30mi/hr, and K= 240veh/mi. We consider a 1-mile length of road with a red trac signal (0 ow) from t= 30s tot= 210s. We de ned two regions: region 1, for vehicles [0 ;50], with 0% automated vehicles, and region 2, for vehicles (50 ;1), with 100% automated vehicles. Therefore, capacities for regions 1 and 2 are 1800vph and 4800vph, respectively. The initial condition c(0;x) was set at a density of 40veh/mi, with c(0;x) = 40kx. Upstream entering ow was 1200veh/hr for t2[0;60) then 2400veh/hr fort2[60;1). These parameters are not intended to suggest representative values for automated vehicle ow. Rather, the purpose is to provide a demonstration of the multiclass Newell's method for solving N(t;x). The solution method described in Section 4.2 was implemented in Java, using IBM's CPLEX 12.6 to solve linear programs. At a time interval of 1s and a spatial interval of 0.01mi, computing the cumulative count map for t2[0;360]s andx2[0;1]mi required 41.48s on a desktop computer with an Intel Core i5{8600K at 3.60GHz with 16GB of memory. Comput- ing the density map required 77.34s. Density was approximated by k(t;x)N(t;x)N(t;x+) . The cumulative counts are shown in Figure 6(a) and the densities are shown in Figure 6(b). There are several interesting patterns of note. First, from t= 60s tot= 360s, entering density is 80veh/mi, which exceeds the critical density of region 1 (60veh/mi) but is well within the uncongested regime of region 2. Consequently, the shockwave separating the queue from the entering ow changes speed at (155s ;0:79mi). Downstream, and after the trac signal turns green at t= 210s, the queue initially dissipates at a rate of 1800veh/hr, with a density of 60veh/mi, until the 50th vehicle has passed. Thereafter, the queue dissipates at the rate of 4800veh/hr, or a density of 160veh/mi, due to the higher capacity of region 2. 20 (a)N(t;x) (b)k(t;x) Figure 6: Illustration of computed cumulative count map and density map from multiclass trac ow scenario 5 Multiclass link transmission model We now apply the multiclass Newell's method into a multiclass version of Yperman et al. (2005)'s link transmission model for dynamic network loading. This paper focuses on de- veloping a link transmission model for a mixture of human and AV ows with a changing ow-density relationship. Dynamic trac assignment for a single-class of vehicles is well- established in the literature, and the methods in this paper can be used with existing node models (Tamp ere et al., 2011) and dynamic trac assignment algorithms (Levin et al., 2015). For more information on dynamic network loading, we refer the reader to Yperman (2007) which describes the single-class link transmission model for link ow modeling as well as its connection to node models. In Section 6, we will demonstrate a dynamic trac assignment using the multiclass link transmission model on a city network, combined with standard node models for trac signals, merges, and diverges. As with other link models, the multiclass link transmission model must calculate sending and receiving ows, which are then used by the node models to calculate entering and exiting ow and the upstream and downstream ends of each link. The link transmission model stores the upstream and downstream boundary conditions up to time t. Cumulative count values within the link are not stored. 5.1 Class proportion regions Consider discretized time with time step  t. The class proportion regions are de ned sep- arately for each link. The class proportion region boundaries could be chosen in several ways. To work with both continuous and discrete ow, we de ne the region boundaries to 21 coincide with time steps. Flow entering between time tandt+ tis de ned as one class proportion region, with boundaries N(t;0) andN(t+ t;0). We assume that Nm(t;0) and Nm(t+ t;0) are known, then we can calculate the average entering ow between [ t;t+ t) as qm(t;0)Nm(t+ t;0)Nm(t;0) t(51) We then create a new class proportion region applying to the vehicles entering between t andt+ twith proportions ofqm(t;0)P m02Mqm0(t;0). This implicitly assumes that class proportions are uniformly mixed between tandt+ t. If desired, a higher resolution of class proportion regions could be implemented. However, when working with continuous ows, class proportions at the upstream end of each link between time tandt+ tare assumed to become uniformly distributed. This is not a problem at centroids, but violates FIFO at the upstream ends of internal links. Although there are fewer mixing points in this LTM than the multiclass CTM (Levin and Boyles, 2016), mixing still occurs. This is actually a general problem with dynamic network loading of multi-commodity ow (Carey et al., 2014). For other papers, the commodities are typically paths or destinations, and vehicle class becomes an additional identi er in this paper. Bar-Gera and Carey (2017) studied methods of ensuring FIFO holds in continuous dynamic network loading models. Although addressing this issue is outside the scope of this paper, it is worth studying in future work. 5.2 Sending ow The sending ow is the maximum ow that could exit in one time step, and the receiving ow is the maximum ow that could enter in one time step. The sending ow S(t) can be calculated by S(t) =N(t+ t;L)N(t;L) (52) whereN(t+ t;L) is the maximum value found based on constraints within the link. (The actual ow may be less than the sending ow depending on downstream constraints from the node model.) Although N(t;L) is stored, N(t+ t;L) must be calculated. Since L is the downstream boundary, the only characteristic that can reach a boundary condition is the uncongested characteristic, which always travels at speed vregardless of the class proportions. Therefore, N(t+ t;L) =N t+ tL v;0 (53) Due to equation (53), onlyL vtime steps of the upstream cumulative count need to be stored. Sending ow is also constrained by capacity, which is not a constant value because it depends on the class proportions. Given class proportion regions 0 :::I, the rst class proportion region of consideration, i[, is the rst region active on the time interval [ t;t+t): i[= arg max i2[0;I]fbi1N(t;L)g (54) 22 Then the capacity constraint ~Qcan be constructed using the following algorithm. Add up Qi[of vehicles until either the class region boundary bi[or the time step  tis reached. If the class region boundary is reached, then continue with the next region i[+ 1 until 1 time step has elapsed. This calculation is de ned by Algorithm 1. ~Qand ~ are state accumulation variables. ~bis the region lower bound, which could be bi1orN(t;L). For each class region i, that region can apply for up to =bi~b Qitime (after which the class region transitions to i+ 1). If the remaining time  t~ > , then addQi to the capacity variable ~Q, update ~b and ~ , and continue with region i+ 1. Otherwise, add Qi(t~ ) capacity and terminate. The sending ow calculation can then be summarized by S(t) = min ~Q;N t+ tL v;0 N(t;L) (55) Algorithm 1 Calculating capacity for sending ow 1:~Q 0 2:~b N(t;L) 3:~ 0 4:fori=i[toIdo 5: =bi~b Qi 6:if~ + tthen 7: ~Q ~Q+Qi(t~ ) 8: return ~Q 9:else 10: ~Q ~Q+Qi 11: ~ ~ + 12: ~b bi 13: end if 14:end for 5.3 Receiving ow Like the sending ow, the receiving ow R(t) is calculated as R(t) =N(t+ t;0)N(t;0) (56) N(t+ t;0) requires tracing a congested characteristic, which has variable speed. The method for tracing this characteristic is de ned in Section 4.2. However, there are some useful implementation results to be discussed here. First, we prove in Proposition 5 that onlyL w[time steps of downstream cumulative counts need to be stored, which determines the memory requirements. Second, observe that if a class region iends with boundary vehicle bisuch thatN t+ tL w[;L bi, then class region iis no longer relevant to calculating N(t+ t;0), and can be discarded from memory. The capacity constraint for the receiving ow can be estimated based on the class proportions of upstream link sending ows. 23 Proposition 5. Downstream cumulative counts must be stored for the previousL w[ttime, wherew[is the minimum congested wave speed. Proof. To calculate N(t+ t;0), the multiclass Newell's method follows a congested char- acteristic of speeds wito a timet+ tIP i=1Tisuch thatIP i=1Tiwi=L. Sincewiw[,IP i=1Ti can be upper-bounded byIP i=1TiL w[therefore t+ tIX i=1Tit+ tL w[(57) which requires looking back at mostL w[tin time. 5.4 Iterative algorithm for receiving ow The repeated use of solvers to calculate N(t+t;0) can limit the computational performance of large dynamic network loading models. In this section we present a polynomial time algorithm based on iteratively solving a system of linear equations. The multiclass link transmission model bene ts from certain simpli cations that are not true in the general multiclass Newell's method. First, since only the sending and receiving ows are calculated, N(t+ t;x) is calculated for x=Lorx= 0. When x=L, the solution is simply looking backwards along an uncongested wave of speed v. Whenx= 0, the solution involves looking backwards along a congested wave with varying speed. However, due to the construction of the region boundaries, N(t+t) is in region I+ 1 where region Ihas boundary vehicle bI= N(t;0). The second simpli cation is therefore that the class boundary region is known with certainty when calculating N(t+t;0). This suggests that N(t+t;0) can be calculated by iteratively solving a system of linear equations. The iteration occurs because the intermediate points (ti;xi) satisfying N(ti;xi) =bicould be determined either by congested or uncongested conditions. The rst iteration can assume that all ( ti;xi) points are determined by congested conditions, then revise them to use uncongested conditions if they result in the contradiction thatN(ti;xi) =biyetN(ti;xi)c tixi v;0 . We assume here that I3. IfI < 3, then the method of Section 4.1 may be used instead. For the rst iteration, assume that congested conditions determine all points ( ti;xi). Therefore, N(t+ t;0) is determined by following a congested characteristic back to the downstream boundary condition c(t;L). This results in the following system of equations: c(t0;L) +K(Lx1) =b1 (58a) K(xixi+1) =bi+1bi 8i2[1;I1] (58b) w1(t1t0) =x1L (58c) wi+1(ti+1ti) =xi+1xi8i2[1;I1] (58d) wI+1(t+ ttI) = 0xI (58e) 24 The simpli ed solution involves removing equations (58b), (58c), and (58d) to obtain a system of equations in only three variables: t0,tI, andxI. These equations can be removed because the di erence in cumulative count between boundary vehicles bi+1andbiis constant, soK(xixi+1) is also constant. Since wi+1is given, the di erence ti+1tican be determined as ti+1ti=xi+1xi wi+1=bi+1bi Kwi+1(59) Equations (58a) and (58b) can be combined into c(t0;L) +K(LxI) =bI (60) Equations (58b) and (59) can be used to nd tIt0as tIt0=b1c(t0;L) Kw 1+IX i=2bibi1 Kwi(61) Whenc(t0;L) is linear (as it is for LTM), equations (58e), (60), and (61) form a system of three equations with variables t0,tI, andxI. Assume that c(t0;L) =c0+q0t0. Solving this system results in the solution t0=bIc0KL+KwI+1(t+t) KwI+1+c0b1 Kw1IP i=2bibi1 Kwi q0 KwI+1+ 1q0 Kw1(62) withtIandxIfound through equations (61) and (60), respectively. Using equations (58c) and (59), all intermediate points ( ti;xi) can be found. The prior solution to ( ti;xi) points assumed that all points were determined by congested conditions. It is easy to check whether any given ( ti;xi) point should be determined by uncongested conditions instead: c tixi v;0 <bi. Let ^iindicate the furthest vehicle position determined by uncongested conditions, i.e. ^i= maxn i2[1;I] :c tixi v;0 <bio (63) Then the true value of ( t^i;x^i) is on a line of slope vfrom a point t0such thatc(t0;0) =b^i. Givenc(t;0), thatt0can be uniquely determined and can be hereafter treated as a constant. This works because the cumulative count does not change along a uncongested characteristic and is therefore independent of the length of that characteristic. We now have a revised system of linear equations: v(t^it0) =x^i0 (64a) 0xI=wI+1(t+ ttI) (64b) K(xixi+1) =bi+1bi 8i2h ^i;I1i (64c) xi+1xi=wi+1(ti+1ti) 8i2h ^i;I1i (64d) 25 As before, we want to simplify equations (64) to have only variables t^i,x^i,tI, andxI. Equations (64c) can be simpli ed to K(x^ixI) =bIb^i (65) Using equation (59), equation (64d) can be simpli ed to tIt^i=IX i=^i+1bibi1 Kwi(66) Equations (64a), (64b), (65), and (66) now form a system of four equations with four vari- ables. The solution is t^i=bIb^i K+vt0+wI+1(t+ t)wI+1IP i=^i+1bibi1 Kwi v+wI+1(67) withtI,xI, andx^ifound by equations (66), (64b), and (65), respectively. Intermediate (ti;xi) points can be calculated by equations (64c) and (64d). As before, if some of these intermediate points should be instead determined by uncongested conditions, calculate a new ^iusing equation (63) and solve equations (64) again. The algorithm to calculate N(t+t;0) using iterative systems of linear equations is summarized in Algorithm 2. Algorithm 2 Iterative algorithm for calculating the receiving ow for the multiclass LTM 1:Set^i 1 2:ifI+ 1^i1then 3: SolveN(t;0) using Newell (1993)'s method (1 region) 4:else ifI+ 1^i2then 5: SolveN(t;0) using 2-region multiclass Newell's method (Section 4.1) 6:end if 7:if^i= 1then 8: Solve equations (58e), (60), and (61) for N(t;0) 9:else 10: Solve equations (64a), (64b), (65), and (66) for N(t;0) 11:end if 12:if i2[1;I] :c tixi v;0 <bi =?then 13: returnbI+KxI 14:else 15: ^i max i2[1;I] :c tixi v;0 <bi 16: go to line 2 17:end if 26 Figure 7: Downtown Austin city network 6 Numerical results from dynamic trac assignment We implemented the multiclass LTM in a dynamic trac assignment model that discretizes ow and tracks the location, path, and driver type (legacy or automated) of individual vehicles with a time step of  t= 15s. By discretizing vehicles, we can reduce the FIFO issues discussed in Section 5.1. The sending and receiving ows calculated by the multiclass LTM are used to determine the propagation of vehicles along links and through nodes. This same software package (albeit with a di erent link ow model) has been used for numerical results in previous studies (e.g. Levin and Boyles, 2016; Levin et al., 2017). We conducted a demonstration of dynamic trac assignment using the multiclass LTM on the downtown Austin city network, shown in Figure 7. This network was constructed by the Network Modeling Center at The University of Texas at Austin for project work for the City of Austin, and was calibrated to match observations from 2010. The network includes 171 zones, 546 nodes, and 1247 links, and has a demand of around 62,836 vehicles over 2 hours. 6.1 Convergence properties of dynamic trac assignment We solved dynamic trac assignment using the method of successive averages. Normal convergence of dynamic trac assignment was observed with the multiclass LTM in use. For instance, Figure 8 shows the gap with respect to iteration for a scenario with 50% AVs. 27 0102030405060708090100 0 5 10 15 20 25 30Gap (%) IterationFigure 8: Convergence of dynamic trac assignment using the multiclass LTM This gap is de ned as gap =TSTTSPTT TSTT(68) where TSTT is the total system travel time and SPTT is the shortest path travel time. A gap of under 2% was achieved after 30 iterations, which is typical for this network and for dynamic trac assignment (Levin et al., 2015). We note that due to the discontinuities in travel times caused by LTM, the existence and uniqueness properties of dynamic user equilibrium cannot be analytically established. The average computation time per iteration was 19.2s on a desktop with an Intel Core i5-9400 processor clocked at 2.90 GHz with 24.0GB of memory. Most of the additional computational e ort of Algorithm 2 compared with LTM occurs when the number of class proportion regions is large. E orts to reduce the number of regions whenever possible are important for achieving a lower computation time. For example, one new region is usually created for each link per time step. However, if the entering ow is zero, then a new region does not have to be created for that link. If the entering ow for a link has a AV proportion identical to the last region, then a new region also does not need to be created there. 6.2 E ects of AV market penetration on trac congestion Next, we evaluated the sensitivity of network travel times with respect to AV market pene- tration. Link capacities and congested wave speeds were scaled according to equation (50). To retain as much of the original network parameters as possible, link capacities were not set equal to the output of equation (50b). Instead, we calculated QforpAV iand with 0 AV 28 market penetration, then scaled the calibrated link capacities proportionally. Formally, if ~Qais the link capacity of link a, andQ pAV i is the capacity calculated by equation (50b), then the new link capacity for region iis chosen as ~QaQ(pAV i) Q(0). Congested wave speeds were then determined by the triangular shape of the ow-density relationship, with a xed jam density. We used reaction times of HV= 1s and AV= 0:5s. For each AV market pene- tration, 10 Monte Carlo simulations were performed to reduce the e ects of stochasticity in the demand. Because ow is discretized, each vehicle is randomly determined to be either a legacy or automated vehicle with probability given by the overall AV market penetration. The results presented here assume that intersection capacity increases with the AV market penetration as speci ed by the ow-density relationship (50). It is possible that intersection capacity may be higher for AVs due to signal-free intersection controls (Dresner and Stone, 2004), but such controls were not included in these results. It is also possible that AVs might instead behave more cautiously at intersections, resulting in lower intersection capacities than predicted here. These results are also highly dependent on the shape of the ow-density relationship for AVs, and it is possible that AVs might reduce capacity instead (Shladover et al., 2012; Calvert et al., 2017; James et al., 2019). Figure 9 shows the average travel time as a function of the city-wide AV market pen- etration. Error bars show the standard deviation of travel times from the Monte Carlo simulations. Since legacy and automated vehicles share all links, the average travel times do not vary based on the driver type. Interestingly, the average travel times remained rela- tively constant between 10% and 50% AVs. We note that Braess (1968) and Daganzo (1998) demonstrated that increases in capacity could result in higher travel times due to sel sh route choice, which was actually demonstrated for cooperative adaptive cruise control by Melson et al. (2018). The lack of improvement could also be due to intersection capacity limitations and queue spillback. Although in these results, the intersection capacity increases with the time- and space-dependent AV market penetration, queue spillback could reduce the ow improvements achieved by capacity increases at intersections. At 60% and higher AV market penetrations, travel times decreased consistently. Of course, the e ects that AVs have on congestion depend both on the network and the legacy and automated vehicle driving parameters. The results in Figure 9 are compared with other results in the literature. For freeways, Shi and Prevedouros (2016) and Yu et al. (2019) observed a similar pattern as shown in Figure 9: speeds increased only slightly at low market penetrations but had an increasing impact at higher market penetrations. However, they predicted more continuous e ects than shown in Figure 9. Intersection controls could be causing the di erence. Ghiasi et al. (2017) predicted even larger increases in capacity at low market penetrations, but without network results. Patel et al. (2016) found that congestion decreased steadily with increasing AV proportions even at low market penetrations, but Levin and Boyles (2016) observed that trac signals were a major bottleneck for the downtown Austin network. 29 02468101214161820 0 10 20 30 40 50 60 70 80 90 100Average travel time (min) AV proportion (%)Figure 9: E ects of AV market penetration on travel times 7 Conclusions This paper extended Newell (1993)'s method for a multiclass kinematic wave theory. Like the original, the proposed method nds exact solutions under triangular ow-density rela- tionships. We started by exploring the multiclass kinematic wave theory, then recognizing that when vehicles obey rst-in- rst-out behavior, that the class-speci c cumulative counts Nm(t;x) can be adequately described by a single aggregate cumulative count N(t;x) and class proportion regions that move with the trac ow. We wrote a corresponding mul- ticlass Lax-Hopf formula (Claudel and Bayen, 2010a). Restricting ourselves to triangular ow-density relationships with varying capacities and congested wave speeds, we developed a linear program to implement a multiclass Newell (1993)'s method to nd exact solutions to the multiclass kinematic wave theory. With the ultimate goal of a multiclass link transmission model, we developed a more ecient iterative algorithm for calculating receiving ows. The multiclass link transmission model was demonstrated on the downtown Austin city network. Dynamic trac assignment was observed to converge normally, and computation times on this city network were easily manageable. Numerical results show a non-linear decrease in overall congestion due to greater AV market penetration. There are many future opportunities for this research. The multiclass LTM could be used for more accurate city-scale predictions of trac ow and congestion with varying proportions of automated vehicles to prepare for their future use on public roads. The presented solution method is fairly complex, and it may be possible to use the problem structure to achieve a more ecient solution method. The ow-density relationship for 30 autonomous vehicles may not follow a triangular shape, so extending the results in Sections 4 and 5 to more general ow-density relationships would be valuable. It may be possible to use Lagrangian coordinates and/or variational theory (Daganzo, 2006; Leclercq et al., 2007; Laval and Leclercq, 2013) to develop a more ecient solution method. The predictions of this model should also be veri ed against microsimulation. Multiclass ow that does not strictly adhere to the rst-in- rst-out property on links is more challenging to model, and could be the subject of future work. Although the standard LTM (Yperman et al., 2005) has been used in system optimal dynamic trac assignment (Levin, 2017; Chakraborty et al., 2018), the complexity of the multiclass Newell's method makes formulating the system optimal dynamic trac assignment problem challenging for multiclass ow. Finally, other multiclass kinematic wave theory problems might be able to adapt the concepts and methods presented here. Acknowledgements The authors gratefully acknowledge the support of the National Science Foundation, Award No. 1935514.