[发明专利]考虑新能源汽车渗透的环境影响评估软件

权利要求书

1.本发明采用R语言开发了一种考虑新能源汽车渗透的环境影响评估软件，该软件的算法框架具体包括以下步骤：

步骤一：定义传统内燃机车(ICEV)、插电式电动汽车(PEV)和混合动力电动汽车(HEV)的排放率，将其表示为行驶速度的函数；

步骤二：使用吸收马尔可夫链模型表示ICEV、PEV和HEV的市场占有率演化过程，其中PEV是吸收状态，表示消费者选择PEV后不再会选择其他类型的汽车；

步骤三：通过交通产生、交通分布、交通方式划分和交通流分配的反馈迭代达到交通系统平衡，计算出平衡状态时的路段交通流量和通行时间；

步骤四：使用各路段的长度和行驶时间计算各路段的平均速度；

步骤五：评估交通系统的环境影响。

2.对于权利1中的步骤三，采用如下的具体计算过程：

步骤1：输入交通网络信息，包括网络结构、路段长度、路段通行能力等；

步骤2：给定交通需求O_i，通过平均分布得到初始化交通分布矩阵设置n＝0，表示迭代次数；

步骤3：通过Frank-Wolfe算法基于用户均衡将交通分布矩阵分配给交通网络，以计算每个路段a上的交通流量和出行时间，之后，起点i和目的地j之间的最短出行时间，即可以通过Dijkstra算法算出；

步骤4：基于采用目的地选择模型来更新交通分布矩阵

步骤5：利用权重递减的连续平均法(MSA)对交通分布矩阵和求平均

步骤6：使用相对根平方误差(RRSE)检查交通分布矩阵的收敛性

如果满足收敛条件，则转至步骤8，否则令n＝n+1，转至步骤7；

步骤7：通过Frank-Wolfe算法基于用户均衡将交通分布矩阵分配给交通网络，以计算每个路段a上的交通流量和出行时间，之后，起点i和目的地j之间的最短出行时间，即可以通过Dijkstra算法计算，计算结果被反馈到步骤4；

步骤8：输出交通分布矩阵以及路段a上的交通流量v_a和起点i与目的地j之间的出行时间

3.权利1和权利2中的算法以R语言编写，具体归纳如下：

#考虑新能源汽车的渗透对环境影响的评估软件

#步骤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中复制，然后执行

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

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

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

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

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

tge＝3000#总的现状交通需求

d0＝matrix(c(1，1，2，0.2*tge，2，1，3，0.4*tge，3，4，2，0.3*tge，4，4，3，0.1*tge)，ncol＝4，byrow＝T)#初始分配方案

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

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

#自定义的Frank-Wolfe算法函数，注意输入的需求矩阵d形式如d0，交通网络e的形式如上面的e，相对误差0.001；

fw＝function(e，d)

{

#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.na(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.na(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.na(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){

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

}

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.ha(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

lamda

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

VS3＝VS+lamda*(VS2-VS)

VS3

#步骤6：收敛性检验

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

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

}

#步骤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″，″Levelof Service″)

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

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

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对的行程时间，初始假设各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

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

list(result，u)

}

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

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

mlogit＝function(u，sg)

{

d＝numeric(4)

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

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

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

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

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

}

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

cda＝function(e，d0，sg){

k＝3

repeat{

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

k＝k+1

if(sqrt(sum((d1[，4]-d0[，4])^2))/sum(d0[，4])＜0.01)break#满足一定的精度要求就停止

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

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

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

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

}

#print(d1)

d2＝fw(e，d1)

#print(d2)

list(d1，d2)

}

#步骤10：输出交通系统的表现；

ge1＝c(sum(d0[1：2，4])，sum(d0[3：4，4]))#用于检查程序的断点

d3＝cda(e，d0，ge1)#对交通需求ge1和交通分布d01调用前面定义的cda函数

write.csv(d3[[1]]，file＝″交通分布矩阵.csv″)#以csv格式保持到当前工作目录

write.csv(d3[[2]][[1]]，file＝″路段平衡结果.csv″)

write.csv(d3[[2]][[2]]，file＝″OD出行时间.csv″)

getwd()#查看输出文件的保存地址

#步骤11：利用Markov chain模型计算各类型汽车的市场占有率；

u＝c(0.7，0.2，0.1)

p0＝matrix(c(0.8559，0.0810，0.0631，0.0554，0.8647，0.0799，0，0，1)，byrow＝T，nc＝3)#假设的转化率矩阵

p＝p0

Q＝p0[1：2，1：2]

U＝diag(2)

N＝solve(U-Q)

est＝N％*％c(1，1)

n＝20

A＝matrix(0，nc＝3，nr＝n)#初始化矩阵

for(i in(1：n))

{

A[i，]＝u％*％p

p＝p％*％p0

}

A

MS＝rbind(u，A)#把现状市场占有率合并进去，目的为了作图；

MS

#步骤12：画出市场占有率的变化图

par(family＝′serif′)#绘图时使用新罗马字，Times New Roman是Windows下默认的衬线字体；

plot(MS[，1]，type＝″b″，ylim＝c(0，0.8)，xlab＝″The planning years″，ylab＝″Market penetration level″，cex.lab＝1.3)

axis(1，seq(0，21，1))#设定x坐标轴

axis(2，seq(0，0.8，0.1))#设定y坐标轴

lines(MS[，2]，type＝″b″，pch＝2)

lines(MS[，3]，type＝″b″，pch＝3)

legend(19，0.6，c(″ICEV″，″HEV″，″PEV″)，pch＝c(1：3))

#步骤13：计算环境影响评价EIA

EIA＝matrix(0，nc＝4，nr＝n+1)

for(i in(1：(n+1)))

{

#i＝1

#ICEV的排放量

ICEV_CO2＝3158*(e[，6]/d3[[2]][[1]][，3]*60)^(-0.56)*e[，6]*d3[[2]][[1]][，2]*MS[i，1]#注意这里的时间为min，所以乘以60以转换为hour，此时速度为mile/hour

ICEV_yOC＝1.3647*(e[，6]/d3[[2]][[1]][，3]*60)^(-0.679)*e[，6]*d3[[2]][[1]][，2]*MS[i，1]

ICEV_NOX＝2.5376*(e[，6]/d3[[2]][[1]][，3]*60)^(-0.42)*e[，6]*d3[[2]][[1]][，2]*MS[i，1]

#PEV的排放量

PEV_EC＝1.79e-8*(e[，6]/d3[[2]][[1]][，3]*60)^4-4.073e-6*(e[，6]/d3[[2]][[1]][，3]*60)^3+3.654e-4*(e[，6]/d3[[2]][[1]][，3]*60)^2-0.0109*(e[，6]/d3[[2]][[1]][，3]*60)+0.2372#耗电量

PEV_CO2＝893*PEV_EC*e[，6]*d3[[2]][[1]][，2]*MS[i，3]

PEV_NOX＝1.66*PEV_EC*e[，6]*d3[[2]][[1]][，2]*MS[i，3]

PEV_SO2＝3.79*PEV_EC*e[，6]*d3[[2]][[1]][，2]*MS[i，3]

#HEV的排放量

a＝0.6#这算系数

HEV_CO2＝a*(3158*(e[，6]/d3[[2]][[1]][，3]*60)^(-0.56)+893*PEV_EC)*e[，6]*d3[[2]][[1]][，2]*MS[i，2]

HEV_NOX＝a*(2.5376*(e[，6]/d3[[2]][[1]][，3]*60)^(-0.42)+1.66*PEV_EC)*e[，6]*d3[[2]][[1]][，2]*MS[i，2]

HEV_VOC＝a*1.3647*(e[，6]/d3[[2]][[1]][，3]*60)^(-0.679)*e[，6]*d3[[2]][[1]][，2]*MS[i，2]

HEV_SO2＝a*3.79*PEV_EC*e[，6]*d3[[2]][[1]][，2]*MS[i，2]

#总的环境评价

EIA[i，1]＝sum(ICEV_CO2+PEV_CO2+HEV_CO2)#CO2

EIA[i，2]＝sum(ICEV_NOX+PEV_NOX+HEV_NOX)#NOX

EIA[i，3]＝sum(ICEV_VOC+HEV_VOC)#VOC

EIA[i，4]＝sum(PEV_SO2+HEV_SO2)#SO2

}

#步骤14：保存在excel表格中，因为量纲不同而不容易放在一个图中

colnames(EIA)＝c(″CO2″，″NOx″，″VOC″，″SO2″)

EIAkg＝EIA/1000#转换为kg

EIAkg

write.csv(EIAkg，file＝″污染排放.csv″)

getwd()#查看输出文件的保存地址。