发明内容

技术问题：本发明要解决的技术问题是如何控制高速公路网的进口匝道流量，使得下游的入城通道的最大排队等待时间不超过预定值，并使高速公路网的通行能力最大化，这样既然能满足疫情防控的需要，又能最大化公路的运输能力。

技术方案：道路流量在进口匝道处被间接控制，匝道控制是一个具有领导者-跟随者决策结构的 Stackelberg博弈。上层的运营者通过对进口匝道进行控制来最大程度地提高通行量。运营者可以预测但不能控制高速公路用户包括目的地选择和路径选择在内的出行行为，而所有用户均以用户最优的方式做出自己的出行决定。另外，下层用户的决策是在上层决策之后做出的。但是，运营者必须预期用户的行为反应来调整其决策。本发明提出一个双层优化模型来描述运营者和用户之间关系的领导者和跟随者的大小，概念框架如图1所示。其上层是具有排队延迟约束的网络通行量的最大化，而下层是带有排队的交通系统平衡，可描述预测高速公路使用者的出行行为。注意，平衡是出行分布和交通分配之间的反馈过程。本发明包括以下步骤：

(一)带排队的交通系统平衡模型

下层模型是将出行分布和交通分配模型相结合的交通系统均衡。长期以来，一直有批判说，行驶时间在传统的四阶段顺序模型中是不一致的，因为行驶时间实际上是内生确定的。一般而言，根据文献可以采用两种方法来解决这个不一致的问题，以实现交通系统的平衡。一种方法是将几个步骤组合成一个等价的数学规划，可以证明结果具有良好的收敛性和一致性(Oppenheim，1995；Sheffi，1985)。另一种是迭代地反馈顺序模型，直到行驶时间满足一致性标准(Boyce和Zhang，1997；Boyce等，1994)。尽管前者在文献中被普遍采用，但后者在每个步骤上都更加灵活(Lin and Wei，2019；Lin，2019)。因此，这里采用第二种方法，即带反馈的顺序模型。

注意，此处的交通分配不是传统的分配方式，因为此处的交通分配包含防疫封锁线上的排队延迟。排队延迟时间的确定是一个关键问题，而排队论通常是分析排队等待成本的绝佳工具。在大多数交通情形下，两次输入的时间间隔和服务时间由指数分布随机描述。基于此，本发明此处采用基于泊松分布假设的随机排队模型，即到达间隔和服务时间遵循指数分布。具体排队模型的推导是基于排队情况的稳定状态，是在系统运行足够长的时间后实现的。

根据传统的交通流理论(Gartner等，1999)，每个收费站的排队现象可以用M/M/c排队模型描述，其中M代表马尔可夫(或泊松)到达或离开的分布，或等价的指数到达或服务时间分布，c代表服务率相等的并行服务台的数量，每个进口路段上可能有一个或多个并行检查站(即服务台)。假定在一条防疫封锁线上有m个进口，车辆根据泊松过程以预测的流入量λ_i到达每个进口i(i＝1，2，...，m)，并且每个进口i 的c_i个并行检查站具有相同的参数为μ的指数分布服务时间，其中c_iμ＞λ_i。

类似于单一服务设施，在给定的进口i上，最常用的排队状况的度量是预期的排队车辆数量(l_i)和预期的排队延迟时间(d_i)。l_i与d_i之间的关系称为Little公式，公式为l_i＝λ_id_i，该关系在相当普遍的条件下有效。当ρ_i＝λ_i/μ时，表达式l_i可以确定如下：

其中，p_i是无客户在进口i的稳态概率，A^*为进口路段的集合，式(3)为稳态条件。排队等待时间d_i可以根据Little公式通过l_i除以λ_i确定。具体公式为：

进口路段的行驶时间由两部分组成。一部分是由路段交通流量决定的路段行驶时间，另一部分是由交通流和并行检查站数量决定的排队延迟时间。注意，假设队列的物理长度为零，并且没有排队后溢，这意味着路段行驶时间与队列的长度无关，在路段所需的行驶时间仅由其流量决定。即当a∈A^*时，，t_a(v_a，c_a)＝t_a(v_a)+d_a(v_a，c_a)，而当时，t_a(v_a，c_a)＝t_a(v_a)，其中t_a(v_a)是路段行驶时间，d_a(c_a)是等待时间。因此对于给定的检查站部署方案，带排队的交通分配是个常规问题。

总的来说，对于给定出行需求和路网，下层模型具体如图2所示。首先，通过目的地选择生成出行分布矩阵。其中多项式logit模型用于目的地选择，它也被认为是最简单、最实用的离散选择模型。在生成出行分布矩阵之后，通过带排队的用户均衡将出行需求分配到路网中以生成道路交通流。注意，进口路段上的预测流量是排队系统的平均到达率。然后，由Dijkstra算法计算出所有始发地与目的地(OD对)的出行时间，包括路段行驶时间和排队延迟时间。这些路径行驶时间被反馈到多项式logit模型以更新出行分布矩阵。重复此过程，直到出行分布矩阵不再变化为止。达到的这种状态即为交通系统平衡。图2说明了下层模型的反馈过程。其中使用的变量符号定义如下：

q_rs：始发地r和目的地s之间的出行需求；

O_r：区域r内的出行需求；

S_r：从始发地r出发的出行者的目的地集合；

β_s：出行者对目的地s的偏好；

t_rs：始发地r和目的地s之间的路径行驶时间；

A：网络中的所有路段的集合；

A^*：是进口路段的集合；

β_t：行驶时间t_rs的系数；

v_a：路段a上的交通流；

c_a：路段a上的检查站数；

t_a：在路段a上的行驶时间，它是交通流v_a和检查站数量c_a的函数；

连接始发地r和目的地s的路径k上的交通流；

路段路径关联关系表示为：

(二)双层规划模型

对于一个给定的上层进口匝道控制方案，下层存在一个带排队的交通系统平衡。运营者的目标是通过控制进口匝道流量，使各出口匝道处排队等待时间在可容忍范围内，并使网络通行量最大化。可接受的控制方案应考虑排队延迟约束。注意，给定出口匝道i的服务水平是路段流λ_i的函数，而路段流量可以被进口匝道控制器控制。平均排队时间为确定可接受的进口匝道控制率提供了一个约束条件，即

d_i(λ_i，c_i)≤T，

常数T是决策者指定的期望水平，如T＝3min。注意，d_i是路段流λ_i和检查站数量c_i函数。根据等式(4)可将平均等待时间d_i的约束更具体地确定为

Yang等(1994)提出了一种双层优化方法，其中进口匝道交通控制问题被描述为Stackelberg博弈。受他们研究的启发，本发明也采用了双层规划模型。上层的高速公路运营者旨在以检查站处的排队延迟为约束，使通过所有匝道的总流量最大。运营者可以预测，但无法控制高速公路用户包括目的地选择和路线选择的出行行为，而所有用户都以用户最优的方式做出决定。由此，上层模型的公式确定为：

其中的符号定义如下：

u_r：通过进口匝道的流入量r∈R；

进口匝道的交通需求r∈R；

R：进口匝道的集合，即始发地；

I：出口匝道的集合；

u：所有流入的向量(上层决策变量)。

在上层模型中，目标函数(6)是通过所有进口匝道的总输入流量最大化；约束(7)意味着每个出口匝道处的等待时间都不应超过预定的可接受水平T，注意路段流量λ_i(u)是所有匝道流入量的函数；约束(8)是稳态条件；而约束条件(9)表示每个匝道的流入量应为非负且小于等于对应的交通需求。而要获得路段流量λ_i(u)，需要先解决下层问题。

(三)带排队的交通系统平衡模型算法

要求解所提出的双层规划模型，最好先求解下层的模型，因为下层模型是内嵌在上层模型中的。在固定的出行需求和已建成的路网的情况下，对上层给定的决策，下层将会有一个稳定的交通流模式。需要注意的是，下层是出行分布和带排队的交通分配之间的反馈过程。其中，连续平均算法(MSA)可用于实现系统平衡。初始出行分布矩阵可以通过具有初始化的始发地-目的地(OD)对行驶时间的多项式logit模型产生，然后通过Frank-Wolfe算法将需求分配给路网，可以生成路段交通流和路段行驶时间。此外，还可以通过预测的进口路段交通流来确定排队延迟时间。每个进口路段的广义行驶时间包括路段行驶时间和排队延迟时间。根据Wardrop路线选择的第一原理，也就是所谓的用户均衡，流量在拥挤的网络中自行安排，使OD对之间所有可使用的路径具有相同且最低的成本。因此，接下来可以用Dijkstra算法更新OD对的最短行驶时间，然后将这些时间反馈到多项式logit模型以生成新的出行分布矩阵。但是，该矩阵不能直接分配给路网，因为直接反馈的收敛通常是不能实现的，必须对连续出行分布矩阵求平均值。尽管也有一些常数权重的成功应用，但是这样通常无法保证收敛。因而本发明此处采用权重递减的MSA来更新出行分布矩阵，权重是迭代次数的倒数。将更新后的矩阵进一步分配给路网，迭代过程一直持续到连续矩阵为近似相等为止。收敛通常通过连续出行需求矩阵之间相对误差的平方根来衡量。如果满足预定的误差值，则终止迭代，稳定状态即称为交通系统平衡。由此产生的交通流结果随后进入上层模型中。图3是带排队的系统均衡算法的流程图。详细的MSA算法步骤如下所示：

步骤1输入来自上层的进口匝道控制方案u和建成的路网。

步骤2基于用初始OD对行驶时间初始化出行分布矩阵令n＝1作为迭代次数。

步骤3带排队的交通分配。出行分布矩阵通过Frank-Wolfe算法分配给路网。生成路段行驶流v_a和路段行驶时间t_a。注意，排队延迟时间可以通过进口路段上的流量来确定。

步骤4通过Dijkstra算法更新OD对r-s间(即)的最短行驶时间。

步骤5出行分布。多项式logit模型用于更新出行分布矩阵

步骤6使用递减权重的方式平均出行分布矩阵和

步骤7收敛识别。使用相对误差的平方根检验出行分布矩阵的收敛性

其中，ε是预先确定的误差值。如果满足收敛条件，则终止迭代并转到步骤9，否则转到步骤8。

步骤8使且n：＝n+1，然后转到步骤3。

步骤9输出出行分布矩阵和路段交通流v_a，返回上层模型。

(四)双层规划模型的算法

双层规划问题是众所周知的NP困难问题，很难用经典的优化算法来解决。即使上层和下层都是线性规划，也具有一定的挑战性，更何况上层是非线性整数规划模型。例如，传统的基于梯度的最优边界线收费问题的求解方法，由于存在多个最优解，对于规模较大的问题往往无法收敛。这种失败导致了一种确定最佳收费水平和收费位置的启发式算法的发展。因而，启发式方法，尤其是遗传算法，在文献中有许多成功的应用(Liu等，2013；Shepherd和Sumalee，2004；Sumalee，2004，2007)。尽管启发式算法很耗时，而且没有办法确保全局最优，但它已被证明可以成功解决边界线通行费优化的问题。因此，这里采用了一种基于精英策略的遗传算法。图4为其算法流程图。详细的基于精英策略的遗传算法步骤如下：

步骤1参数初始化。设置遗传算法中使用的参数，包括种群规模M，最大代数Gen，交叉概率 p_c，突变概率p_m，代数符号gen＝1，精英保留概率p_e。请注意，种群规模取决于问题的性质，但通常包含数百个可行的解。

步骤2随机生成一个可行的初始种群。一个染色体是由m个基因组成的一个解，每个基因表示单个进口匝道运行通过的流量。随机生成一条染色体，如果不可行，继续生成直到可行为止。最后产生总计M个可行的染色体，分散在所有可能解的范围。

步骤3选择操作。上层模型的目标函数用作适应性函数以评估种群中所有染色体的性能。注意，这里是为了最大程度地提高网络通行量，由保留概率p_e标注精英群体，适应度最差的p_e群体则被舍弃。

步骤4交叉操作。剩下的(1-p_e)M染色体用于交叉操作，这些父辈染色体随机配对，进行交叉的概率为p_c。如果进行交叉，则确定一个随机基因的位置进行交叉。如果根据上层模型中的约束，新生染色体不可行，则尝试其他基因位置，直到可行为止。这些新生成的染色体通常具有其父代的许多特征。

步骤5变异操作。确定一个随机基因位置在定义域内进行突变，进行突变的概率为p_m。如果新染色体不可行，则尝试其他基因位置直到可行。

步骤6生成下一代群体。在遗传操作之后，仍然存在可行的(1-p_e)M个染色体。添加标注的p_eM 个精英群体以确保种群总数达到M，使当前一代中最好的染色体能够继续保留到下一代不变，从而保证了解的质量不会因遗传迭代而降低。更新代数的表示为gen：＝gen+1。

步骤7终止判断。如果代数达到最大值，即gen≥Gen，则终止迭代过程并输出防疫封锁线的最佳部署方案。否则，转到步骤3。

(五)程序设计

timestart＜-Sys.time()#用于计时的

#步骤1：初始化，按格式输入数据和必要的包；

#1.1加载计算最短路径的包，准备调用dijkstra最短路径算法，注意igraph包首次使用需要安装，然后才能调用。

#install.packages(″igraph″)#安装igraph包

library(igraph)

options(digits＝3)#保留三位小数

#1.2创建图的距离矩阵，包含所有的候选路段。第一列为路段标号(Road)，第二列为路段起点标号(Road origin)，第三列为路段终点标号(Road destination)，第四列为该路段自由流时间(free flow time)，第五列为道路通行能力(capacity)，第六列为道路长度(length)。此处以交通配流中常用的Nguyen-Dupuis网络为例，详细的参数设置可参考程序文档。

#也可以在Excel中复制Table 1，然后执行

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

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，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，14，10，11，6.0，800，4.00，15，11，2，9.0，800，6.00，16，11，3，8. 0，800，5.00，17，12，6，7.0，800，4.00，18，12，7，14.0，800，6.00，19，13，3，11.0，800，7.00)，ncol＝6，byrow＝T)#为避免0出现在分母，用一个充分小的数代表0，这样该路段的阻抗充分大，确保没有流量.

e＝matrix

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

#e#用于检查程序的断点

#1.3输入初始交通需求矩阵d0，第一列为起讫点对的标号(OD pair)，第二列为起点标号(origin)，第三列为终点标号(destination)，第四列为交通需求(demand)。

sg＝c(1000，1000)#各地期待的交通需求，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，byrow＝T)#初始分配方案 colnames(d0)＝c(″OD pair″，″Origin″，″Destination″，″Demand″)

#d0#用于检查程序的断点

#后面用到的参数值

M＝200#如果到达率超出了排队系统的服务能力，则假设其排队时间为大M

el＝c(11，15，16，19)#入城通道(entry link)的集合

miud＝c(1500，1500)#各个出行地的交通需求

#cj＝rep(9，4)#每个入城通道允许的服务台数

cj＝c(9，3，5，5)#目前的安全检查站最大配置能力

mu＝2#每个服务台的平均服务率(每分钟测2台车)

AT＝2#Acceptable level of service.2minutes

#以下是路段更新函数中会用到的参数，原本定义在fw函数内部，但是由于在其他函数中也会用到，所以将其移到此处，定义为全局变量

t0＝e[，4]#自由流时间

c＝e[，5]#道路通行能力

a＝0.15

b＝4

#自定义的Frank-Wolfe算法函数，注意输入的需求矩阵d形式如d0，交通网络e的形式如上面的e，相对误差0.01。注意对于入城通道，其时间

#为路段出行时间和排队时间之和

fw＝function(e，d，s)

{

#1.4根据路径自由流时间计算各个OD对的最短路径和路径流量

g＝add.edges(graph.empty(13)，t(e[，2：3])，weight＝e[，4])#创建图，13为节点的个数，以时间为权重而非路径的长度

b12＝get.shortest.paths(g，from＝″1″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点2的最短路径

b13＝get.shortest.paths(g，from＝″1″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点3的最短路径

b42＝get.shortest.paths(g，from＝″4″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点2的最短路径

b43＝get.shortest.paths(g，from＝″4″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点3的最短路径

#创建一个矩阵，用于保存各个OD对的最短路径和流量

V＝cbind(e[，1])

st0＝numeric(4)#存放初始的各OD对最短行驶时间

colnames(V)＝″Road″

V

#OD对12的最短路径和流量

sp12＝as.vector(b12)#转化为路段标号(Road)

st0[1]＝sum(e[sp12，4])#各路段时间求和

x12＝cbind(e[sp12，1]，rep(d[1，4]，length(sp12)))#路段标号和流量，算法中的迭代起点

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

x12

V＝merge(V，x12，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵

V[is.ha(V)]＝0

V

#OD对13的最短路径和流量

sp13＝as.vector(b13)#转化为路段标号(Road)

st0[2]＝sum(e[sp13，4])#各路段时间求和

x13＝cbind(e[sp13，1]，rep(d[2，4]，length(sp13)))#路段标号和流量，算法中的迭代起点

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

x13

V＝merge(V，x13，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵

V[is.ha(V)]＝0

V

#OD对42的最短路径和流量

sp42＝as.vector(b42)#转化为路段标号(Road)

st0[3]＝sum(e[sp42，4])#各路段时间求和

x42＝cbind(e[sp42，1]，rep(d[3，4]，length(sp42)))#路段标号和流量，算法中的迭代起点

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

x42

V＝merge(V，x42，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵

V[is.na(V)]＝0

V

#OD对43的最短路径和流量

sp43＝as.vector(b43)#转化为路段标号(Road)

st0[4]＝sum(e[sp43，4])#各路段时间求和

x43＝cbind(e[sp43，1]，rep(d[4，4]，length(sp43)))#路段标号和流量，算法中的迭代起点

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

x43

V＝merge(V，x43，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵

V[is.aa(V)]＝0

V

#当所有最短路径上的流量求和，得到初始流量

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

VS

#步骤2：更新各路段的阻抗

#以下的参数原本定义在该子函数内部，但由于其他子函数也会使用，就会出现找不到对象的情况，所以将其放在子函数的外面

#定义为全局变量

#t0＝e[，4]#自由流时间

#c＝e[，5]#道路通行能力

#a＝0.15

#b＝4

tp＝function(v){#s为排队系统的服务台数向量

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

#对于入城通道要加上排队时间

lamda＝v[el]/60#入城通道的平均车流量(从每小时转化为每分钟)

rho＝lamda/mu

td＝numeric(length(lamda))#定义一个变量存放逗留时间

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

td[lamda＝＝0]＝0

#上述为Little公式，通过排队长计算排队时间，此时应该注意，如果流入量为0，

#而分母不能为0，程序会报错，因此对于lamda为0应该特别处理一下

td[rho/s＞＝1]＝M#定义排队时间，将不满足平稳条件的时间赋值为大M

tf[el]＝tf[el]+td

tf

}

k＝1#控制循环次数，达到该次数后强行退出，防止死循环(来回震荡)

repeat{

#步骤3：寻找下一个迭代方向

g2＝add.edges(graph.empty(13)，t(e[，2：3])，weight＝tp(VS))#构造图，13为节点的个数，更新路段阻抗

b12＝get.shortest，paths(g2，from＝″1″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点2的最短路径

b13＝get.shortest.paths(g2，from＝″1″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点3的最短路径

b42＝get，shortest.paths(g2，from＝″4″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点2的最短路径

b43＝get.shortest.paths(g2，from＝″4″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点3的最短路径

#创建一个临时矩阵，用于保存各个OD对的最短路径和流量

V＝cbind(e[，1])

colnames(V)＝″Road″

V

#OD对12的最短路径和流量

sp12＝as.vector(b12)#转化为路段标号(Road)

x12＝cbind(e[sp12，1]，rep(d[1，4]，length(sp12)))#路段标号和流量，算法中的迭代起点

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

x12

V＝merge(V，x12，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵

V[is.na(V)]＝0

V

#OD对13的最短路径和流量

sp13＝as.vector(b13)#转化为路段标号(Road)

x13＝cbind(e[sp13，1]，rep(d[2，4]，length(sp13)))#路段标号和流量，算法中的迭代起点

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

x13

V＝merge(V，x13，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵

V[is.na(V)]＝0

V

#OD对42的最短路径和流量

sp42＝as.vector(b42)#转化为路段标号(Road)

x42＝cbind(e[sp42，1]，rep(d[3，4]，length(sp42)))#路段标号和流量，算法中的迭代起点

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

x42

V＝merge(V，x42，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵

V[is.na(V)]＝0

V

#OD对43的最短路径和流量

sp43＝as.vector(b43)#转化为路段标号(Road)

x43＝cbind(e[sp43，1]，rep(d[4，4]，length(sp43)))#路段标号和流量，算法中的迭代起点

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

x43

V＝merge(V，x43，by＝″Road″，all＝TRUE)#定义V为专门保存迭代起点的矩阵

V[is.na(V)]＝0

V

#当所有最短路径上的流量求和，得到迭代方向

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

VS2

#步骤4：计算迭代步长

step＝function(lamda){

x2＝VS2

x1＝VS

q＝x1+lamda*(x2-x1)

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

}

#lamda＝uniroot(step，c(0，1))$root#注意lamda的取值范围，步长不能太长，uniroot要求两端的函数值符号相反，有的函数不一定满足。采用optimize函数可以确保找到一元函数的最优值。

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

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

#print(lamda)

#步骤5：确定新的迭代起点

VS3＝VS+lamda*(VS2-VS)

VS3

#步骤6：收敛性检验

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

VS＝VS3#如果不满足收敛条件则用新点VS3替代原点VS，如此循环直到收敛

k＝k+1

if(k＝＝100)break#循环达到100次就强行退出

}

#步骤7：输出平衡状态的特征矩阵result和OD行驶时间矩阵u。

#步骤7.1：输出平衡状态各路径的流量、通行时间和速度。

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″)

#步骤7.2：输出各OD行驶时间矩阵u

g＝add.edges(graph.empty(13)，t(e[，2：3])，weight＝result[，3])#创建图，13为节点的个数，result为步骤7生成的矩阵

b 12＝get.shortest.paths(g，from＝″1″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点2的最短路径

b13＝get.shortest.paths(g，from＝″1″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点1到终点3的最短路径

b42＝get.shortest.paths(g，from＝″4″，to＝″2″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点2的最短路径

b43＝get.shortest.paths(g，from＝″4″，to＝″3″，mode＝″out″，output＝″epath″)$epath[[1]]#从起点4到终点3的最短路径

#创建一个行驶时间矩阵，用于保存各个OD对的行程时间，初始假设各OD行程时间为0

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

#OD对12的行程时间

sp12＝as.vector(b12)#转化为路段标号(Road)

u[1，4]＝sum(result[sp12，3])#各路段时间求和

#OD对13的行程时间

sp13＝as.vector(b13)#转化为路段标号(Road)

u[2，4]＝sum(result[sp13，3])#各路段时间求和

#OD对42的行程时间

sp42＝as.vector(b42)#转化为路段标号(Road)

u[3，4]＝sum(result[sp42，3])#各路段时间求和

#OD对43的行程时间

sp43＝as.vector(b43)#转化为路段标号(Road)

u[4，4]＝sum(result[sp43，3])#各路段时间求和

u＝cbind(u，st0)#OD对间无障碍行驶时间st0

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

#以列表的形式输出result矩阵和OD行驶时间矩阵

list(result，u)

}

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

#cj＝rep(1，4)

#fw(e，d0，cj)#用于检查程序的断点

#步骤8：定义目的地选择的多项式Logit函数mlogit，输入为各OD行驶时间时间矩阵u和各地交通需求 sg，输出为新的交通分布矩阵。

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)

}

#步骤9：定义一个给定交通需求d0和sg及交通网络e下综合的交通分布与交通分配交替迭代平衡函数。对于初始交通分布可以求得用户平衡状态时各OD的行驶时间矩阵，用户根据该矩阵重新选择目的地，对于新的交通分布又可以生成新的行驶时间矩阵，该过程一直循环进行，一直到交通分布矩阵不再变化为止。 cda＝function(e，d0，sg，s){

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#满足一定的相对误差要求就停止

d0[，4]＝d0[，4]+(1/k)*(d1[，4]-d0[，4])#这里采用迭代加权法(Method ofSuccessive Average，MSA)，此处采用循环次数的倒数作为权重，随着循环次数的增加而减少。

#print(d0)#用于检查程序的断点

#cat(k)#用于检查程序

if(k＝＝100)break#如果循环次数达到100次但还没有满足精度要求也跳出循环

}

#print(d1)

d2＝fw(e，d1，s)

#print(d2)

#输出入城通道的出行时间和排队时间

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

tf＝t0*(1+a*(v/c)^b)#路段出行时间，注意这里的参数t0，a，c，b在配流子函数fw中定义过，但是由于定义为局部变量，在这个

#子函数中不能直接使用，系统会报错找不到这些对象，解决方法是在此处重新定义，或者请其定义为全局变量，这样各个子函数

#都可以使用了，这里将其定义为全局变量。

lamda＝v[e1]/60#入城通道的平均车流量(从每小时转化为每分钟)

rho＝lamda/mu

td＝line(s，rho)/lamda

td[rho/s＞＝1]＝M#定义排队时间，将不满足平稳条件的时间赋值为大M

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

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

#将列表输出

list(d1，d2，d3)

}

#系统平衡状态

#cda(e，d0，sg，cj)

#cj＝rep(5，4)#测试数据

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

bef＝cda(e，d0，sg，cj)#初始的平衡状态before-after analysis，用于检查程序到此处是否正常运行

bef

##########################求解双层规划的遗传算法#############################

#步骤1：Initialization参数初始化

M2＝500#每代染色体的数目，本题生成可行解较慢，M2不宜太大

Gen＝30#最多代数

pc＝0.5#交叉概率

pm＝0.1#变异概率

pe＝0.1#精英占比

gen＝1#标记代数

#定义计算每个通道逗留时间的函数，后面好调用，输入变量为服务台向量s，服务强度向量rho，输出为队长向量lq。

#需要注意的是要确保输入的s和rho满足rho.s＝rho/s＜1的平衡条件，否则函数结果不正确。

line＝function(s，rho)#计算排队长的子函数

{

rho.s＝rho/s#单个服务台的服务强度

p0＝rep(0，length(s))#初值，顾客数为0的概率

lq＝rep(0，length(s))#初值，平均排队长

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

{

sum1＝0

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

{

sum1＝sum1+rho[k1]^k2/factorial(k2)#factorial函数计算阶乘

}

p0[k1]＝(sum1+rho[k1]^(s[k1])/(factorial(s[k1])*(1-rho.s[k1])))^(-1)#系统中顾客数为0的概率p0

lq[k1]＝rho[k1]^s[k1]*rho.s[k1]*p0[k1]/(factorial(s[k1])*(1-rho.s[k1])^2)#将c(s，p)的公式带入排队长Lq 的公式

}

return(lq+rho)#计算平均队长，而不是排队长

}

#步骤2：Generate a feasible initial population randomly；

pop＝matrix(nrow＝M2，ncol＝length(miud))#用于存放每一代的种群

gep＝matrix(nrow＝Gen，ncol＝length(miud)+2)#用于存放每一代的最优染色体及其表现，包含了代数、染色体、适应度，所以列数为ength(lamda)+2

i＝1#计数器，标记每一代中的染色体

repeat

{

#随机生成一条染色体

chr＝numeric(length(miud))#初始值，这里的chr就是表示决策变量

for(j in 1：length(miud))chr[j]＝sample(seq(1：miud[j])，size＝1)#随机生成一个初始解

#检查其是否符合约束条件，如果不符合的话重新生成新的染色体

bef＝cda(e，d0，chr，cj)#计算平衡流量

lamda＝bef[2]][[1]][c(11，15，16，19)，2]/60#入城通道的平均车流量(从每小时转化为每分钟)

rho＝lamda/mu

td＝numeric(length(lamda))#定义一个变量存放逗留时间

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

td[lamda＝＝0]＝0

#上述为Little公式，通过排队长计算排队时间，此时应该注意，如果流入量为0，

#而分母不能为0，程序会报错，因此对于lamda为0应该特别处理一下

td[rho/cj＞＝1]＝M#定义排队时间，将不满足平稳条件的时间赋值为大M

if(max(td)＞AT)next#确保满足rho.s＝rho/s＜1的平衡条件；任何一个入城通道排队时间不能大于AT

#print(chr)#输出可行解，用于查看中间计算过程的

pop[i，]＝chr

cat(i，chr)#输出可行解，用于查看中间计算过程的

cat(″\n″)#输出换行符

i＝i+1

if(i＞M2)break

}

repeat{

#步骤3：selection operation执行选择操作

per＝matrix(nrow＝M2，ncol＝2)#用于存放每个染色体的适应度函数值，第一列为染色体的编号

for(i in 1：M2)

{

#对每个染色体计算其适应度函数

per[i，]＝c(i，sum(pop[i，]))#找出所有入城通道吞吐量最大的那个作为该分配方案的表现

print(c(gen，per[i，]))#用于监控计算过程

}

#舍弃表现最差的Pe，标记表现最好的Pe，以便后面使用

per2＝per[order(per[，2]，decreasing＝TRUE)，]#该问题是总交通量越大越好，按照从大到小排

per3＝per2[1：(M2*(1-pe))，]#留下1-pe的进行后面的交叉和变异，舍弃最后pe的染色体

per4＝per2[1：(M2*pe)，]#表现最好的精英群体

#保存每一代最好的表现，以便后期绘图

gep[gen，]＝c(gen，pop[per2[1，1]，]，per2[1，2])#注意第一列为代数，最后一列为适应度，中间为本代的最佳染色体

#步骤4：Crossover operation

#进行交叉和变异的染色体

pop2＝pop[sample(per3[，1])，]#sample函数进行随机排序，然后两两成对

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

{

if(runif(1)＞pc)next#交叉概率，如果随机数大于pc就不交叉，小于pc就交叉

k＝＝0#计数器，防止死循环

repeat

{

c1＝sample(length(miud)，1)#选择一个进行单点交叉的位置

c2＝pop2[i，c1]#c2是两个基因交换时的中间变量

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

pop2[i+1，c1]＝c2

#print(pop2[i，])

#循环控制，防止死循环

k＝k+1

if(k＞20)break

#要求交换以后的两个新染色体是可行解，否则继续随机选择位置进行交换

bef＝cda(e，d0，pop2[i，]，cj)#计算平衡流量

lamda＝bef[[2]][[1]][c(11，15，16，19)，2]/60#入城通道的平均车流量(从每小时转化为每分钟)

rho＝lamda/mu

td＝numeric(length(lamda))#定义一个变量存放逗留时间

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

td[lamda＝＝0]＝0

#上述为Little公式，通过排队长计算排队时间，此时应该注意，如果流入量为0，

#而分母不能为0，程序会报错，因此对于lamda为0应该特别处理一下

td[rho/cj＞＝1]＝M#定义排队时间，将不满足平稳条件的时间赋值为大M

if(max(td)＞AT)next#确保满足rho.s＝rho/s＜1的平衡条件；任何一个入城通道排队时间不能大于AT

bef＝cda(e，d0，pop2[i+1，]，cj)#计算平衡流量

lamda＝bef[[2]][[1]][c(11，15，16，19)，2]/60#入城通道的平均车流量(从每小时转化为每分钟)

rho＝lamda/mu

td＝numeric(length(lamda))#定义一个变量存放逗留时间

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

td[lamda＝＝0]＝0

#上述为Little公式，通过排队长计算排队时间，此时应该注意，如果流入量为0，

#而分母不能为0，程序会报错，因此对于lamda为0应该特别处理一下

td[rho/cj＞＝1]＝M#定义排队时间，将不满足平稳条件的时间赋值为大M

if(max(td)＜＝AT)

{print(pop2[i，])#查看计算过程

print(pop2[i+1，])#查看计算过程

break}

#确保满足rho.s＝rho/s＜1的平衡条件；任何一个入城通道排队时间不能大于AT

#注意在进行比较时，&&如果发现左边对象的值为FALSE，那么就不会计算右边的对象了，&& 的这种计算方法叫做短路计算。

}

}

#步骤5：Mutation operation执行变异操作

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

{

if(runif(1)＞pm)next#变异概率，如果随机数大于pm就不变异，否则变异

repeat

{

c3＝sample(length(miud)，1)#选择一个进行单点变异的位置

#将该位置的基因值变为取值范围内的一个随机值

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

print(pop2[i，])

#直到可行解才退出

bef＝cda(e，d0，pop2[i，]，cj)#计算平衡流量

lamda＝bef[[2]][[1]][c(11，15，16，19)，2]/60#入城通道的平均车流量(从每小时转化为每分钟)

rho＝lamda/mu

td＝numeric(length(lamda))#定义一个变量存放逗留时间

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

td[lamda＝＝0]＝0

#上述为Little公式，通过排队长计算排队时间，此时应该注意，如果流入量为0，

#而分母不能为0，程序会报错，因此对于lamda为0应该特别处理一下

td[rho/cj＞＝1]＝M#定义排队时间，将不满足平稳条件的时间赋值为大M

if(max(td)＜＝AT)break

}

}

#步骤6：Generate the next generation population.

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

pop＝pop3#返回开始执行迭代

gen＝gen+1

#步骤7：Termination judgment

if(gen＞Gen)break

}

#步骤8：输出需要的计算结果

gep#每代的表现，第一列为代数，最后一列为总的交通吞吐量，中间为决策变量

write.csv(gep，file＝″基因演化过程.csv″)#以csv格式保持到当前工作目录

#使用gep的最后一行信息，输出最优的计算结果

bef＝cda(e，d0，gep[Gen，2：3]，cj)#计算平衡流量

bef#显示计算结果

gep[Gen，2：3]#网络最大通行能力

lamda＝bef[[2]][[1]][c(11，15，16，19)，2]/60#入城通道的平均车流量(从每小时转化为每分钟)

rho＝lamda/mu

rho/cj

td＝numeric(length(lamda))#定义一个变量存放逗留时间

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

td[lamda＝＝0]＝0

#上述为Little公式，通过排队长计算排队时间，此时应该注意，如果流入量为0，

#而分母不能为0，程序会报错，因此对于lamda为0应该特别处理一下

td[rho/cj＞＝1]＝M#定义排队时间，将不满足平稳条件的时间赋值为大M

td

###计算程序的运行时间

timeend＜-Sys.time()

runningtime＜-timeend-timestart

print(runningtime)#输出运行时间

设置城市防疫封锁线是阻断病毒传播的有效手段，出行的人们只有经过检测后才能通过封锁线。本发明提出了一种通过控制高速公路网进口匝道流量，以使防疫封锁线上的各入城通道的最大排队等待时间不超过预定值，并最大化公路网通行能力的软件。本发明能够有效避免防疫封锁线上的过长排队等待时间，可以帮助网络管理者确定有效的进口匝道流量控制方案。