CN110472353B

CN110472353B - Traffic network design method based on user utility maximization

Info

Publication number: CN110472353B
Application number: CN201910771290.2A
Authority: CN
Inventors: 林宏志; 尹璐
Original assignee: Southeast University
Current assignee: Southeast University
Priority date: 2019-08-20
Filing date: 2019-08-20
Publication date: 2023-07-18
Anticipated expiration: 2039-08-20
Also published as: CN110472353A

Abstract

The invention provides a traffic network design method based on user utility maximization. Algorithmically, this is a two-layer model, with the upper layer being used to maximize the utility of the network user at a given investment budget, and the lower layer being the feedback process between the Nested Logit model and the user equalization model. The travel demand prediction adopts a discrete selection model, so that the utility of the network user can be measured. However, the solution of the bilayer model is very challenging, the global optimal solution is difficult to obtain by the traditional optimization algorithm, and the software designs a simulated annealing algorithm to find the optimal solution. From the development language, the R language is free of charge, and the matrix language algorithm is high in efficiency, so that the software is easy to use for scientific researchers and industries. Simulation researches show that the invention can effectively find the traffic network design with the maximum utility of network users under a given investment budget.

Description

Traffic network design method based on user utility maximization

Technical field:

the invention provides a traffic network design method based on user utility maximization, and belongs to the technical field of traffic engineering.

The background technology is as follows:

user-centric mobile services are one of the emerging mobile services in the future. Unlike traditional vehicle-centric mobile services, it is claimed a mobile solution that is centered on improving the user experience. While there are many ways in which the user experience can be measured, user utility is the most common indicator.

The consumer remains widely used in the design of traffic networks for elastic travel needs. Yang and Bell ^[1] Consumer residuals are used as an objective function of road network design. Since then, this measurement method has been widely used. Yang and Meng ^[2] A two-layer planning model is proposed to determine the optimal throughput and toll for newly built roads, thereby maximizing customer and producer residuals. Szeto and Lo ^[3，4] Three different government road network design strategies were compared, wherein network benefits were measured in terms of the total amount of discounts left by the consumer. Chen et al ^[5] A genetic algorithm is proposed to solve the problem of construction-operation-handover (BOT) network design when demand is uncertain, where the governmental goal is to maximize the desired consumer surplus and minimize its variance. Ukkusuri and Patil ^[6] ConsiderTo the uncertainty and resilience of the requirements, a consumer residual maximization model of the network design problem is established. Jiang and Szeto ^[7] Decision makers are considered to aim at increasing consumer residuals and reducing health costs when designing time-varying traffic networks. Zhang et al ^[8] Aiming at the road network design problem, a double-layer planning model is provided, wherein the upper layer is used for realizing the residual maximization of consumers, and the lower layer adopts random user balance of elastic requirements.

User utility is based on the individual's perception of a product or service. In general, the expected maximum utility (log sum) is an indicator of consumer utility in a Logit selection model. Although discrete selection models are often used in transportation, log sum is rarely used in project assessment ^[9] . Dong et al ^[10] Consumer utility based utility is suggested to replace the traditional consumer residue. Bhandari et al ^[11] When introducing subway systems in de, logsum is used to estimate the utility of the transfer. Geurs et al ^[12] The expected maximum utility is used as an accurate measure of the benefit of traffic users. Tillema et al ^[13] It is suggested to evaluate the impact of urban congestion charging from a log sum-based consumer utility perspective. Nahmas-Biran and Shibatan ^[14] The proposal can use logsum as a measure of user utility in traffic item assessment. In addition, fairness should also be considered so that resources are allocated more reasonably. Van Wee ^[15] The advantages and disadvantages of the logsum method are discussed and more research into this problem is desired, especially with respect to fairness and social rejection problems. Javanmardi et al ^[16] The service network design problem of introducing utility functions into objective functions is studied.

Reference is made to:

[1]Yang，H.，Bell，M.G.H.(1998)Models and algorithms for road network design：a review and some new developments[J].Transport Reviews 18，257-278.

[2]Yang，H.，Meng，Q.(2000)Highway pricing and capacity choice in a road network under a build-operate-transfer scheme[J].Transp.Res.Pt.A-Policy Pract.34，207-222.

[3]Szeto，W.Y.，Lo，H.K.(2005)Strategies for road network design over time：Robustness under uncertainty[J].Transportmetrica 1，47-63.

[4]Szeto，W.Y.，Lo，H.K.(2008)Time-dependent transport network improvement and tolling strategies[J].Transp.Res.Pt.A-Policy Pract.42，376-391.

[5]Chen，A.，Subprasom，K.，Ji，Z.W.(2006)A simulation-based multi-objective genetic algorithm(SMOGA)procedure for BOT network design problem[J].Optim.Eng.7，225-247.

[6]Ukkusuri，S.V.，Patil，G.(2009)Multi-period transportation network design under demand uncertainty[J].Transp.Res.Pt.B-Methodol.43，625-642.

[7]Jiang，Y.，Szeto，W.Y,(2015)Time-dependent transportation network design that considers health cost[J].TransportnetricaA 11，74-101.

[8]Zhang，X.，Wang，H.，Wang，W.(2015)Bi-level programming and algorithms for stochastic network with elastic demand[J].Transport 30，117-128.

[9]de Jong，G.，Daly，A.，Pieters，M.，van der Hoorn，T.(2007)The logsum as an evaluation measure：Review of the literature and new results[J].Transp.Res.Pt.A-Policy Pract.41，874-889.

[10]Dong，X.J.，Ben-Akiva，M.E.，Bowman，J.L.，Walker，J.L.(2006)Moving from trip-based to activity-based measures of accessibility[J].Transp.Res.Pt.A-Policy Pract.40，163-180.

[11]Bhandari，K.，Kato，H.，Hayashi，Y.(2009)Economic and Equity Evaluation of Delhi Metro[J].International Journal of Urban Sciences 13，187-203.

[12]Geurs，K.，Zondag，B.，de Jong，G.，de Bok，M.(2010)Accessibility appraisal of land-use/transport policy strategies：More than just adding up travel-time savings[J].Transport.Res.Part D-Transport.Environ.15，382-393.

[13]Tillema，T.，Verhoef，E.，van Wee，B.，van Amelsfort，D.(2011)Evaluating the effects of urban congestion pricing：geographical accessibility versus social surplus[J].Transp.Plan.Technol.34，669-689.

[14]Nahmias-Biran，B.H.，Shiftan，Y.(2016)Towards a more equitable distribution of resources：Using activity-based models and subjective well-being measures in transport project evaluation[J].Transp.Res.Pt.A-Policy Pract.94，672-684.

[15]Van Wee，B.(2016)Accessible accessibility research challenges[J].Journal of Transport Geography 51，9-16.

[16]Javanmardi，S.，Hosseini-Nasab，H.，Mostafaeipour，A.，Fakhrzad，M.，Khademizare，H.(2017)Developing a New Algorithm for a Utility-based Network Design Problem with Elastic Demand[J].Int.J.Eng.30，758-767.

the invention comprises the following steps:

technical problems: policy makers and researchers generally want an analysis technique that is less dependent on set point analysis, but has more behavioural mechanisms. Therefore, the utility centering on the user is maximized as an evaluation index of the traffic network design, and a new traffic network design method is provided.

The technical scheme is as follows: the invention provides a traffic network design method based on user utility maximization, which is a double-layer model, and specifically comprises the following steps:

(one) concrete method for establishing upper layer model

According to the random utility theory, the Nested Logit model probability can be written as the product of two Logit probabilities. The first layer is a binary choice of whether to go out. The second layer is to make multiple selections in the destination of the trip. The utility obtained from destination selection s by decision maker n residing at departure point r is:

wherein the method comprises the steps ofDepending on the variables describing whether or not to travel, +.>Depending on the variable describing destination s, +.>Is the extremum distribution of independent same distribution. Usually, is->And->Are designated as linear parameters. There may be many interpretation variables, depending on the availability of the data. However, OD travel time is recognized as the most fundamental one. Therefore, it can be added with->The expression is as follows:

wherein beta is ₀ Is a constant term of the number of the items,is of coefficient beta ₁ The shortest travel time between OD versus rs, +.>Is a vector of other observable variables having coefficients β.

The probability of selecting a destination s is the product of two probabilities, namely the probability of going out and the probability of determining that the destination s is selected after going out:

wherein the method comprises the steps ofIs the marginal probability of going out, +.>Is the conditional probability of selecting the destination s after determining the trip. The marginal conditional probability and the conditional probability take the form of a binomial logic and a polynomial logic, respectively, which can be expressed as:

wherein the method comprises the steps of

And S is _r Is the set of destinations of origin r.Is the desired maximum utility that decision maker n derives from the destination selection. />Because its expression is also called logsum term. Lambda is commonly referred to as the logsum coefficient. Lambda reflects the degree of independence between destinations, with greater lambda exhibiting greater independence and less correlation. Lambda is in the range of 0,1]. When λ=1 means that the choice between destinations is completely independent. In this case, the Nested Logit model can be reduced to a standard Logit model.

Note that if an individual n is considered to be representative n, the subscript n may be ignored. Potential travel demand O given origin r _r The travel demand at the departure point r can be expressed as:

Q _r ＝Q _r p _rT ， (7)

wherein p is _rT Obtained from equation (4). The travel demand between OD and rs can be definedThe method comprises the following steps:

wherein p is _rs Obtained from equation (3). Note that at this time the OD travel distribution matrixThe matrix is measured by a person and is not generally directly assigned to the road network, and is converted into units of passenger cars (pcu) using a vehicle loading factor mu. Thus, the travel distribution matrix in units of passenger cars can be expressed as:

q _rs ＝O _r p _rs /μ. (9)

where μ can be determined by investigation.

Since the goal of the traffic system is to provide as much service as possible to the user, the policy goal is to maximize the user's utility on the outbound journey. According to binomial selection formula (4) indicating whether to travel, the expected travel utility of the individual n residing at the departure point r can be expressed as:

since the absolute level of utility is not measurable, it can be divided by the marginal utility θ of revenue _n Thereby obtaining the outgoing utility of the user n residing at the departure place r:

this division converts units of utility into currency. Assuming that the individual residing in the traffic analysis region r has the same utility as the representative individual n, the subscript n is ignored. Consumer utility sum CW for resident out of living area r _r Equal toAnd the total travel demand Q of the district _r Product of two times：

While the user utility of the system is the set of all origin areas, expressed as

Combined solving travel demand Q _r The mathematical programming of the upper layer can be written as:

wherein B is the investment budget, g _a Is a cost function of the road segment a,which is the maximum increment allowed on road segment a. Other symbols are consistent with the previous definitions. The OD travel time is implicitly determined by the underlying model.

(II) specific method for building lower model

The cyclic feedback process of the underlying model is shown in fig. 2. In the Nested Logit model, from the initial solution of the OD travel time, a travel distribution matrix can be generated. And then, distributing the travel demands among the ODs to the road network by using a user balance model. The road travel time can thus be obtained. According to the Wardrop principle, a shortest path algorithm (typically Dijkstra algorithm) can be used to calculate the new OD travel time. For a continuous average, the new OD travel time is a weighted sum of the previous two OD travel times. The latest OD travel time is then entered into the new logic model. The iterative process will continue until the OD travel time for successive iterations is substantially equal. In other words, the OD travel time is internally determined and consistent.

The detailed steps of the continuous averaging Method (MSA) are as follows:

step 1, initializing. The feasible solution of OD travel time is expressed asIn addition, let i=1, which represents the number of iterations.

Step 2 combines the travel-destination selection modes. Travel distribution matrixIs generated by a Nested Logit model, and the parameters can be obtained in advance through demonstration research.

And 3, carrying out traffic distribution by using a user balance model. Given a traffic network, calculating the travel time of a road section by adopting a Frank-Wolfe algorithm

And 4, solving the shortest path travel time by using a Dijkstra algorithm. According to the Wardrop principle, the shortest path travel time can be used as a new travel time

Step 5, updating the OD travel time. MSA method averaging with decreasing weight (i.e. inverse of iteration number i)And

step 6 uses a predetermined relative root error (RRSE) epsilon to check the convergence of the OD travel time:

if the convergence condition is met, go to step 7. Otherwise, let i=i+1 andand then goes to step 2.

Step 7, ending the feedback process. Output consistent OD travel timeAnd travel distribution matrix->

(III) specific steps of simulated annealing algorithm

For complex bilayer problems, traditional methods are difficult to work with due to their non-linearity, non-microability, non-convexity, etc. In this case, the heuristic approach, although computationally demanding, can be used as an alternative tool to solve some of the difficulties described above and to obtain an approximate optimal solution. The invention employs a Simulated Annealing (SA) algorithm that is particularly suited for bilayer models. The algorithm flow chart is shown in fig. 3. The specific steps of the algorithm are as follows:

step 0, initializing: given a feasible solution set Ω, the parameters 0 < δ < 1.0 and the integer L ₀ ，M ₀ And T _s (stop temperature). In addition, let k=0 be the outer loop control parameter, k ₁ =0 is the inner loop control parameter.

Step 1, searching an initial temperature and an initial point: uniformly randomly generating M on a feasible solution set Ω ₀ Points, with z ⁽ⁱ⁾ (i＝1，2，...，M ₀ ) To represent. For each z ⁽ⁱ⁾ Using an underlying model to obtain a sum z ⁽ⁱ⁾ Related OD travel time and then calculate corresponding upper layer performanceI.e. the user utility. Next, calculate +.>The variance of (a) is the initial temperature T ⁽⁰⁾ While the initial point y ⁽⁰⁾ E.OMEGA.then being optimal->The point to which the value corresponds.

Step 2, checking the stopping condition of the outer loop: if T ^(k) ＜T _s Then stop. Otherwise, go to step 3.

Step 3, checking the stop condition of the internal circulation: if k is ₁ ＞L ₀ N (where N is the number of decision variables in the upper model), then go to step 6. Otherwise, go to step 4.

Step 4, generating a feasible point: at the position ofGenerates a small random disturbance in the vicinity of (a) to generate a feasible point x.

Step 5Metropolis criterion: the upper model is user utility maximization, which is different from the standard minimization problem. If it isMake->Let k be ₁ ：＝k ₁ +1, then go to step 3. Otherwise, ifGreater than [0, 1]]Inner one random value, let ∈ ->Let k be ₁ ：＝k ₁ +1, then go to step 3; otherwise, let->Let k be ₁ ：＝k ₁ +1, then go to step 3.

Step 6 cooling schedule: calculating a target valueStandard deviation of (c) with sigma (T) ^(k) ) And (3) representing. The cooling temperature was set as follows:

and let k: the number of times of =k+1,k ₁ =0 then goes to step 2.

Program design method

Traffic network design method based on user utility maximization #)

Elastic traffic demand (connected logic): generating an equalized OD travel time

# step 1: initializing. Data and necessary packets are entered in a format.

#1.1 load the compute shortest path packet, prepare to invoke dijkstra shortest path algorithm, note that iggraph packet needs to be installed for the first time, and then can invoke.

# install packages ("iggraph") the first time an iggraph package needs to be installed at run-time

library(igraph)

options(digits＝3)

#1.2 creates an information matrix of the network, comprising all segments. The first column is Road section label (Road), the second column is Road section origin label (Road origin), the third column is Road section destination label (Road destination), the fourth column is free flow time (free flow time), the fifth column is Road traffic capacity (capacity), and the sixth column is Road length (length). For example, the Nguyen-Dupuis network commonly used in traffic research, reference is made here to the program document for detailed parameter settings.

# can also be replicated in Excel and then executed

#e＝read.delim(″clipboard″，header＝F)

e＝matrix(c(1，1，5，7.0，800，4.00，2，1，12，9.0，400，6.00，3，4，5，9.0，200，5.00，4，4，9，12.0，800，8.00，5，5，6，3.0，350，2.00，6，5，9，9.0，400，5.00，7，6，7，5.0，800，3.00，8，6，10，13.0，250，8.00，9，7，8，5.0，250，3.00，10，7，11，9.0，300，6.00，11，8，2，9.0，550，5.00，12，9，10，10.0，550，6.00，13，9，13，9.0，600，5.00，14，10，11，6.0，700，4.00，15，11，2，9.0，500，6.00，16，11，3，8.0，300，5.00，17，12，6，7.0，200，4.00，18，12，8，14.0，400，6.00，19，13，3，11.0，600，7.00)，ncol＝6，byrow＝TRUE)

colnames(e)＝c(″Road″，″Road origin″，″Road destination″，″Free Time″，″Road capacity″，″Road length″)

# e# used to check program breakpoints

#1.3 the traffic demand matrix d0 is entered, the first column is the sign (OD pair) of the origin-destination pair, the second column is the origin sign (origin), the third column is the destination sign (destination), and the fourth column is the traffic demand (demand).

sg=c (2000, 3000) # potential traffic occurrence for each cell, note that there are only two origins in the Nguyen-Dupuis network.

d0 =matrix (c (1, 2,0.5×sg [1],2,1,3,0.5×sg [1],3,4,2,0.5×sg [2],4, 3,0.5×sg [2 ]), ncol=4, bypass=true) # initial traffic distribution scheme, where the average distribution is used

colnames(d0)＝c(″OD pair″，″Origin″，″Destination″，″Demand″)

#d0# is used to check program breakpoints

Traffic distribution function: the function of the Frank-Wolfe algorithm for traffic distribution is customized, note that the input demand matrix d is in the form of d0, the traffic network e is in the form of e, and the relative error (RRSE) of iteration stop is 0.01.

fw＝function(e，d)

{

#1.4 computing the shortest Path and Path traffic for each OD pair based on Path free flow time

g=add.edges (graph.empty (13), t (e [,2:3 ]), weight=e [,4 ])# creates a graph, 13 is the number of nodes, and time is taken as weight instead of path length

b12 Get short paths (g, from= "1", to= "2", mode= "out", output= "epath") and $ epath [ [1] ] # shortest path from start 1 to end 2

b13 Get short paths (g, from= "1", to= "3", mode= "out", output= "epath") and $ epath [ [1] ] # shortest path from start 1 to end 3

b42 Get short paths (g, from= "4", to= "2", mode= "out", output= "epath") and $ epath [ [1] ] # shortest path from start 4 to end 2

b43 Get short paths (g, from= "4", to= "3", mode= "out", output= "epath") and $ epath [ [1] ] # shortest path from start 4 to end 3

# create a matrix for preserving shortest paths and traffic for each OD pair

V＝cbind(e[，1])

st0=numeric (4) # stores the shortest running time of each OD versus free stream

colnames(V)＝″Road″

V

Shortest path and traffic for the # OD pair 12

sp12 = as.vector (b 12) # translates to Road segment label (Road)

sum of time of each section st0[1] =sum (e [ sp12,4 ])#

x12=cbind (e [ sp12,1], rep (d [1,4], length (sp 12))) # section label and traffic, iteration start in algorithm

colnames(x12)＝c(″Road″，″V12″)

x12

V=merge (V, x12, by= "Road", all=true) # defines V as a matrix that exclusively holds the start of an iteration

V[is.na(V)]＝0

V

Shortest path and traffic for # OD pair 13

sp13 = as.vector (b 13) # translates to Road segment label (Road)

sum of time of each section st0[2] =sum (e [ sp13,4 ])#

x13=cbind (e [ sp13,1], rep (d [2,4], length (sp 13))) # section label and traffic, iteration start in algorithm

colnames(x13)＝c(″Road″，″V13″)

x13

V=merge (V, x13, by= "Road", all=true) # defines V as a matrix that exclusively holds the start of an iteration

V[is.na(V)]＝0

V

Shortest path and traffic for the # OD pair 42

sp42 = as.vector (b 42) # translates to Road segment label (Road)

sum of time of each section st0[3] =sum (e [ sp42,4 ])#

x42=cbind (e [ sp42,1], rep (d [3,4], length (sp 42))) # section label and traffic, iteration start in algorithm

colnames(x42)＝c(″Road″，″V42″)

x42

V=merge (V, x42, by= "Road", all=true) # defines V as a matrix that exclusively holds the start of an iteration

V[is.na(V)]＝0

V

Shortest path and traffic for # OD pair 43

sp43 = as.vector (b 43) # translates to Road segment label (Road)

sum of time of each section st0[4] =sum (e [ sp43,4 ])#

x43=cbind (e [ sp43,1], rep (d [4,4], length (sp 43))) # section label and traffic, iteration start in algorithm

colnames(x43)＝c(″Road″，″V43″)

x43

V=merge (V, x43, by= "Road", all=true) # defines V as a matrix that exclusively holds the start of an iteration

V[is.na(V)]＝0

V

# when the flows on all shortest paths are summed, the initial flow is obtained

VS＝rowSums(V[，seq(ncol(V)-3，ncol(V))])

VS

# step 2: updating impedance of each road segment

t0=e, 4] # free stream time

c=e, 5- # road traffic capacity

a＝0.15

b＝4

tp＝function(v){

t0*(1+a*(v/c)^b)

}

repeat{

# step 3: searching for the next iteration direction

g2 The image is constructed by =add.edges (graph.empty (13), t (e [,2:3 ]), weight=tp (VS)) # and 13 is the number of nodes, and the road section impedance update function defined in step 2 is called

b12 Get short paths (g 2, from= "1", to= "2", mode= "out", output= "epath") and $ epath [ [1] ] # shortest path from start 1 to end 2

b13 Get short paths (g 2, from= "1", to= "3", mode= "out", output= "epath") and $ epath [ [1] ] # shortest path from start 1 to end 3

b42 Get short paths (g 2, from= "4", to= "2", mode= "out", output= "epath") for the shortest path $epath [ [1] ] # from start point 4 to end point 2

b43 Get short paths (g 2, from= "4", to= "3", mode= "out", output= "epath") and $ epath [ [1] ] # shortest path from start 4 to end 3

# create a temporary matrix for preserving shortest paths and traffic for each OD pair

V＝cbind(e[，1])

colnames(V)＝″Road″

V

Shortest path and traffic for the # OD pair 12

sp12 = as.vector (b 12) # translates to Road segment label (Road)

colnames(x12)＝c(″Road″，″V12″)

x12

V[is.na(V)]＝0

V

Shortest path and traffic for # OD pair 13

sp13 = as.vector (b 13) # translates to Road segment label (Road)

colnames(x13)＝c(″Road″，″V13″)

x13

V[is.na(V)]＝0

V

Shortest path and traffic for the # OD pair 42

sp42 = as.vector (b 42) # translates to Road segment label (Road)

colnames(x42)＝c(″Road″，″V42″)

x42

V[is.na(V)]＝0

V

Shortest path and traffic for # OD pair 43

sp43 = as.vector (b 43) # translates to Road segment label (Road)

x43=cbind (e [ sp43,1], rep (d [4,4],1ength (sp 43))) # section label and traffic, iteration start in algorithm

colnames(x43)＝c(″Road″，″V43″)

x43

V[is.na(V)]＝0

V

# when the flows on all shortest paths are summed, the iteration direction is obtained

VS2＝rowSums(V[，seq(ncol(V)-3，ncol(V))])

VS2

# step 4: and calculating an iteration step. And calling an optimize function to perform one-dimensional search.

step＝function(1amda){

x2＝VS2

x1＝VS

q＝x1+lamda*(x2-x1)

sum((x2-x1)*tp(q))

}

# lambda=uniroot (step, c (0, 1)) $ root# notes that the range of values of lambda cannot be too long for the step size, the sign of the function values at both ends of the uniroot requirement c (0, 1) is opposite, and some functions are not necessarily satisfied. The use of the optimize function ensures that the optimal value of the unary function is found.

g＝function(lamda){step(lamda)^2}

lamda＝optimize(g，c(0，1))$minimum

lamda

# step 5: determining a new iteration start

VS3＝VS+lamda*(VS2-VS)

VS3

# step 6: convergence test

if ((sqrt (sum ((VS 3-VS)/(2)))/sum (VS)) < 0.01) break# satisfies the convergence condition and exits

Vs=vs3# if the convergence condition is not met then replacing the origin VS with the new point VS3, and cycling through until convergence

}

# step 7: and outputting a result matrix result and an OD driving time matrix u of the balanced state.

# step 7.1: and outputting the flow, the transit time and the speed of each path in the equilibrium state.

result＝cbind(e[，1]，round(VS，0)，tp(VS)，e[，6]/(tp(VS)/60)，e[，5]，round(VS，0)/e[，5])

colnames(result)＝c(″Road″，″Volume″，″Time″，″Speed″，″Road Capacity″，″Level of Service″)

# step 7.2: outputting each OD driving time matrix u

g=add.edges (graph.empty (13), t (e [,2:3 ]), weight=result [,3 ])# creates a graph, 13 is the number of nodes, and result is the matrix generated in step 7

Creating a travel time matrix for storing the travel time of each OD pair, initially assuming each OD travel time to be 0

u＝matrix(c(1，1，2，0，2，1，3，0，3，4，2，0，4，4，3，0)，ncol＝4，byrow＝TRUE)

Travel time of # OD pair 12

sp12 = as.vector (b 12) # translates to Road segment label (Road)

u [1,4] = sum (result [ sp12,3 ])# summation of each section time

Travel time of # OD pair 13

sp13 = as.vector (b 13) # translates to Road segment label (Road)

u [2,4] = sum (result [ sp13,3 ])# summation of each section time

Travel time of # OD pair 42

sp42 = as.vector (b 42) # translates to Road segment label (Road)

u [3,4] = sum (result [ sp42,3 ])# sum of each section time

Travel time of # OD pair 43

sp43 = as.vector (b 43) # translates to Road segment label (Road)

u [4,4] = sum (result [ sp43,3 ])# summation of each section time

u=cbind (u, st 0) #od inter-pair free-running time st0

colnames(u)＝c(″Road″，″Origin″，″Destination″，″Shortest travel time″，″Free travel Time″)

# output result matrix and OD travel time matrix in list form

list(result，u)

}

# fw (e, d 0) # for checking program break points

# step 8: traffic distribution function: the custom destination selected polynomial Logit function mLogit is input as each OD driving time matrix u, traffic demand sg of each place, existing employment number ee, increased employment number ie and output as new traffic distribution matrix.

nk = 0.8# used Logit parameters

sg1＝NULL

mlogit＝function(u，sg，ee，ie)

{

d＝numeric(4)

L1=log (exp ((-0.1×u [1,4] +log (ee [1] +ie [1 ]))/nk) +exp ((1-0.1×u [2,4] +log (ee [2] +ie [2 ]))/nk)) # log sum of cell 1

L4=log (exp ((-0.1×u [3,4] +log (ee [1] +ie [1 ]))/nk) +exp ((1-0.1×u [4,4] +log (ee [2] +ie [2 ]))/nk)) # log sum of cell 4

sg1[1]＝sg[1]*exp(-1+nk*L1)/(1+exp(-1+nk*L1))

sg1[2]＝sg[2]*exp(-1+nk*L4)/(1+exp(-1+nk*L4))

d[1]＝sg1[1]*exp((-0.1*u[1，4]+log(ee[1]+ie[1]))/nk)/(exp((-0.1*u[1，4]+log(ee[1]+ie[1]))/nk)+exp((1-0.1*u[2，4]+log(ee[2]+ie[2]))/nk))

d[2]＝sg1[1]*exp((1-0.1*u[2，4]+log(ee[2]+ie[2]))/nk)/(exp((-0.1*u[1，4]+log(ee[1]+ie[1]))/nk)+exp((1-0.1*u[2，4]+log(ee[2]+ie[2]))/nk))

d[3]＝sg1[2]*exp((-0.1*u[3，4]+log(ee[1]+ie[1]))/nk)/(exp((-0.1*u[3，4]+log(ee[1]+ie[1]))/nk)+exp((1-0.1*u[4，4]+log(ee[2]+ie[2]))/nk))

d[4]＝sg1[2]*exp((1-0.1*u[4，4]+log(ee[2]+ie[2]))/nk)/(exp((-0.1*u[3，4]+log(ee[1]+ie[1]))/nk)+exp((1-0.1*u[4，4]+log(ee[2]+ie[2]))/nk))

cbind(u[，1：3]，d)

}

# step 9: defining a traffic network e, elastic traffic demands d0 and sg, and alternately iterating and balancing functions of traffic distribution and traffic distribution under the destination employment number ee and the increment ie. For the initial traffic distribution, the travel time matrix of each OD in the user equilibrium state can be obtained, the user reselects the destination according to the matrix, and for the new traffic distribution, a new travel time matrix can be generated, and the process is circulated until the traffic distribution matrix is not changed any more.

cda＝function(e，d0，sg，ee，ie){

k＝0

repeat{

d1＝mlogit(fw(e，d0)[[2]]，sg，ee，ie)

k＝k+1

if (sqrt (sum ((d 1, 4-d 0, 4)) 2))/sum (d 0, 4) less than 0.01) break# stops when meeting a certain precision requirement, and the relative error is 0.01

d0, 4=d0, 4+ (1/k) ×d1, 4-d0, 4 # an iterative weighting method (Method of Successive Average, MSA) is used here, where the inverse of the number of loops is used as the weight, decreasing with increasing number of loops.

# print (d 1) # for checking break points of program

# print (k) # for checking program

if (k= =100) break# jumps out of the loop if the number of loops reaches 100 but the accuracy requirement has not been met

}

#print(d1)

d2＝fw(e，d1)

u＝d2[[2]]

sw1＝-1+nk*L1

sw4＝-1+nk*L4

d4＝cbind(d1，as.vector(c(sw1，sw1，sw4，sw4)))

#print(d2)

list (d 4, d 2) # outputs traffic distribution d1 and traffic network representation d2

}

The next section is for the inspection procedure, i.e. the current system performance without increasing employment.

# ee=c (5, 5) # current employment population

# ie=c (0, 0) # no increase employment

Selection of allocation scheme with minimum system travel time for #d3=cda (e, d0, sg, ee, ie)

The # -sg [1] (ee+ie)% ] exp (-0.01 x d3[ [2] ] [ [2] ] [1:2,4 ])+ sg [2] (ee+ie)% ] exp (-0.01 x d3[ [2] ] [ [2] ] [3:4,4 ])# outputs the network reachability at this time.

ee=c (50, 50) # current employment population

ie=c (0, 0) # increase employment, if no employment is increased, ie=c (0, 0)

swm＝cda(e，d0，sg，ee，ie)

sum(swm[[1]][1：2，4])#Origin 1

sum(swm[[1]][3：4，4])#Origin 4

social total effect of sw0= swm [1] +, 4]% swm [1] +, 5# primary state

Step # 10: simuated annealing algorithm

Installation of the packages generating Dirichlets is carried out by # instruments packages ("DIRECT"), which are first installed and again not needed

library (DIRECT) # load packet

set (12345) # sets a random seed so that the same random number is generated each time

nsd=100# represents the number of dirichlet samples, and for each comparison, the optimal one is found

Number of candidate segments of candidate for candidate=7 #

inv=2000# investment budget

fd=null# stores a viable Dirichlet allocation scheme, with the aim of finding decision variables after determining the optimal performance scheme

st=numeric (nsd) # social utility (social wellfire) for storing each Dirichlet viable solution

The start-up time of the # calculation program-! The following is carried out The foregoing is to define the function and not take up computation time.

timestart＜-Sys.time()

for(iin(1：nsd))

{

rd=rdichlet (1, rep) (1, handle)) # generates a Dirichlet random distribution

len=inv×rd/(0.3×c (5,8,2,3,5,5,7))# increase value of traffic capacity of alternative road section

If len is not a viable allocation scheme, then regenerating an allocation scheme until a viable scheme is generated, 400 is the maximum increase in capacity set

while(max(len)＞400)

{

#seed()

rd=rdichlet (1, rep) (1, handle)) # regenerates a Dirichlet random distribution

}

fd＝rbind(fd，len)

# below is the construction of a new network map e for each allocation model

Ren=cbind (c (3, 4,5,9, 13, 16, 19), as.vector (len)) # an update scheme

colnales (Ren) =c ("Road", "Enhancement") # and the link number are combined together

Ren=merge (e, ren, by= "Road", all=true) # merge

Ren2, 7 is. NA (Ren 2, 7 ]) =0# replaces NA with 0, since NA participates in the calculation as a result of NA

Ren2, 5=ren 2, 5+ren 2, 7# is added to the original capacity to become an improved capacity

e1 =ren2 [,1:6] # a new network e1

d3 =cda (e 1, d0, sg, ee, ie) # calls the cda function defined previously for elastic traffic demand and traffic distribution d0

Total utility of the # solving system

st[i]＝d3[[1]][，4]％*％d3[[1]][，5]

print ("ith," social utility of "feasible solution is" round (st [ i ], 2))) # looks at intermediate results

}

Below # is simulated annealing method

Step 11 #: determining an initial temperature and an initial point

Te=null# notes that T cannot be used, that is, the system reservation name, indicates logical true

Te[2]＝var(st)#temperatrue

Ren=cbind (c (3, 4,5,9, 13, 16, 19), fd [ whish. Max (st), ]) # outer loop, optimal update scheme, initial point

swt =c (sw 0, max (st)) # outer circulation for storing changes in social utility

Ts＝1#stop temperature

The k=2# controls the outer loop, since the first column road number is added, here 1 is added

k1 Control inner loop =1 #

L0=20# length of markov chain, 20 times the number of decision variables

Ad=0.02 # percentage of original investment as field

One parameter used by par=0.5# cooling formula

Step # 12: internal and external circulation

repeat{

if(Te[k]＜Ts)break

Ren_in=matrix (0, nr=handover, nc=l0×handover) # stores the internal cyclic scheme ren_in [,1:2]

social utility of swt _in=null# deposit internal circulation

Ren_in[，1]＝Ren[，k]

swt_in[1]＝swt[k]

for(k1in(2：(L0*candicate)))

{

#k1＝3

repeat

{

rp0＝runif(candicate，min＝-1，max＝1)

rp1＝rp0-mean(rp0)

len=inv+ad+1/(0.3×c (5,8,2,3,5,5,7)) # change value of traffic capacity of alternative road segment

x＝Ren_in[，k1-1]+len

if(min(x)＞0&max(x)＜400)break

}

Ren1＝cbind(Ren[，1]，x)

colnales (Ren 1) =c ("Road", "Enhancement") # and the link number are combined together

Ren=merge (e, ren1, by= "Road", all=true) # merge

e1 =ren2 [,1:6] # optimal network e

# computing social utility for New networks

F1＝d3[[1]][，4]％*％d3[[1]][，5]

#Metroplis Rule

if(F1＞swt_in[k1-1]){

Ren_in[，k1]＝round(Ren1[，2]，2)

swt_in[k1]＝F1

}else if((exp(F1-swt_in[k1-1])/Te[k])＞runif(1)){

Ren_in[，k1]＝round(Ren1[，2]，2)

swt_in[k1]＝F1

}else{

Ren_in[，k1]＝Ren_in[，k1-1]

swt_in[k1]＝swt_in[k1-1]

}

See intermediate results #)

print (paste ("th", k-1, "outer loops", k1, "social utility of an inner loop feasible scheme is" swt _in [ k1 ])

}

#Cooling

sw_sd＝sd(swt_in)

Te[k+1]＝Te[k]*(1+Te[k]*log(1+par)/(3*sw_sd))^(-1)

Ren＝cbind(Ren，Ren_in[，ncol(Ren_in)])

swt[k+1]＝swt_in[length(swt_in)]

k＝k+1

}

# step 13: output the final required result

#Ren1＝Ren

Ren1＝Ren[，c(1，ncol(Ren))]

Ren=merge (e, ren1, by= "Road", all=true) # merge

e1 =ren2 [,1:6] # optimal network e

# computing social utility for New networks

F1＝d3[[1]][，4]％*％d3[[1]][，5]

F1# output social utility

write csv (d 3[ [1] ], file= "traffic distribution and social utility matrix. Csv") # is held to the current working catalog in csv format

write. Csv (d 3[ 2] ] [1] ], file= "road section equalization result. Csv") # is held to the current working directory in csv format

write. Csv (d 3[ [2] ] [ [2] ], file= "OD travel time. Csv")

write. Csv (Te, file= "Cooling. Csv")

getwd () # looks at the save address of the output file

Runtime of # # computing program

timeend＜-Sys.time()

runningtime＜-timeend-timestart

print (runningtime) # output runtime.

The beneficial effects are that: the invention develops a new software for traffic network design by using R language, wherein the upper layer aims to maximize the user utility, and the lower layer is an iterative feedback process between a Nested Logit model and a user balance model, wherein the OD travel time is determined by endogenous generation. Since the bilayer model is generally considered a strong NP-complete problem, conventional optimization algorithms are computationally difficult to implement. Therefore, the invention provides a meta-heuristic simulated annealing algorithm to solve the problem. In addition, the software has the following technical characteristics: 1) Open source free: the R language is open source software, and the traffic planning and management software is mostly foreign software, so that the price is high; 2) Easy to operate: the foreign commercial software has complex general interface and complicated application, is difficult to master by common engineers and planners, and is easy to operate; 3) The matrix algorithm is fast and efficient: the foreign business software is mostly developed by adopting the C language, is relatively low in efficiency when the algorithm is realized, and greatly quickens the calculation process by adopting the matrix language. Research shows that the software can be used for traffic network design with maximized user utility.

Description of the drawings:

FIG. 1 is a diagram of a two-layer model structure

FIG. 2 is a cyclic feedback process of an underlying model

FIG. 3 is a flow chart of a simulated annealing algorithm

FIG. 4 is an Nguyen-Dupuis network for simulation studies

The specific embodiment is as follows:

the Nguyen-Dupis network shown in fig. 4 is widely used for method validation in various traffic studies. Table 1 shows road segment parameters including free-flow travel time, road segment traffic capacity, and road segment length. There are two departure areas 1 and 4 and two destination areas 2 and 3 in the Nguyen-Dupis network. Assume that during peak hours, the potential travel demands for departure areas 1 and 4 are 2000pcu/h and 3000pcu/h, respectively. That is, O ₁ =2000 pcu/h and O ₄ ＝3000pcu/h。

TABLE 1 road segment parameters for Nguyen-Dupuis network

There are many explanatory variables in the new logic model of travel-destination combination selection. However, as a simulation study we consider only the most important variable, namely OD travel time t _rs And employee number em of destination s _s . For example, the system utility Y of a representative individual residing at the departure 1 ₁₂ And Y ₁₃ Is assumed to be:

the parameters may be given by a empirical calculation. Typically, the coefficient of OD travel time is negative, as it has a negative effect. Assume that the staff numbers are em respectively ₂ =5ten thousand em ₃ =5ten thousand. The expected travel utility can be reduced to:

V ₁ ＝-1+λln(exp(Y ₁₂ /λ)+exp(Y ₁₃ /λ)) (19)

wherein the utility left in the home (not going out) is 0, i.e., the reference value is 0.W (W) ₁ Expressed as a constant-1, i.e., the outbound intrinsic utility is negative. Let lambda be 0.8.

In trip distribution, the study uses a classical user equalization method that combines the road segment performance function with the equalization state. One common road segment performance function given by the united states highway office (BPR) is as follows:

wherein t is _a (v _a ，c _a ) Is given a traffic volume v _a The traffic capacity of road section is c _a Is a function of the impedance over the segment a of (a),is the free flow impedance of segment a. Alpha and beta are flow/delay coefficients that can be calibrated by demonstration. The BPR is given empirical values of α and β of 0.15 and 4.0, respectively, which are also used in the simulation studies herein.

The convergence criterion for MSA is 0.01, i.e. epsilon=0.01. For a given road network design, a single stable solution can be converged, and the OD travel time matrix of the solution is consistent. Parameters used in simulated annealing algorithms typically affect the calculation time and the calculation result. They are set to δ=0.5, l ₀ ＝20，M ₀ =100 and T _s =1. Assuming a road section investment function of

g _a (Δc _a )＝0.3×Δc _a ×l _a ，a∈A (21)

Wherein l _a Is the length of road segment a (in kilometers). Furthermore, the investment budget is determined as b=2000. The program is written in the open source language R. The final OD distribution matrix is shown in table 2. The results show that travel requirements at departure 1 and departure 4 are only 50.55% and 49.13% respectively. That is, almost halfIs limited. Therefore, we need more investment to improve the user experience. Note that in the routine, the vehicle loading coefficient μ is set to 1, and the marginal utility of the income θ is also set to 1. The optimum network design and road segment operating conditions are shown in table 3. This means that in order to maximize the user utility, road segments 3, 9, 16 need to be serviced while road segments 4,5, 13 need to have increased road traffic capacity.

TABLE 2 travel distribution matrix

TABLE 3 optimal traffic network design and road segment run status

/>

Claims

1. A traffic network design method based on user utility maximization is characterized in that:

(1) Representing the traffic network design problem as a double-layer planning problem with a leader-follower decision structure, wherein the upper layer problem is the utility maximization of network users, and the lower layer problem is the traffic system balancing problem;

(2) The upper layer problem is that the utility of network users is maximized by adjusting the traffic network design, the measure of the utility of the users is based on a Nested logic model, wherein the Nested logic model is adopted for travel generation and destination selection, and the utility of the network users is measured through a utility function of the Nested logic model;

(3) The selection probability of the Nested Logit model is written as the product of two Logit probabilities, the first layer is two choices of whether to go out or not, the second layer is to make multiple choices in the destination of going out, and the utility obtained by the decision maker n living in the departure place r from the destination choice s is expressed as:

wherein the method comprises the steps ofFor the variables of whether to go out, ->For describing the variables of destination s +.>Is an extremum distribution which is independently and uniformly distributed and is toThe expression is as follows:

wherein beta is ₀ Is a constant term of the number of the items,is of coefficient beta ₁ The shortest travel time between origin-destination OD pair rs, +.>Is a vector of observable variables having coefficients β;

wherein the method comprises the steps ofIs the marginal probability of going out, +.>The conditional probability of selecting the destination s after the travel is determined, and the marginal conditional probability and the conditional probability respectively adopt the forms of a binomial Logit and a polynomial Logit, which are respectively expressed as:

wherein the method comprises the steps of

S _r Is the set of destinations of the origin r,is the expected maximum utility of decision maker n from destination selection, +.>Since its expression is also called logsum term, λ is called logsum coefficient, λ reflects the degree of independence between destinations, and a larger λ shows greater independence and smaller correlation, and a range of λ is [0,1]When λ=1 means that the choice between destinations is completely independent, in which case the new logic model is reduced to a standard logic model;

assuming the overall behavioural mechanisms are the same, individual n is typically negligibleSubscript n, potential travel demand O given origin r _r The travel demand at the departure point r is expressed as:

Q _r ＝O _r p _rT ， (7)

wherein p is _rT From equation (4), the travel requirement between origin-destination OD pair rs is defined as:

wherein p is _rs Obtained from equation (3) at which the OD travel distribution matrixThe matrix is measured by people, cannot be directly distributed to the road network, and needs to be converted into a unit of a standard passenger car pcu by using a vehicle loading coefficient mu, so that a travel distribution matrix of the unit of the standard passenger car is expressed as:

q _rs ＝O _r p _rs /μ， (9)

wherein μ is determined by investigation;

the policy goal is to increase the utility of the outgoing journey, and according to binomial selection formula (4) representing whether to travel, the expected travel utility of the individual n residing at the departure place r is expressed as:

since the absolute level of utility is not measurable, it is divided by the marginal utility θ of revenue _n Thereby obtaining the outgoing utility of the user n residing at the departure place r:

equation (11) converts units of utility into currency, assuming residence in traffic minutesThe individual of the analysis region r has the same utility as the representative individual n, and the subscript n is ignored, and the sum CW of the consumer utility of the resident of the living region r is calculated _r Equal toAnd the total travel demand Q of the district _r The product of:

the total utility of the system is then the set of all departure areas, expressed as:

combined solving travel demand Q _r The mathematical formulation of the upper layer is:

wherein B is the investment budget, g _a Is a cost function of the road segment a,the OD travel time is implicitly determined by the lower model, which is the maximum increment allowed on the road section a;

(4) The lower layer problem is traffic system balancing, which is a feedback iteration process between a Nested logic model and a user balancing model, wherein in the Nested logic model, starting from an initial solution of OD travel time, a travel distribution matrix is generated, then travel demands among ODs are distributed to a traffic network by using the user balancing model, so that road section travel time is obtained, then a shortest path algorithm is used for calculating new OD travel time, a continuous average method is adopted, the new OD travel time is expressed as a weighted sum of the OD travel times of the former two times, then the latest OD travel time is input into the Nested logic model, and the cycle iteration is continued until the OD travel time for iteration is equal, and the traffic system balancing state is achieved at the moment;

(5) The double layer planning problem is solved using a simulated annealing algorithm.