阅读说明：本技术 疫情防控导向的高速公路网进口匝道流量控制软件 (Expressway network inlet ramp flow control software for epidemic situation prevention and control guiding ) 是由 林宏志 于璐 于 2021-06-17 设计创作，主要内容包括：本发明涉及一种疫情防控导向的高速公路网进口匝道流量控制软件,属于交通工程技术领域。防疫封锁线是阻断病毒传播的有效手段,然而,如何控制进口匝道的车流量,以使下游的防疫封锁线上各入城通道的最大排队等待时间不超过预定值成为一个亟待解决的技术问题。本发明采用了一个双层规划模型,其上层模型为公路网进口匝道的优化控制,下层模型为交通系统平衡,其中排队等待时间由排队论计算得出。设计了一种启发式算法求解该双层规划模型,其下层为连续平均算法,而上层为带精英策略的遗传算法。采用免费开源的R语言开发了计算软件。具体的实施表明本发明对于管理者确定进口匝道的最佳控制策略十分有效。(The invention relates to epidemic situation prevention and control guiding expressway network entrance ramp flow control software, and belongs to the technical field of traffic engineering. However, how to control the traffic flow of the entrance ramp so that the maximum queuing waiting time of each urban-entering channel on the downstream epidemic prevention blocking line does not exceed a preset value becomes a technical problem to be solved urgently. The invention adopts a double-layer planning model, the upper layer model is the optimization control of the ramp at the entrance of the road network, the lower layer model is the traffic system balance, and the queuing waiting time is calculated by the queuing theory. A heuristic algorithm is designed to solve the double-layer planning model, the lower layer of the model is a continuous average algorithm, and the upper layer of the model is a genetic algorithm with an elite strategy. Computing software was developed in the free-sourced R language. The specific implementation shows that the method is very effective for determining the optimal control strategy of the ramp of the inlet by a manager.)

1. The invention provides a highway network inlet ramp flow control method for epidemic prevention and control guidance, which is characterized by comprising the following technical steps:

(1) a double-layer planning model is adopted, wherein the upper layer is a network manager, and the lower layer is a traffic traveler;

(2) the upper layer model is a nonlinear programming model, the maximum queuing waiting time of each city-entering channel on an anti-epidemic blocking line is taken as a constraint condition, the maximum traffic capacity of a highway network is taken as a target, the flow control of an entrance ramp is taken as a decision variable, and the queuing waiting time of the vehicles on each city-entering channel is calculated by an M/M/c queuing model, wherein a first symbol M represents that the arrival of the vehicles obeys Poisson flow, a second symbol M represents that the service time of a single inspection station obeys negative exponential distribution, and a third symbol c represents the number of the anti-epidemic inspection stations arranged on the single city-entering channel;

(3) the arrival flow of vehicles on each urban entering channel is determined by a lower layer model, the lower layer model is a combined model of travel distribution and traffic distribution, the travel distribution and the traffic distribution achieve traffic system balance through circulation feedback, the travel distribution adopts a multinominal Logit model, the traffic distribution adopts a user balance model with queuing, the queuing waiting time of travelers at an epidemic prevention inspection station is regarded as a part of generalized travel cost, and the travelers can minimize the generalized travel cost of individuals through destination selection and path selection;

(4) a heuristic algorithm is designed to solve, wherein a continuous average algorithm is adopted for a lower model, and a genetic algorithm with an elite strategy is adopted for an upper model.

(5) And writing a program for realizing the calculation by adopting an R language.

The technical field is as follows:

the invention provides epidemic prevention and control guiding expressway network entrance ramp flow control software, and belongs to the technical field of traffic engineering.

Background art:

an epidemic prevention lockout line is a restriction on people who enter and exit a designated geographical area (e.g., a community, city or region), and only people who are qualified through body temperature detection can pass through the lockout line. This has proven to be an effective way to prevent the spread of infectious viruses to other cities. However, queuing at security checkpoints on the lockout line is reported to be too long and waiting costs are too high. Considering that the traffic flow of all highways including the exit ramp is determined by the upstream ramp controller, the invention designs a method and an algorithm for optimizing the control scheme of the entrance ramp of the urban highway network, so that the queuing waiting time of the exit road sections surrounding the epidemic prevention blocking line is in a tolerable range.

There are many ways to measure the performance of a blocking inline queuing system. The most important of these is the queuing delay cost, i.e., vehicle waiting time, which is also an important indicator in traffic network analysis. There are generally two ways to determine queuing delay costs in a traffic network. Yang and Yagar (1994) propose a method for implicitly determining queuing times. They established a two-level planning model for the traffic control problem of the entrance ramp, wherein the upper layer is to determine the flow of the entrance ramp for minimizing the total travel time of the system, and the lower layer is a traffic balance model for describing the queuing and waiting phenomenon at the ramp. It can be demonstrated that the segment queuing time corresponds exactly to the lagrangian multiplier associated with the segment capacity constraint. Yang and Yagar (1995) further extended the two-layer model to optimize signal timing in congested road networks, a study that adequately accounts for queuing and congestion on saturated road segments. In addition to traffic control, Yang and Lam (1996) develop a two-level model for determining the toll modes of bottleneck road sections such as tunnels and bridges. The underlying problem is a network equalization model with queuing that describes the path selection behavior of users in queuing and congestion. The upper layer problem is to determine road tolls for operators, so that the operators can optimize the performance of the system after fully considering the routing behavior of users. Previous studies were conducted with fixed origin to destination travel requirements, and Yang and Bell (1997) extended them further to the network equalization model with queuing with elastic requirements. In this case, the representation of the queuing delay time is not explicit and is implicitly determined as the lagrangian multiplier associated with the segment capacity constraint in the underlying problem. Only if this road capacity is reached, queues will form, otherwise the road travel time will only depend on the flow. Another method of determining queuing delay is the bottleneck model of Vickrey (1969), which uses deterministic queuing theory to explicitly derive queuing delay times, aiming to solve the issue of departure time selection for commuters on roads that are bottleneck-constrained during early rush hours. It can model the formation and dissipation of queues of upstream road segments at bottlenecks in a straightforward and easy-to-handle manner and thus can serve as a baseline model for dynamic representations of traffic congestion during rush hour periods. The classical bottleneck model indicates that the travel time on each road segment is the sum of the free stream time and the queuing delay time before the downstream bottleneck. Over the last 50 years, significant progress has been made in the study of bottleneck models. Through the bottleneck model, scholars have more and more knowledge about the characteristic of rush hour traffic queuing. Li et al (2020) have recently been well summarized on bottleneck models and have shown promise in their future research directions where more details can be understood by interested readers.

Although the queuing delay at the lockline is similar to that at turnstile controllers, signalized intersections, and bottlenecks, it has different characteristics, which provides potential opportunities for further innovation. First, the detection time for each vehicle is random, since this time depends on the number of uncertain passengers in the vehicle. Thus, conventional deterministic queuing models may result in large deviations of the results from the actual values, thereby limiting the practical application of the model. The invention particularly adopts a random queuing model, and can model the uncertain detection time of each inspection station. Second, the locklines emphasize the desired level of queuing delay. Although queuing has been well studied in traffic networks, no method has been studied for achieving a desired level of queuing delay in quarantine lines. Finally, existing research has focused primarily on travel behavior analysis and demand-side strategies (especially congestion pricing), with little focus on supply-side strategies (Li et al, 2020). In summary, the present invention is directed to an application software for optimizing the flow control of the ramp entrance of the urban expressway network for epidemic prevention and control.

The invention content is as follows:

the technical problem is as follows: the invention aims to solve the technical problem of how to control the flow of the entrance ramp of the expressway network, so that the maximum queuing waiting time of downstream urban entering channels does not exceed a preset value, and the traffic capacity of the expressway network is maximized, thereby meeting the requirement of epidemic situation prevention and control and maximizing the transportation capacity of highways.

The technical scheme is as follows: road traffic is indirectly controlled at the entrance ramp, which is a Stackelberg game with a leader-follower decision structure. The operator at the upper level controls the entrance ramp to maximize the traffic. The operator can predict but not control the travel behaviour of highway users including destination selection and route selection, and all users make their own travel decisions in a user-optimal way. In addition, the decision of the lower layer user is made after the upper layer decision. However, the operator must adjust his decision in anticipation of the user's behavioral response. The invention proposes a two-layer optimization model to describe the size of the leader and follower of the relationship between the operator and the user, and the conceptual framework is shown in fig. 1. The upper layer is the maximization of network traffic with queuing delay constraint, and the lower layer is the balance of traffic system with queuing, which can describe the predicted travel behavior of highway users. Note that balancing is a feedback process between travel distribution and traffic distribution. The invention comprises the following steps:

traffic system balance model with queuing

The lower model is a traffic system balance combining travel distribution and traffic distribution models. It has long been stated that travel times are inconsistent in conventional four-phase sequential models because travel times are actually determined endogenously. Generally, two approaches can be taken according to the literature to solve this problem of inconsistency to achieve a balance of the traffic system. One approach is to combine several steps into an equivalent mathematical plan that can prove good convergence and consistency of results (Oppenheim, 1995; Sheffi, 1985). The other is to iteratively feed back sequential models until the travel time meets the consistency criterion (Boyce and Zhang, 1997; Boyce et al, 1994). While the former is commonly used in the literature, the latter is more flexible at each step (Lin and Wei, 2019; Lin, 2019). Therefore, the second approach, sequential model with feedback, is used here.

Note that traffic distribution here is not a traditional distribution because it involves queuing delays on the quarantine line. The determination of queuing delay time is a critical issue, and queuing theory is generally an excellent tool for analyzing queuing waiting costs. In most traffic situations, the time interval and service time of the two inputs is randomly described by an exponential distribution. Based on this, the present invention here employs a stochastic queuing model based on the Poisson distribution assumption, i.e. the inter-arrival and service time follow an exponential distribution. The derivation of the specific queuing model is based on the steady state of the queuing situation, and is implemented after the system runs for a long enough time.

According to conventional traffic flow theory (Gartner et al, 1999), the queuing phenomenon at each toll station can be described by an M/c queuing model, where M represents the distribution of markov (or poisson) arrivals or departures, or equivalently exponential arrival or service time distribution, c represents the number of parallel service desks with equal service rates, and there may be one or more parallel inspection desks (i.e., service desks) on each entry leg. Assuming that there are m entrances on one quarantine line, the vehicle is predicted to flow in lambda according to the poisson process_iTo each inlet i (i ═ 1, 2.., m), and c of each inlet i_iSeveral parallel checkpoints have the same parameter, mu, of exponentially distributed service time, where c_iμ＞λ_i。

Similar to a single service facility, the most common measure of queuing conditions at a given import i is the expected number of vehicles in line (l)_i) And expected queuing delay time (d)_i)。l_iAnd d_iThe relationship between them is called Little formula, the formula is l_i＝λ_id_iThis relationship is valid under fairly common conditions. When rho_i＝λ_iWhen/. mu.h, expression l_iThe following can be determined:

wherein p is_iIs the steady state probability of no client at import i, A^*For the set of entry segments, equation (3) is the steady state condition. Queue waiting time d_iCan pass through l according to Little formula_iDivided by λ_iAnd (4) determining. The concrete formula is as follows:

the travel time of the entrance section consists of two parts. One part is link travel time determined by link traffic flow and the other part is queue delay time determined by traffic flow and the number of parallel checkpoints. Note that assuming that the physical length of the queue is zero and there is no queue back overflow, this means that the road segment travel time is independent of the length of the queue, the travel time required at the road segment being determined only by its traffic. I.e. when a ∈ A^*Time, t_a(v_a，c_a)＝t_a(v_a)+d_a(v_a，c_a) When in factWhen t is_a(v_a，c_a)＝t_a(v_a) Wherein t is_a(v_a) Is the time of travel of the road section, d_a(c_a) Is the latency. Traffic allocation with queuing is therefore a conventional problem for a given checkpoint deployment scenario.

Overall, the underlying model is detailed in fig. 2 for a given travel demand and road network. First, a travel distribution matrix is generated by destination selection. Where a polynomial logit model is used for destination selection, it is also considered to be the simplest, most practical discrete selection model. After the travel distribution matrix is generated, travel demands are distributed into a road network through user balance with queuing to generate road traffic flow. Note that the predicted traffic on the ingress leg is the average arrival rate of the queuing system. Then, travel times for all origins and destinations (OD pairs), including link travel time and queuing delay time, are calculated by Dijkstra's algorithm. These path travel times are fed back to the polynomial logit model to update the travel distribution matrix. This process is repeated until the travel distribution matrix is no longer changing. This state is achieved as traffic system balance. Fig. 2 illustrates the feedback process of the underlying model. The variable symbols used therein are defined as follows:

q_rs: travel demand between origin r and destination s;

O_r: travel demands within region r;

S_r: a set of destinations of travelers from origin r;

β_s: a traveler's preference for destination s;

t_rs: the route travel time between the origin r and the destination s;

a: a set of all segments in the network;

A^*: is a collection of import road segments;

β_t: travel time t_rsThe coefficient of (a);

v_a: traffic flow on road segment a;

c_a: number of checkpoints on road segment a;

t_a: travel time on road section a, which is traffic flow v_aAnd number of inspection stations c_aA function of (a);

a traffic flow on a route k connecting the origin r and the destination s;

the link and path association is expressed as:

(II) double-layer planning model

For a given upper level on-ramp control scheme, there is a traffic system balance with queuing at the lower level. The goal of the operator is to control the ingress ramp traffic to keep the queue latency at each egress ramp within a tolerable range and maximize network traffic. An acceptable control scheme should take into account queuing delay constraints. Note that the service level for a given off-ramp i is the segment flow λ_iAnd the segment flow may be controlled by an inlet turnstile controller. The average queuing time provides a constraint for determining an acceptable rate of on-ramp control, i.e.

d_i(λ_i，c_i)≤T，

The constant T is a desired level specified by the decision maker, e.g., T ═ 3 min. Note that d_iIs a road section stream lambda_iAnd number of inspection stations c_iA function. The average latency d can be calculated according to equation (4)_iIs more specifically determined as

Yang et al (1994) proposed a two-level optimization method in which the on-ramp traffic control problem is described as the Stackelberg game. Inspired by their research, the present invention also employs a two-tier planning model. The highway operator at the upper level aims to maximize the total flow through all ramps, subject to queuing delays at the checkpoints. The operator can predict, but has no control over, the travel behavior of the highway users, including destination selection and route selection, and all users make decisions in a user-optimal manner. Thus, the formula of the upper model is determined as:

wherein the symbols are defined as follows:

u_r: the inflow through the inlet ramp R ∈ R;

the traffic demand R of the entrance ramp belongs to R;

r: a set of on-ramp, i.e. origin;

i: a set of exit ramps;

u: all incoming vectors (upper layer decision variables).

In the upper model, the objective function (6) is the maximization of the total input flow through all the inlet ramps; constraint (7) means that the waiting time at each exit ramp should not exceed a predetermined acceptable level T, noting the road section traffic λ_i(u) is a function of all ramp inflows; the constraint (8) is a steady state condition; and the constraint (9) indicates that the inflow of each ramp should be non-negative and less than or equal to the corresponding traffic demand. To obtain the road section flow lambda_i(u), the lower layer problem needs to be solved first.

Traffic system balance model algorithm with queuing

To solve the proposed two-layer planning model, it is preferable to solve the model of the lower layer first, because the model of the lower layer is embedded in the model of the upper layer. Given the decision of the upper layer, the lower layer will have a steady traffic flow pattern, given fixed travel requirements and the established road network. It should be noted that the lower layer is the feedback process between the travel distribution and the traffic distribution with queuing. Among other things, a continuous averaging algorithm (MSA) may be used to achieve system balancing. The initial travel distribution matrix may be generated by a polynomial logic model with initialized origin-destination (OD) versus travel time, and then assigning the demand to the road network by Frank-Wolfe algorithm, and the link traffic flow and the link travel time may be generated. In addition, the queuing delay time may also be determined by the predicted traffic flow of the intake section. The generalized travel time for each of the import road segments includes a road segment travel time and a queuing delay time. According to the first principle of Wardrop routing, known as user balancing, traffic is self-arranged in a congested network with the same and lowest cost for all available paths between OD pairs. Thus, the shortest travel times of the OD pairs can then be updated with Dijkstra's algorithm and then fed back to the polynomial logit model to generate a new travel distribution matrix. However, this matrix cannot be directly assigned to the road network, since convergence of direct feedback is usually not possible, and the continuous travel distribution matrix must be averaged. Although there are also some successful applications of constant weights, this usually does not guarantee convergence. The invention here therefore updates the travel distribution matrix with MSA with decreasing weights, which are the inverse of the number of iterations. And further distributing the updated matrix to the road network, and continuing the iterative process until the continuous matrixes are approximately equal. Convergence is usually measured by the square root of the relative error between the demand matrices for successive rows. If a predetermined error value is met, the iteration is terminated and the steady state is referred to as traffic system balancing. The traffic flow results thus produced are then entered into the upper layer model. Fig. 3 is a flow chart of a system equalization algorithm with queuing. The detailed MSA algorithm steps are as follows:

step 1, inputting an entrance ramp control scheme u from an upper layer and a built road network.

Step 2 based on the initial OD vs. travel timeInitializing trip distribution matrixLet n equal to 1 as the number of iterations.

And 3, traffic distribution with queuing. Travel distribution matrixAnd distributing the road networks by using a Frank-Wolfe algorithm. Generating a road section travel flow v_aAnd a link travel time t_a. Note that the queuing delay time may be determined by the flow on the ingress path segment.

Step 4, updating the space between OD and r-s (namely: between OD and r-s) through a Dijkstra algorithm) The shortest travel time.

And 5, distributing the rows. Polynomial logit model for updating travel distribution matrix

Step 6, averaging the travel distribution matrix by using a decreasing weight modeAnd

and 7, convergence identification. Checking convergence of travel distribution matrix using square root of relative error

Where ε is a predetermined error value. If the convergence condition is satisfied, the iteration is terminated and goes to step 9, otherwise to step 8.

Step 8 ofAnd n: n +1 and then go to step 3.

Step 9 output travel distribution matrixAnd road section traffic flow v_aAnd returning to the upper layer model.

Algorithm of (IV) two-layer planning model

The two-layer programming problem is a well-known NP-hard problem that is difficult to solve with classical optimization algorithms. Even if the upper and lower layers are both linear programming, there is a challenge, and even more, the upper layer is a non-linear integer programming model. For example, in the conventional solution method for the optimal boundary line charging problem based on the gradient, the problem with a large scale often cannot be converged due to the existence of a plurality of optimal solutions. This failure has led to the development of a heuristic algorithm that determines the optimal charge level and charge location. Thus, heuristic methods, especially genetic algorithms, have many successful applications in the literature (Liu et al, 2013; Shepherd and Sumalee, 2004; Sumalee, 2004, 2007). While the heuristic algorithm is time consuming and has no way to ensure global optimality, it has proven successful in solving the problem of borderline toll optimization. Therefore, a genetic algorithm based on the elite strategy is used here. Fig. 4 is a flowchart of the algorithm. The detailed genetic algorithm based on the elite strategy comprises the following steps:

step 1 parameter initialization. Setting genetic algorithmsParameters used in (1) include population size M, maximum algebra Gen, and crossover probability p_cProbability of mutation p_mAlgebraic sign gen ═ 1, elite retention probability p_e. Note that the population size depends on the nature of the problem, but typically contains hundreds of viable solutions.

Step 2 randomly generating a feasible initial population. A chromosome is a solution consisting of m genes, each representing the flow through which a single ingress ramp runs. A chromosome is randomly generated and if not feasible, generation continues until feasible. A total of M feasible chromosomes are generated, scattered over the range of all possible solutions.

And 3, selecting operation. The objective function of the upper model is used as an adaptive function to evaluate the performance of all chromosomes in the population. Note that here, in order to maximize the throughput of the network, the probability p is retained_eLabeling elite population, worst fitness p_eThe population is discarded.

And 4, performing cross operation. (1-p) remaining_e) M chromosomes are used for crossover operations, the parent chromosomes are randomly paired, and the probability of crossover is p_c. If the crossover is performed, the position of a random gene is determined for crossover. If the new chromosome is not feasible according to the constraints in the upper model, other gene locations are tried until feasible. These newly generated chromosomes generally have many of the characteristics of their parents.

And 5, performing mutation operation. Determining a random gene position to mutate in a defined domain, the probability of mutating being p_m. If the new chromosome is not feasible, other gene locations are tried until feasible.

Step 6 generates the next generation population. After genetic manipulation, there is still a possibility that (1-p) is available_e) M chromosomes. Adding annotated p_eM elite groups are used for ensuring that the total number of the populations reaches M, so that the best chromosome in the current generation can be continuously kept unchanged until the next generation, and the understanding quality is ensured not to be reduced due to genetic iteration. The representation of the update algebra is gen: gen + 1.

Step 7 ends the determination. If the algebra reaches the maximum value, namely Gen is more than or equal to Gen, the iteration process is terminated and the optimal deployment scheme of the epidemic prevention locking line is output. Otherwise, go to step 3.

(V) Programming

time for timing

# step 1: initializing, inputting data and necessary packets according to formats;

#1.1 loads the packet for calculating the shortest path, prepares to invoke the dijkstra shortest path algorithm, and notes that the igraph packet needs to be installed for the first use and then can be invoked.

Packs ("igraph") # install igraph pack

library(igraph)

options (digits ═ 3) # reserved three decimal places

#1.2 creates a distance matrix for the graph that contains all of the candidate segments. The first column is a Road segment number (Road), the second column is a Road starting point number (Road origin), the third column is a Road end point number (Road destination), the fourth column is the free flow time of the Road segment, the fifth column is Road traffic capacity (capacity), and the sixth column is Road length (length). Here, taking the Nguyen-Dupuis network commonly used in traffic distribution as an example, the detailed parameter settings can refer to the program document.

Table 1 may also be copied in Excel and then executed

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

When the number of dots is equal to three, e ═ matrix (c (1, 1, 5, 7.0, 800, 4.00, 2, 1, 12, 9.0, 800, 6.00, 3, 4, 5, 9.0, 800, 5.00, 4, 4, 9, 12.0, 800, 8.00, 5, 6, 3.0, 800, 2.00, 6, 5, 9, 9.0, 800, 5.00, 7, 6, 7, 5.0, 800, 3.00, 8, 6, 10, 13.0, 800, 8.00, 9, 7, 8, 5.0, 800, 3.00, 10, 7, 11, 9.0, 800, 6.00, 11, 8, 2, 9.0, 800, 5.00, 12, 9, 10, 10.0, 800, 6.00, 13, 9, 13, 9.0, 800, 5.00, 10, 10.0, 10, 6.00, 13, 9, 13, 9.0, 10.00, 10, 6.00, 6, 13, 9, 9.0, 10, 6, 6.00, 6, 13, 7, 6, 15, 6, 15, 7, 6, 10, 6, 7, 6, 10, 7, 13, 10, 6, 13, 6, 7, 10, 7, 9, 10, 9, 10, 9, 10, 7, 9, 10, 9, 0, 7, 9, 0, 10, 9, 10, 9, 10, 9, 0, 9, 10, 9, 0, 9, 10, 9, 0, 9, 0, 9, 10, 9, 10, 9, 0, 10, 9, 10, 9, 10, 9, 10, 9, 10, 9, the impedance of the section is thus sufficiently large to ensure that there is no flow.

e＝matrix

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

# e # for checking program breakpoints

#1.3 inputs the initial traffic demand matrix d0, with the first column being the origin point pair designation (OD pair), the second column being the origin point designation (origin), the third column being the destination point designation (destination), and the fourth column being the traffic demand (demand).

sg-c (1000 ) # expected traffic needs in various regions, origin 1 and origin 4;

d0 ═ matrix (c (1, 1, 2, 0.5 × sg [1], 2, 1, 3, 0.5 × sg [1], 3, 4, 2, 0.5 × sg [2], 4, 4, 3, 0.5 × sg [2]), ncol ═ 4, bynow ═ T) # initial distribution scheme colnamees (d0) ═ c ("paod" original "," Origin "," Destination "," Demand ")

# d0# used for checking program breakpoints

Parameter values used later in the text

M200 # if the arrival rate exceeds the service capability of the queuing system, the queuing time is assumed to be large M

C (11, 15, 16, 19) # set of entrenched channels (entry links)

Traffic demand of each trip place of miud ═ c (1500 ) #

Allowed number of service stations per inbound channel of # cj (9, 4) #

cj ═ c (9, 3, 5, 5) # currently maximum configuration capability of security checkpoints

Average service rate per service desk mu 2# (2 cars per minute)

AT＝2#Acceptable level of service.2minutes

Hereinafter, the parameters used in the link update function are originally defined within the fw function, but are also used in other functions, and are therefore defined as global variables

t0 ═ e [, 4] # free flow time

Traffic capacity of e, 5# road

a＝0.15

b＝4

The # self-defined Frank-Wolfe algorithm function is characterized in that the form of an input demand matrix d is d0, the form of a traffic network e is e, and the relative error is 0.01. Note that for the inbound channel, its time

# is the sum of trip time and queuing time of road section

fw＝function(e，d，s)

{

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

graph is created with g (add)

Short path (g, from 1 ", to 2", mode out ", output ═ ept [ [1] ] # shortest path from start point 1 to end point 2

Short path (g, from 1 ", to 3", mode out ", output ═ ept [ [1] ] # shortest path from start point 1 to end point 3

Short path (g, from ═ 4 ═ to ═ 2 ═ mode ═ out ═ output ═ eat [ [1] ] # shortest path from start point 4 to end point 2

Short path (g, from ═ 4 ═ to ═ 3 ═ mode ═ out ═ output ═ eat [ [1] ] # shortest path from start point 4 to end point 3

# creates a matrix for holding the shortest paths and flows for each OD pair

V＝cbind(e[，1])

st0 number stores initial OD for shortest travel time

colnames(V)＝″Road″

V

Shortest path and traffic of # OD to 12

Vector (b12) # as sp12 converted to Road segment designation (Road)

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

x12 ═ cbind (e [ sp12, 1], rep (d [1, 4], length (sp12))) # road segment label and flow, iteration start in the algorithm

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

x12

V-merge (V, x12, by "Road", all TRUE) # defines V as a matrix that exclusively holds iteration start points

V[is.ha(V)]＝0

V

Shortest path and traffic of # OD pair 13

Vector (b13) # as sp13 converted to Road segment designation (Road)

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

x13 ═ cbind (e [ sp13, 1], rep (d [2, 4], length (sp13))) # road segment label and flow, iteration start in the algorithm

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

x13

V-merge (V, x13, by "Road", all TRUE) # defines V as a matrix that exclusively holds iteration start points

V[is.ha(V)]＝0

V

Shortest path and traffic of # OD pair 42

Vector (b42) # as sp42 converted to Road segment designation (Road)

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

x42 ═ cbind (e [ sp42, 1], rep (d [3, 4], length (sp42))) # road segment label and flow, iteration start in the algorithm

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

x42

V-merge (V, x42, by "Road", all TRUE) # defines V as a matrix that exclusively holds iteration start points

V[is.na(V)]＝0

V

Shortest path and traffic of # OD pair 43

Vector (b43) # as sp43 converted to Road segment designation (Road)

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

x43 ═ cbind (e [ sp43, 1], rep (d [4, 4], length (sp43))) # road segment number and flow, iteration start in the algorithm

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

x43

V-merge (V, x43, by "Road", all TRUE) # defines V as a matrix that exclusively holds iteration start points

V[is.aa(V)]＝0

V

When the flow on all shortest paths is summed, the initial flow is obtained

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

VS

# step 2: updating the impedance of each road segment

The parameters below # are originally defined inside the sub-function, but since other sub-functions are also used, the object cannot be found, so that the parameters are placed outside the sub-function

# is defined as a Global variable

# t0[, 4] # free flow time

Traffic capacity of e, 5# road

#a＝0.15

#b＝4

tp { # s is the number of servers vector for the queuing system

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

Adding queue time to the city entering channel

Average vehicle flow rate for lambda/v [ el ]/60# urban road (from hourly to minute)

rho＝lamda/mu

td ═ numeric (length (lamda)) # defines a variable residence time

td[lamda！＝0]＝line(s，rho)[lamda！＝0]/lamda[lamda！＝0]

td[lamda＝＝0]＝0

# above is the Little formula, the queuing time is calculated by the queuing length, and it should be noted at this time that, if the inflow is 0,

if # is not 0, the program will report an error, so that a value of lambda 0 should be specially processed

td [ rho/s > ═ 1] ═ M #, defines the queuing time, and assigns the time not meeting the stability condition as large M

tf[el]＝tf[el]+td

tf

}

Controlling the cycle number by 1#, forcibly quitting after the cycle number is reached, and preventing endless circulation (oscillation back and forth)

repeat{

# step 3: find the next iteration direction

g2 (add) is graph (13), t (e, 2: 3), weight tp (vs) # structure diagram, 13 is the number of nodes, and the link impedance is updated

Short, path (g2, from ═ 1 ", to ═ 2", mode ═ out ", output ═ epth"), shortest path from start point 1 to end point 2, #

Short path (g2, from 1 ", to 3", mode out ", output ═ overcoat") and shortest path from start point 1 to end point 3, #

b42 ═ get, short paths (g2, from ═ 4 ", to ═ 2", mode ═ out ", output ═ epth"), shortest path from start point 4 to end point 2[ [1] ] #

Short path (g2, from ═ 4 ═ to ═ 3 ═ mode ═ out ═ output ═ oath [ [1] ] # from start point 4 to end point 3

# creates a temporary matrix for holding the shortest paths and flows for each OD pair

V＝cbind(e[，1])

colnames(V)＝″Road″

V

Shortest path and traffic of # OD to 12

Vector (b12) # as sp12 converted to Road segment designation (Road)

x12 ═ cbind (e [ sp12, 1], rep (d [1, 4], length (sp12))) # road segment label and flow, iteration start in the algorithm

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

x12

V-merge (V, x12, by "Road", all TRUE) # defines V as a matrix that exclusively holds iteration start points

V[is.na(V)]＝0

V

Shortest path and traffic of # OD pair 13

Vector (b13) # as sp13 converted to Road segment designation (Road)

x13 ═ cbind (e [ sp13, 1], rep (d [2, 4], length (sp13))) # road segment label and flow, iteration start in the algorithm

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

x13

V-merge (V, x13, by "Road", all TRUE) # defines V as a matrix that exclusively holds iteration start points

V[is.na(V)]＝0

V

Shortest path and traffic of # OD pair 42

Vector (b42) # as sp42 converted to Road segment designation (Road)

x42 ═ cbind (e [ sp42, 1], rep (d [3, 4], length (sp42))) # road segment label and flow, iteration start in the algorithm

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

x42

V-merge (V, x42, by "Road", all TRUE) # defines V as a matrix that exclusively holds iteration start points

V[is.na(V)]＝0

V

Shortest path and traffic of # OD pair 43

Vector (b43) # as sp43 converted to Road segment designation (Road)

x43 ═ cbind (e [ sp43, 1], rep (d [4, 4], length (sp43))) # road segment number and flow, iteration start in the algorithm

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

x43

V-merge (V, x43, by "Road", all TRUE) # defines V as a matrix that exclusively holds iteration start points

V[is.na(V)]＝0

V

When the flow on all the shortest paths is summed, the iteration direction is obtained

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

VS2

# step 4: calculating an iteration step size

step＝function(lamda){

x2＝VS2

x1＝VS

q＝x1+lamda*(x2-x1)

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

}

And # lamda ═ uniroot (step, c (0, 1)) $ root # takes note of the value range of lamda, the step length cannot be too long, uniroot requires that the signs of function values at two ends are opposite, and some functions are not necessarily satisfied. The use of the optimize function ensures that the optimal value of the univariate function is found.

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

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

#print(lamda)

# step 5: determining a new iteration start

VS3＝VS+lamda*(VS2-VS)

VS3

# step 6: convergence test

if((sqrt(sum((VS3-VS)^2))/sum(VS))＜0.01)break

VS-VS 3#, if the convergence condition is not satisfied, replace the original point VS with the new point VS3, and cycle until convergence

k＝k+1

if (k) break # loop reaches 100 times, it is forced to exit

}

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

# step 7.1: and outputting the flow, the transit time and the speed of each path in the balanced 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 running time matrix u

Creating a graph with g (add)

Short path (g, from 1, to 2, mode out, output, exhaust) shortest path from start point 1 to end point 2

Short path (g, from 1 ", to 3", mode out ", output ═ ept [ [1] ] # shortest path from start point 1 to end point 3

Short path (g, from ═ 4 ═ to ═ 2 ═ mode ═ out ═ output ═ eat [ [1] ] # shortest path from start point 4 to end point 2

Short path (g, from ═ 4 ═ to ═ 3 ═ mode ═ out ═ output ═ eat [ [1] ] # shortest path from start point 4 to end point 3

# creates a travel time matrix for holding the travel time for each OD pair, initially assuming each OD travel time is 0

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

Travel time of # OD to 12

Vector (b12) # as sp12 converted to Road segment designation (Road)

Sum of u [1, 4] ═ sum (result [ sp12, 3]) # segment times

Travel time of # OD to 13

Vector (b13) # as sp13 converted to Road segment designation (Road)

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

Travel time of # OD vs. 42

Vector (b42) # as sp42 converted to Road segment designation (Road)

Sum of u [3, 4] ═ sum (result [ sp42, 3]) # each segment time

Travel time of # OD pair 43

Vector (b43) # as sp43 converted to Road segment designation (Road)

Sum of u [4, 4] ═ sum (result [ sp43, 3]) # each segment time

Obstacle-free travel time st0 between u-cbind (u, st0) # OD pairs

colnames(u)＝c(″OD″，″Origin″，″Destination″，″Path time″，″Free time″)

# outputs the result matrix and the OD travel time matrix in the form of a list

list(result，u)

}

#cj＝c(3，5，4，5)

#cj＝rep(1，4)

# fw (e, d0, cj) # used to check the breakpoint of a program

# step 8: and defining a polynomial Logit function mlogit selected by a destination, inputting each OD driving time matrix u and each region traffic demand sg, and outputting a new traffic distribution matrix.

mlogit＝function(u，sg)

{

d＝numeric(4)

d[1]＝sg[1]*exp(0.5-0.1*u[1，4])/(exp(0.5-0.1*u[1，4])+exp(-0.1*u[2，4]))

d[2]＝sg[1]*exp(-0.1*u[2，4])/(exp(0.5-0.1*u[1，4])+exp(-0.1*u[2，4]))

d[3]＝sg[2]*exp(0.5-0.1*u[3，4])/(exp(0.5-0.1*u[3，4])+exp(-0.1*u[4，4]))

d[4]＝sg[2]*exp(-0.1*u[4，4])/(exp(0.5-0.1*u[3，4])+exp(-0.1*u[4，4]))

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

}

# step 9: an integrated traffic distribution and traffic distribution alternating iterative balancing function under the given traffic demands d0 and sg and the traffic network e is defined. And obtaining a driving time matrix of each OD in the user balance state for the initial traffic distribution, reselecting the destination by the user according to the matrix, and generating a new driving time matrix for the new traffic distribution, wherein the process is circulated until the traffic distribution matrix is not changed. Function (e, d0, sg, s) containing Chinese character

k＝2

repeat{

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

#print(d1)

k＝k+1

if (sqrt (sum ((d1[, 4] -d0[, 4]) < 2))/sqrt (sum ((d0[, 4]) < 2)) < 0.05) break # meets a certain relative error requirement, and then the system stops working

d0[, 4] - [ d0[, 4] + (1/k) × (d1[, 4] -d0[, 4]) # here uses an iterative weighting Method (MSA), here using the inverse of the number of cycles as the weight, decreasing as the number of cycles increases.

Print (d0) # used to check the breakpoint of a program

# cat (k) # used for inspection procedure

if the loop number reaches 100 times but the precision requirement is not met, the loop is also jumped out

}

#print(d1)

d2＝fw(e，d1，s)

#print(d2)

Travel time and queuing time of # input and output city channel

v＝d2[[1]][，2]

t0 (1+ a (v/c) b) # link travel time, note that here the parameters t0, a, c, b are defined in the distribution subfunction fw, but since they are defined as local variables, here the parameters t0, a, c, b are defined as local variables

The # subfunction can not be used directly, the system can report the error and can not find the objects, the solution method is redefined here, or the system is required to be defined as the global variable, thus each subfunction

# may be used and is defined herein as a global variable.

Average vehicle flow rate for lambda ═ v [ e1]/60# (from hourly to minute)

rho＝lamda/mu

td＝line(s，rho)/lamda

td [ rho/s > ═ 1] ═ M #, defines the queuing time, and assigns the time not meeting the stability condition as large M

d3＝cbind(el，tf[el]，td，rho/s，s)

colnames(d3)＝c(″Entry link″，″Travel time″，″Delay time″，″Steady condtion″，″checkpoints″)

# output the list

list(d1，d2，d3)

}

# System equilibrium State

#cda(e，d0，sg，cj)

Test data of # cj ═ rep (5, 4) #

#cj＝c(3，5，4，5)

bef cda (e, d0, sg, cj) # initial equilibrium state before-after analysis for checking whether the program is running normally up to this point

bef

The # # # # # # # # # # # # # # # #### # #################################################################################################### # of the double-layer programming is solved

# step 1: initialization of Initialization parameter

The number of chromosomes of each generation of M2-500 # is slow, and M2 is not suitable to be too large

Gen 30# maximum algebra

pc ═ 0.5# crossover probability

probability of variation of pm ═ 0.1#

pe-0.1 # elite ratio

Mark generation 1#

# defines the function for calculating the residence time for each channel, called later, with the input variables being the helpdesk vector s, the service strength vector rho, and the output being the queue length vector lq.

Care is taken to ensure that the input s and rho satisfy the equilibrium condition rho.s-rho/s < 1, otherwise the function results are incorrect.

line is a subfunction for calculating queue length

{

rho.s ═ rho/s # single service desk service strength

probability of 0 for p0 ═ rep (0, length (s)) # initial value and customer count

Initial value of lq ═ rep (0, length (s)) # and average queue length

for(k1 in(1：length(s)))

{

sum1＝0

for(k2 in 0：(s[k1]-1))

{

Factorial function calculation factorial 1 sum1+ rho [ k1] < k 2/[ factor (k2) # factor

}

The probability p0 that the number of customers is 0 in a p0[ k1] (sum1+ rho [ k1] (s [ k1])/(factor (s [ k1]) x (1-rho. s [ k1 ]))))))) (1) # system

Lq [ k1] ═ rho [ k1] < s [ k1] < rho.s [ k1] < p0[ k1 ]/(factual (s [ k1]) < 1-rho.s [ k1]) < 2 >) # the formula for c (s, p) is substituted into the formula for the queuing length Lq

}

return (lq + rho) # calculates the average queue length, not the queueing length

}

# step 2: generating a reactive initial distribution random;

the pop (M2) is used to store the population of each generation

gep (narrow is Gen, ncol is length (fair) +2) # is used to store the optimum chromosome and its expression for each generation, including generation, chromosome and fitness, so the number of rows is ength (lamda) +2

i-1 # counter for marking chromosomes in each generation

repeat

{

# random Generation of one chromosome

The initial value of chr ═ numeric (length (loud) # where chr represents the decision variable

for (j in 1: length (loud)) chr [ j ] ═ sample (seq (1: loud [ j ]), size ═ 1) # generates an initial solution randomly

# checking whether it meets the constraints and if not, regenerating a new chromosome

bef-cda (e, d0, chr, cj) # calculating balance flow rate

lamda-bef [2] ] [1] ] [ c (11, 15, 16, 19), 2]/60# average traffic flow (from hourly to per minute) into urban passageways

rho＝lamda/mu

td ═ numeric (length (lamda)) # defines a variable residence time

td[lamda！＝0]＝line(cj，rho)[lamda！＝0]/lamda[lamda！＝0]

td[lamda＝＝0]＝0

# above is the Little formula, the queuing time is calculated by the queuing length, and it should be noted at this time that, if the inflow is 0,

if # is not 0, the program will report an error, so that a value of lambda 0 should be specially processed

td [ rho/cj > ═ 1] ═ M #, defines the queuing time, and assigns the time not meeting the stability condition as large M

if (max (td) > AT) next # ensures that the equilibrium condition rho.s ═ rho/s < 1 is met; queuing time of any one in-city channel cannot be larger than AT

# print (chr) # output feasible solution for viewing intermediate calculation

pop[i，]＝chr

cat (i, chr) # outputs feasible solutions for viewing intermediate calculation processes

cat ("\ n") # output line feed character

i＝i+1

if(i＞M2)break

}

repeat{

# step 3: selection operation performs selection operation

Where per M2, ncol 2# is used to store fitness function values for each chromosome, and the first column is the chromosome number

for(i in 1：M2)

{

# calculate its fitness function for each chromosome

per [ i, ] ═ c (i, sum (pop [ i, ])) # finds the one with the highest throughput of all inbound channels as the representation of the allocation scheme

print (c (gen, per [ i ]) # is used to monitor the calculation process

}

# discard the worst-performing Pe and mark the best-performing Pe for later use

per2 per (per, 2, trading TRUE), the question is that the larger the total traffic, the better, in large to small rows

per3 ═ per2[ 1: (M2 x (1-pe)) ] # leaving behind the 1-pe crossover and mutation that was made later, leaving the chromosome of the last pe discarded

per4 ═ per2[ 1: (M2 pe) ] # the elite population with the best performance

# keep the best performance of each generation for later drawing

gep [ gen ], c (gen, pop [ per2[1, 1], ], per2[1, 2]) # note that the first column is the generation number, the last column is the fitness, and the middle is the best chromosome of the generation

# step 4: crossover operation

# chromosomes that are crossed and mutated

pop2 ═ pop [ sample (per3[, 1]), ] # sample functions are randomly ordered, and then pairwise

for(i in seq(from＝1，to＝M2*(1-pe)，by＝2))

{

if (runif (1) > pc) next # cross probability, if the random number is larger than pc, the random number is not crossed, and if the random number is smaller than pc, the random number is crossed

Counter with 0 ═ k #, preventing dead cycles

repeat

{

c1 sample (1) # selects a position for single-point crossing

c2 ═ pop2[ i, c1] # c2 is an intermediate variable in the case of two gene exchanges

pop2[i，c1]＝pop2[i+1，c1]

pop2[i+1，c1]＝c2

#print(pop2[i，])

# circulation control, preventing dead circulation

k＝k+1

if(k＞20)break

# requires that two new chromosomes after crossover are a viable solution, otherwise positions continue to be randomly selected for crossover

bef Cda (e, d0, pop2[ i, j, cj) # calculating the equilibrium flow

lamda ═ bef [ [2] ] [ [1] ] [ c (11, 15, 16, 19), 2]/60#, average vehicle flow (from hourly to minute) into the tunnel

rho＝lamda/mu

td ═ numeric (length (lamda)) # defines a variable residence time

td[lamda！＝0]＝line(cj，rho)[lamda！＝0]/lamda[lamda！＝0]

td[lamda＝＝0]＝0

# above is the Little formula, the queuing time is calculated by the queuing length, and it should be noted at this time that, if the inflow is 0,

if # is not 0, the program will report an error, so that a value of lambda 0 should be specially processed

td [ rho/cj > ═ 1] ═ M #, defines the queuing time, and assigns the time not meeting the stability condition as large M

if (max (td) > AT) next # ensures that the equilibrium condition rho.s ═ rho/s < 1 is met; queuing time of any one in-city channel cannot be larger than AT

bef Cda (e, d0, pop2[ i + 1], cj) # calculating the equilibrium flow

lamda ═ bef [ [2] ] [ [1] ] [ c (11, 15, 16, 19), 2]/60#, average vehicle flow (from hourly to minute) into the tunnel

rho＝lamda/mu

td ═ numeric (length (lamda)) # defines a variable residence time

td[lamda！＝0]＝line(cj，rho)[lamda！＝0]/lamda[lamda！＝0]

td[lamda＝＝0]＝0

# above is the Little formula, the queuing time is calculated by the queuing length, and it should be noted at this time that, if the inflow is 0,

if # is not 0, the program will report an error, so that a value of lambda 0 should be specially processed

td [ rho/cj > ═ 1] ═ M #, defines the queuing time, and assigns the time not meeting the stability condition as large M

if(max(td)＜＝AT)

{ print (pop2[ i ]) # View calculation Process

print (pop2[ i + 1]) # View calculation Process

break}

# ensures that the equilibrium condition of rho.s ═ rho/s < 1 is satisfied; queuing time of any one in-city channel cannot be larger than AT

Note that when the comparison is made, & & if the value of the left object is found to be FALSE, then the right object is not calculated, and this calculation method of & & is called short circuit calculation.

}

}

# step 5: mutation operation

for(i in 1：(M2*(1-pe)))

{

if (runif (1) > pm) next # mutation probability, if the random number is larger than pm, the next # mutation probability is not changed, otherwise, the next # mutation probability is not changed

repeat

{

c3 sample (1) # selects a site for single point mutation

Changing the gene value of the position into a random value in a value range

pop2[i，c3]＝sample(miud[c3]，1)

print(pop2[i，])

# Exit until feasible solution

bef Cda (e, d0, pop2[ i, j, cj) # calculating the equilibrium flow

lamda ═ bef [ [2] ] [ [1] ] [ c (11, 15, 16, 19), 2]/60#, average vehicle flow (from hourly to minute) into the tunnel

rho＝lamda/mu

td ═ numeric (length (lamda)) # defines a variable residence time

td[lamda！＝0]＝line(cj，rho)[lamda！＝0]/lamda[lamda！＝0]

td[lamda＝＝0]＝0

# above is the Little formula, the queuing time is calculated by the queuing length, and it should be noted at this time that, if the inflow is 0,

if # is not 0, the program will report an error, so that a value of lambda 0 should be specially processed

td [ rho/cj > ═ 1] ═ M #, defines the queuing time, and assigns the time not meeting the stability condition as large M

if(max(td)＜＝AT)break

}

}

# step 6: the next generation position is generated.

pop3＝rbind(pop2，pop[per4[，1]，])

pop3# returns to start executing iterations

gen＝gen+1

# step 7: termination judgment

if(gen＞Gen)break

}

# step 8: outputting the required calculation result

gep # Performance of each generation, first row is algebra, last row is total traffic throughput, and decision variables are in the middle

Csv (gep, file ═ gene evolution process. csv ") # is kept in csv format to the current working directory

# uses the last line information of gep to output the optimal calculation result

bef cda (e, d0, gep [ Gen, 2: 3], cj) # calculating the balanced flow

bef # shows the calculation results

gep [ Gen, 2: 3] # network maximum traffic capacity

lamda ═ bef [ [2] ] [ [1] ] [ c (11, 15, 16, 19), 2]/60#, average vehicle flow (from hourly to minute) into the tunnel

rho＝lamda/mu

rho/cj

td ═ numeric (length (lamda)) # defines a variable residence time

td[lamda！＝0]＝line(cj，rho)[lamda！＝0]/lamda[lamda！＝0]

td[lamda＝＝0]＝0

# above is the Little formula, the queuing time is calculated by the queuing length, and it should be noted at this time that, if the inflow is 0,

if # is not 0, the program will report an error, so that a value of lambda 0 should be specially processed

td [ rho/cj > ═ 1] ═ M #, defines the queuing time, and assigns the time not meeting the stability condition as large M

td

Run time of # calculation program

timeend＜-Sys.time()

runningtime＜-timeend-timestart

print (running time) # output runtime

The arrangement of the urban epidemic prevention locking line is an effective means for blocking virus propagation, and people going out can pass through the locking line only after detection. The invention provides software for controlling the traffic flow of an entrance ramp of a highway network so as to ensure that the maximum queuing waiting time of each urban entering channel on an epidemic prevention locking line does not exceed a preset value and maximize the traffic capacity of the highway network. The method can effectively avoid overlong queuing waiting time on the epidemic prevention blocking line, and can help a network manager to determine an effective control scheme of the entrance ramp flow.

Description of the drawings:

FIG. 1 is a conceptual framework of a two-tier planning model;

FIG. 2 is an iterative process of the underlying model;

FIG. 3 is a flow chart of the MSA algorithm;

FIG. 4 is a flow chart of a genetic algorithm based on the elite strategy;

FIG. 5 is a Nguyen-Dupuis road network.

The specific implementation mode is as follows:

in order to verify the validity of the proposed models and algorithms, a specific implementation was made. The Nguyen-Dupuis road network as shown in fig. 5 is widely used in traffic engineering to verify various methods. The characteristics of the road segments, including free-stream travel time, road segment capacity and road segment length, are shown in table 1.

TABLE 1 road segment characteristics of Nguyen-Dupuis road network

There are two origins and two destinations in the Nguyen-Dupuis network. The maximum driving demands at the origin 1 and 4 are 1500veh/h and 1500veh/h, respectively, i.e.Andeach origin has a turnstile controller to control the travel demand into the highway network. All existing road segments are marked 1 to 19 and the destination areas are 2 and 3, so the exit road segments to the destination are 11, 15, 16, 19. The problem is to determine the ramp flow control rate for each origin so that the downstream exit turnsUnder the premise that the queuing waiting time is tolerable, the network traffic is maximized.

The parameters used in the underlying model are summarized below. In addition, the polynomial logic model for destination selection is simplified to

Wherein, beta_sIs a traveler preference for destination s, β_tIs a coefficient of OD versus the travel time of the path between r-s. Beta is a_sAnd beta_tCan be calibrated based on actual data, where β is set₂＝0.5，β₃＝0，β_t-0.1. That is, the traveler preference on the destination area 2 is 0.5 and the traveler preference on the destination area 3 is 0, i.e., the traveler habitually prefers the destination 2. The travel time coefficient is-0.1, which means that travel time is of negative utility. In addition, the BPR road section impedance function is adopted to embody the crowding effect in flow distribution, and the formula is as follows:

whereinIs the free-stream travel time on road segment a; α and β are capacity/delay coefficients that can be calibrated by actual data, and conventionally, α is set to 0.15 and β is set to 4. e.g. of the type_aIs the traffic capacity of the section a. Other symbols are consistent with the foregoing definitions. The convergence criterion for MSA is set to be 0.01. For the established road network, a stable traffic system can be realized.

The lower model is nested in the upper model. The parameters used in the upper layer model are listed below. The population scale M is 500; maximum algebra Gen is 30; proportion of elite population p_e0.1; cross probability p_e0.1; probability of variation p_m0.5. Although these parameters are often used in genetic algorithms, they are also valuedThe parameters have to be adjusted to find a reasonable setting for the problem. Assume that the maximum acceptable latency T for each exit ramp is 2 min. The number of test checkpoints for the exit segments 11, 15, 16, 19 is 9, 3, 5, respectively, depending on their physical condition.

Note that traffic volume is typically in hours (veh/h) and arrival rate is typically in minutes (veh/min), so unit conversion is required. The average service rate of a single checkpoint is assumed to be 2veh/min, i.e. the checkpoint examines two vehicles per minute on average. The running time is 2.8h when the computer is programmed by using a popular open source language R3.6.3 in a personal computer provided with an Intel Core i7-4790CPU @3.60 GHz. The maximum network traffic is 2366 veh/h. The on-ramp control rate at origin 1 is 1385veh/h and the control rate at origin 2 is 981 veh/h. At this time, the specific network performance is shown in table 2.

TABLE 2 queuing Performance with maximum network traffic

The results show that the traffic flow is different for different exit road sections. The traffic flow of the outlet road section 11 is maximum and is 1025 veh/h; and the outlet section 15 has the smallest flow, only 319 veh/h. The waiting time for each exit segment is also different, but none of the predetermined waiting times is exceeded. The longest waiting time of the road section 15 is 1.660 min; the waiting time for the section 19 is the shortest, 0.866 min. Although an approximately optimal solution can be given based on genetic algorithms, care should be taken to select optimal parameters, i.e. algebra, population number, cross probability and mutation probability.