CN110489871A - 考虑新能源汽车渗透的环境影响评估软件 - Google Patents
考虑新能源汽车渗透的环境影响评估软件 Download PDFInfo
- Publication number
- CN110489871A CN110489871A CN201910771372.7A CN201910771372A CN110489871A CN 110489871 A CN110489871 A CN 110489871A CN 201910771372 A CN201910771372 A CN 201910771372A CN 110489871 A CN110489871 A CN 110489871A
- Authority
- CN
- China
- Prior art keywords
- road
- section
- matrix
- starting point
- flow
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 230000007613 environmental effect Effects 0.000 title claims abstract description 18
- 230000008595 infiltration Effects 0.000 title claims abstract description 12
- 238000001764 infiltration Methods 0.000 title claims abstract description 12
- 230000006870 function Effects 0.000 claims abstract description 33
- 238000000034 method Methods 0.000 claims abstract description 31
- 230000008569 process Effects 0.000 claims abstract description 15
- 238000002485 combustion reaction Methods 0.000 claims abstract description 6
- 230000008140 language development Effects 0.000 claims abstract description 4
- 239000011159 matrix material Substances 0.000 claims description 97
- 238000004422 calculation algorithm Methods 0.000 claims description 31
- 230000002123 temporal effect Effects 0.000 claims description 16
- 230000035515 penetration Effects 0.000 claims description 7
- 230000000694 effects Effects 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 5
- 238000011156 evaluation Methods 0.000 claims description 4
- 238000004321 preservation Methods 0.000 claims description 4
- 230000003137 locomotive effect Effects 0.000 claims description 3
- 238000012360 testing method Methods 0.000 claims description 3
- 238000012935 Averaging Methods 0.000 claims description 2
- 101001037160 Xenopus laevis Homeobox protein Hox-D1 Proteins 0.000 claims description 2
- 238000004364 calculation method Methods 0.000 claims description 2
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 230000003247 decreasing effect Effects 0.000 claims description 2
- 238000010586 diagram Methods 0.000 claims description 2
- 239000003344 environmental pollutant Substances 0.000 claims description 2
- 230000010429 evolutionary process Effects 0.000 claims description 2
- 231100000719 pollutant Toxicity 0.000 claims description 2
- 241001416181 Axis axis Species 0.000 claims 1
- 238000012546 transfer Methods 0.000 abstract description 10
- 238000005265 energy consumption Methods 0.000 description 5
- 230000007704 transition Effects 0.000 description 4
- 238000010521 absorption reaction Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000005431 greenhouse gas Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 239000003245 coal Substances 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000003912 environmental pollution Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000037361 pathway Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
-
- G—PHYSICS
- G08—SIGNALLING
- G08G—TRAFFIC CONTROL SYSTEMS
- G08G1/00—Traffic control systems for road vehicles
- G08G1/01—Detecting movement of traffic to be counted or controlled
Landscapes
- Business, Economics & Management (AREA)
- Engineering & Computer Science (AREA)
- Human Resources & Organizations (AREA)
- General Physics & Mathematics (AREA)
- Physics & Mathematics (AREA)
- Economics (AREA)
- Strategic Management (AREA)
- Development Economics (AREA)
- Game Theory and Decision Science (AREA)
- Entrepreneurship & Innovation (AREA)
- Marketing (AREA)
- Operations Research (AREA)
- Quality & Reliability (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- Theoretical Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
新能源汽车(主要包括插电式电动汽车(PEV)和混合动力电动汽车(HEV))正在逐步替代传统的内燃机汽车(ICEV)。由于ICEV、PEV和HEV在排放上的巨大差异，传统的交通环境影响评估软件正面临挑战。本发明使用R语言开发了一种考虑新能源汽车渗透的环境影响评估方法。首先，定义了三种类型车辆的不同排放率，将其表示为路段速度的函数。其次，采用吸收马尔可夫链模型对这三种车辆类型之间的转移过程进行建模，预测其市场占有率的变化。然后，采用带反馈的四阶段模型预测道路的交通流量和行驶时间。最后，根据各路段的平均行驶速度、路段长度、车辆的市场占有率水平进行环境影响评估。本软件对于评估新能源汽车的渗透对环境的影响十分有用。
Description
技术领域：
本发明采用R语言开发了一种考虑新能源汽车渗透的环境影响评估软件，属于交通工程技术领域。
背景技术：
交通领域的温室气体排放占全球总量的很大一部分，这使传统内燃机汽车(ICEV)长期受到人们的批评。最近，新能源汽车，包括插电式电动汽车(PEV)和混合动力电动汽车(HEV)，在交通中的使用对各个领域产生了重大影响。毫无疑问，新能源汽车是一个很好的环境解决方案。然而，传统的交通规划忽略了ICEV、PEV和HEV之间的巨大差异[1]。为了研究新能源汽车作为可持续交通方式的潜在影响，这项新技术必须被整合到交通网络中进行环境影响评价。
尽管传统的交通环境影响评估软件得到了广泛研究，但随着新能源汽车的渗透，其环境影响评估还未进行深入研究。随着新能源汽车成为无处不在的交通工具，Duell等[2]将PEV的能源消耗纳入网络设计决策中。他们设计了一个多目标规划问题，以最小化系统能源消耗和系统总出行时间。Gardner等[3]设计了一个双层框架，在出行需求变化的情况下，评估PEV在环境污染和能源消耗方面的影响。Jiang和Xie[4]为带有交通方式和路线选择的混合网络均衡问题提出了一种凸优化模型。其中，内燃机汽车和新能源汽车在行驶里程和出行成本构成方面有所区别。He等[5]发现由于行驶里程有限、充电站稀缺以及可能较长的电池充电或替换时间，PEV驾驶员的路线选择行为与传统ICEV驾驶员的路线选择行为不同。Papargyri等[6]利用仿真程序分析了电动汽车在未来十年中对温室气体减排的潜在贡献。Krause等[7]建立了贝叶斯信念网络来评估2030年德国轿车的可能特征，包括不同车辆类型的市场份额(包括ICEV、HEV和PEV)、二氧化碳排放量和用户成本。Javid和Nejat[8]评估了车辆车型的组成和排放水平的变化。但是，排放量不是在网络层次上进行评估。Ma等[9]定义了新能源汽车的环境成本并提出了随机用户均衡模型来描述驾驶员的路线选择行为。他们还展示了将新能源汽车引入交通网络中的环境效益，并揭示了新能源汽车数量与整个交通网络的环境成本之间的关系。Duel等[10]介绍了将出行者的行为和机动车的能源消耗纳入评估过程来促进新能源汽车发展的方法。
参考文献：
[1]Mitropoulos L K，Prevedouros P D.Incorporating sustainabilityassessment in transportation planning：an urban transportation vehicle-basedapproach[J].Transportation Planning and Technology，2016，39(5)：439-463.
[2]Duell M，Gardner L，Waller S T.Multiobjective Traffic Network DesignAccounting for Plug-in Electric Vehicle Energy Consumption[C].Proceedings ofProceedings of the 92nd Annual Meeting of the Transportation Research Board，Washington，D.C.，2013.
[3]Gardner L M，Duell M，Waller S T.A framework for evaluating the roleof electric vehicles in transportation network infrastructure under traveldemand variability[J].Transportation Research Part A：Policy and Practice，2013，49：76-90.
[4]Jiang N，Xie C.Computing and Analyzing Mixed Equilibrium NetworkFlows with Gasoline and Electric Vehicles[J].Comput.-Aided CivilInfrastruct.Eng.，2014，29(8)：626-641.
[5]He F，Yin Y，Lawphongpanich S.Network equilibrium models withbattery electric vehicles[J].Transportation Research Part B：Methodological，2014，67：306-319.
[6]Papargyri E，Kanaroglou P S，Photis Y N.Electric vehicles andtraffic related pollution reduction：a simulation model for Hamilton，Ontario，Canada[J].Glob.Nest.J.，2014，16(4)：753-761.
[7]Krause J，Small M J，Haas A，Jaeger C C.An expert-based bayesianassessment of 2030 German new vehicle CO2 emissions and related costs[J].Transp.Policy，2016，52：197-208.
[8]Javid R J，Nejat A.A comprehensive model of regional electricvehicle adoption and penetration[J].Transp.Policy，2017，54：30-42.
[9]Ma J，Cheng L，Li D W，Tu Q.Stochastic Electric Vehicle NetworkConsidering Environmental Costs[J].Sustainability，2018，10(8).
[10]Duell M，Gardner L M，Waller S T.Policy implications ofincorporating distance constrained electric vehicles into the traffic networkdesign problem[J].Transportation Letters-the International Journal ofTransportation Research，2018，10(3)：144-158.
发明内容：
技术问题：以往的软件至少存在三个问题。首先，随着电动汽车在交通网络中的渗透，其环境影响评估方法尚未被研究。其次，没有讨论不同车型车辆市场渗透水平的动态变化。本发明中，吸收马尔可夫链模型用于预测市场结构的变化，因为PEV被认为是不会转移到其他车辆类型的吸收状态。最后，在以往的环境影响评估中，只有交通分配用于交通规划，这是非常局限的。因此，本发明提出了带有反馈的四阶段模型，以实现交通系统均衡，充分反映出行者的决策行为。
技术方案：本发明采用R语言开发了一种考虑新能源汽车渗透的环境影响评估软件，该发明具体包括以下步骤：
(一)总体技术路线
步骤一：定义传统内燃机车(ICEV)、插电式电动汽车(PEV)和混合动力电动汽车(HEV)的排放率，将其表示为行驶速度的函数；
步骤二：使用吸收马尔可夫链模型表示ICEV、PEV和HEV的市场占有率演化过程，其中PEV是吸收状态，表示消费者选择PEV后不再会选择其他类型的汽车；
步骤三：通过交通产生、交通分布、交通方式划分和交通流分配的反馈迭代达到交通系统平衡，计算出平衡状态时的路段交通流量和通行时间；
步骤四：使用各路段的长度和行驶时间计算各路段的平均速度；
步骤五：评估交通系统的环境影响。
(二)步骤一的具体计算过程
对环境影响的分析需要ICEV、PEV和HEV的排放率。ICEV的排放与许多因素有关，包括交通方式、发动机类型、驾驶习惯、车辆速度等。但许多研究人员研究排放与平均车速之间的关系，因为其他因素是多样的，并且难以测量。这里，排放与速度的关系是从MOVES2010a(机动车辆排放模拟器)中得到的。每辆车每英里产生的排放(g/mi)通常拟合为路段a上平均行驶速度sa的幂函数。请注意，排放率随着速度的增加而降低。它们表示如下：
CO2 ICEV(sa)＝3158sa -0.56 (1)
VOCICEV(sa)＝1.3647sa -0.679 (2)
NOx ICEV(sa)＝2.5376sa -0.42 (3)
此外，路段a的平均速度sa定义为每单位时间所行驶的距离。用公式表示为
其中，la以英里为单位，是路段a的长度；ta(va)以小时为单位，是路段出行时间，是以pcu/h为单位的路段流量va的函数，以说明拥堵效应。因此，路段a上每个ICEV产生的排放(以克为单位)可以表示为路段流量va的函数，如下所示：
虽然PEV不直接产生尾气排放，但它们仍然间接地通过发电站污染大气，特别是燃煤电站。为了说明PEV对环境的影响，利用车辆消耗的总能源量乘以平均电站排放率(g/kWh)。计算结果是为了提供电动车辆消耗的能源而产生的总排放量。本方法使用从特斯拉汽车公司获得的数据，将每个PEV的每英里能源消耗率(kWh/mi)拟合为路段a上平均行驶速度sa的多项式函数：
ECPEV(sa)＝1.79e-8sa 4-4.073e-6sa 3+3.654e-4sa 2-0.0109sa+0.2372 (8)
同样，路段平均速度sa可以由公式(4)替换以考虑拥塞效应。因此，路段a上每个PEV的能源消耗(kWh)可以表示为路段流量va的函数，如下所示：
本发明使用北美环境合作委员会在分析中公布的CO2、NOx和SO2的平均排放率，分别为893g/kWh、1.66g/kWh、3.79g/kWh[3]。最后，路段a上每个PEV以克为单位产生的排放可以表示为路段流量va的函数，如下所示：
CO2 PEV(va)＝893·ECPEV(va) (10)
NOx PEV(va)＝1.66·ECPEV(va) (11)
SO2 PEV(va)＝3.79·ECPEV(va) (12)
HEV是一种将传统ICEV系统与PEV系统相结合的混合动力车辆。现代的HEV利用效率提升技术，使其比同等规模的ICEV产生更少的排放。由于HEV是ICEV和PEV的组合，因此为了简单起见，使用一个因子α来整合它们的排放。因此，路段a上每个HEV产生的以克为单位的排放也可以表示为路段流量va的函数，如下所示：
CO2 HEV(va)＝α·[CO2 ICEV(va)+CO2 PEV(va)] (13)
NOx HEV(va)＝α·[NOx ICEV(va)+NOx PEV(va)] (14)
VOCHEV(va)＝α·VOCICEV(va) (15)
SO2 HEV(va)＝α·SO2 PEV(va) (16)
(三)步骤二的具体计算过程
从长远来看，车主可以将车辆类型从一种替换为另一种。假设PEV是一种吸收类型，这意味着如果车主选择了PEV，他/她将不会转移到其他类型。三种车辆类型的转移过程如图1所示。因此，三种车型的市场渗透水平演化可以用吸收马尔可夫链(AMC)模型来表示。通过研究不同状态的初始概率和状态之间的转移概率来确定不同状态的市场渗透率，从而达到预测未来的目的。
考虑一组状态S＝{s1，s2，...，sr}，演变过程从这些状态之一开始，并连续从一个状态转移到另一个状态。每次转移都称为一步。如果当前处于状态si，则在下一步骤，以转移概率pij移动到状态sj，并且该概率不取决于当前状态之前的状态。该过程也可以保持其所处的状态，并且这以概率pii发生。设P＝[pij]为马尔可夫链的转移矩阵。矩阵Pn的第ij项给出了马尔可夫链从状态si开始，在n步之后处于状态sj的概率。注意到，Pn的每个行向量的每个元素非负且和为1。设u是表示初始分布的概率向量。那么，在n步之后链处于状态si的概率是下面向量中的第i项
u(n)＝uPn (17)
对于AMC模型，对状态重新编号，以便使过渡状态排在前面。如果前q个状态是过渡状态且最后的r个状态是吸收状态，则转移矩阵将具有以下标准形式
其中I是r行r列的单位矩阵，0是r行q列的零矩阵，R是q行r列的非零矩阵，Q是q行q列的矩阵。对于吸收马尔可夫链P，矩阵N＝(U-Q)-1为矩阵P的基本矩阵，其中U是q行q列的单位矩阵。
在吸收马尔可夫链中，过程被吸收的概率是1。假设链在状态si开始，设ti为链变为吸收状态之前的预期转移步骤数，并且令t为第i个元素为ti的列向量。那么
t＝Nc， (18)
其中，c是列向量，其所有的元素都是1。
(四)步骤三的具体计算过程
步骤三是带反馈的四阶段模型，算法流程图如图2所示，具体的计算过程如下：
步骤1：输入交通网络信息，包括网络结构、路段长度、路段通行能力等·；
步骤2：给定交通需求Oi，通过平均分布得到初始化交通分布矩阵设置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上的交通流量va和起点i与目的地j之间的出行时间
(五)程序设计
#考虑新能源汽车的渗透对环境影响的评估软件
#步骤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″，″Road capacity″，″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.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
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″，″Level of 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_VOC＝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()#查看输出文件的保存地址。
有益效果：目前，新能源汽车的引入是缓解交通污染的有效途径。然而，这对传统的环境影响评估提出了巨大的挑战，因为新能源汽车具有与传统内燃机车不同的特征。在交通网络层面对新能源汽车的渗透进行环境影响评估对于短期和长期的新能源汽车政策制定都非常有用。此外，本软件具有以下技术特点：1)开源免费：R语言是开源软件，而交通规划与管理软件多是国外软件，价格昂贵；2)易于操作：国外商业软件普遍界面复杂、应用繁琐，难以为一般的工程师和规划师掌握，而本软件易于操作；3)矩阵算法快捷高效：国外商业软件多采用C语言开发，在进行算法实现时比较低效，采用矩阵语言大大加快了计算过程。
附图说明：
图1是三种车辆类型间的转移过程
图2是算法流程图
图3是Nguyen-Dupuis测试网络
图4是模拟的三种车型的转移过程
图5是各类型车辆市场占有率水平的动态变化趋势
具体实施方式：
如图3所示的Nguyen-Dupuis交通网络广泛用于交通研究以验证各种方法。路段参数包括自由流动出行时间，路段通行能力和路段长度，如表1所示。
表1.Nguyen-Dupuis网络的路段参数
Nguyen-Dupuis网络中有两个起点区域1和4以及两个目的地区域2和3。假设在峰值时间内起点1和4的交通发生分别为1800pcu/h和1200pcu/h。也就是说，O1＝1800pcu/h，O4＝1200pcu/h。在交通分布步骤中，众所周知，传统的重力模型中不包括许多具有显著解释力的关键变量。其中，最有影响力的是出行者目的地偏好。例如，出行者通常喜欢传统的目的地而不是新开发的地区。因此，目的地选择采用具有出行者偏好的多项式logit模型。尽管存在多种解释变量，为简便起见，反馈过程中的目的地选择模型被简化为
其中βj是出行者对目的地j上的偏好，βt是O-D对ij之间的径路出行时间系数。可以使用经验数据校准βj和βt的值。在这里我们设β2＝0、β3＝1和βt＝-0.1。也就是说，目的地区域2上的出行者偏好是0，而目的地区域3是1，这意味着出行者通常更喜欢目的地3。出行时间系数是-0.1，这意味着出行时间是负效用。
在交通流分配步骤中，本发明采用经典的用户均衡方法，该方法将路段性能函数应用于均衡状态中。常用的路段性能函数由美国公路局(BPR)开发，其表示如下：
其中ta(va)是具有交通流量va的给定路段a的阻抗函数；ca是路段通行能力；是路段a的自由流阻抗；α和β是可以凭经验校准的延迟系数。对于α和β，传统的BPR值分别为0.15和4.0，这也用于我们的仿真研究中。因此，我们可以通过Frank-Wolfe算法获得路段交通流量va。
收敛标准RRSE设定为0.01，即ε＝0.01。通过使用上述参数，可以收敛得到单个稳定的解，这个解具有一致的出行时间/成本和交通分布矩阵。此外，可以使用Dijkstra算法计算i和j之间的最短路径出行时间。给定交通网络下的路段性能如表2所示。
表2.给定交通网络下的路段性能
三种车型的市场渗透水平是动态的，随着时间的推移而变化。车主可以继续使用一种类型的车辆或转移到其他类型的车辆。注意到PEV是一种吸收状态。这可以用吸收马尔可夫链模型表示。转移概率可能受到交通政策的影响，例如新能源汽车的补贴。依次对于ICEV、HEV和PEV，假设一步转移矩阵P是
吸收马尔可夫链模型中假设的一步转移过程如图4所示。因此，功能矩阵N为
根据公式(18)，在链变为吸收状态前，期望的转移步骤次数(即：年)为
这意味着目前的ICEV所有者预计14.4年后变为吸收状态，目前的HEV所有者预计13.3年后变为吸收状态。请注意，这里的期望值并不意味着大约十五年后将不存在ICEV和HEV。
此外，假设目前的市场份额是u＝(0.7，0.2，0.1)，这意味着ICEV目前占据了大部分的汽车市场。根据公式(17)，图5中显示了未来20年市场渗透水平的动态变化。
请注意，假设HEV和PEV出行者的行为与传统的ICEV出行者相同。目前，他们在行为方面与ICEV出行者没有区别。在目的地选择、路线选择等方面，它们与传统车辆的出行者的行为相同。因此，假设车辆在每个路段上的比例与其市场渗透水平一致。然而，三种车辆类型，ICEV、PEV和HEV，通过其先前定义的排放率来区分。虽然据报道，HEV可以减少高达90％的烟雾污染物排放，并将二氧化碳排放量减少一半，但这里的因子α设定为0.6。此外，通过交通系统均衡能够得到如表2所示的路段交通量和路段出行时间。因此，可以在交通网络中随时间评估环境影响，如表3所示。
表3.随时间的战略环境影响
很明显，随着传统的ICEV逐渐被放弃，CO2、NOx和VOC的排放不断下降。它们在最后一年的排放量分别占基准年的43.8％，53.9％和19.8％。然而，随着电动汽车成为主流，SO2的排放量在增加。在最后一年SO2的排放量是基准年的3.913倍。虽然最初几年发生了巨大的变化，但在最后几年环境影响的变化开始变得稳定，这被称为稳定状态。
Claims (3)
1.本发明采用R语言开发了一种考虑新能源汽车渗透的环境影响评估软件，该软件的算法框架具体包括以下步骤：
步骤一：定义传统内燃机车(ICEV)、插电式电动汽车(PEV)和混合动力电动汽车(HEV)的排放率，将其表示为行驶速度的函数；
步骤二：使用吸收马尔可夫链模型表示ICEV、PEV和HEV的市场占有率演化过程，其中PEV是吸收状态，表示消费者选择PEV后不再会选择其他类型的汽车；
步骤三：通过交通产生、交通分布、交通方式划分和交通流分配的反馈迭代达到交通系统平衡，计算出平衡状态时的路段交通流量和通行时间；
步骤四：使用各路段的长度和行驶时间计算各路段的平均速度；
步骤五：评估交通系统的环境影响。
2.对于权利1中的步骤三，采用如下的具体计算过程：
步骤1：输入交通网络信息，包括网络结构、路段长度、路段通行能力等；
步骤2：给定交通需求Oi，通过平均分布得到初始化交通分布矩阵设置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上的交通流量va和起点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()#查看输出文件的保存地址。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910771372.7A CN110489871A (zh) | 2019-08-20 | 2019-08-20 | 考虑新能源汽车渗透的环境影响评估软件 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910771372.7A CN110489871A (zh) | 2019-08-20 | 2019-08-20 | 考虑新能源汽车渗透的环境影响评估软件 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN110489871A true CN110489871A (zh) | 2019-11-22 |
Family
ID=68551606
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910771372.7A Pending CN110489871A (zh) | 2019-08-20 | 2019-08-20 | 考虑新能源汽车渗透的环境影响评估软件 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110489871A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111275965A (zh) * | 2020-01-20 | 2020-06-12 | 交通运输部科学研究院 | 一种基于互联网大数据的实时交通仿真分析系统与方法 |
| CN112613699A (zh) * | 2020-12-01 | 2021-04-06 | 吉林大学 | 基于汽车大数据的初始碳排放权分配方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102054222A (zh) * | 2010-12-10 | 2011-05-11 | 吉林大学 | 基于居民出行分析的城市机动车排放量化方法 |
| CN104408924A (zh) * | 2014-12-04 | 2015-03-11 | 深圳北航新兴产业技术研究院 | 一种基于耦合隐马尔可夫模型的城市道路异常交通流检测方法 |
| CN108647475A (zh) * | 2018-05-25 | 2018-10-12 | 东南大学 | 基于负载均衡的城市离散交通网络设计r语言实现方法 |
| CN108804801A (zh) * | 2018-05-25 | 2018-11-13 | 东南大学 | 基于目标配流的城市离散交通网络设计r语言实现方法 |
| CN108876042A (zh) * | 2018-06-08 | 2018-11-23 | 东南大学 | 新型交通分布与交通流分配组合模型的r语言处理方法 |
| CN109409950A (zh) * | 2018-10-12 | 2019-03-01 | 东南大学 | 基于交通安全的最优收费定价软件 |
| CN110472354A (zh) * | 2019-08-20 | 2019-11-19 | 东南大学 | 一种新能源汽车渗透的环境影响评估方法 |
-
2019
- 2019-08-20 CN CN201910771372.7A patent/CN110489871A/zh active Pending
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102054222A (zh) * | 2010-12-10 | 2011-05-11 | 吉林大学 | 基于居民出行分析的城市机动车排放量化方法 |
| CN104408924A (zh) * | 2014-12-04 | 2015-03-11 | 深圳北航新兴产业技术研究院 | 一种基于耦合隐马尔可夫模型的城市道路异常交通流检测方法 |
| CN108647475A (zh) * | 2018-05-25 | 2018-10-12 | 东南大学 | 基于负载均衡的城市离散交通网络设计r语言实现方法 |
| CN108804801A (zh) * | 2018-05-25 | 2018-11-13 | 东南大学 | 基于目标配流的城市离散交通网络设计r语言实现方法 |
| CN108876042A (zh) * | 2018-06-08 | 2018-11-23 | 东南大学 | 新型交通分布与交通流分配组合模型的r语言处理方法 |
| CN109409950A (zh) * | 2018-10-12 | 2019-03-01 | 东南大学 | 基于交通安全的最优收费定价软件 |
| CN110472354A (zh) * | 2019-08-20 | 2019-11-19 | 东南大学 | 一种新能源汽车渗透的环境影响评估方法 |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111275965A (zh) * | 2020-01-20 | 2020-06-12 | 交通运输部科学研究院 | 一种基于互联网大数据的实时交通仿真分析系统与方法 |
| CN111275965B (zh) * | 2020-01-20 | 2021-02-05 | 交通运输部科学研究院 | 一种基于互联网大数据的实时交通仿真分析系统与方法 |
| CN112613699A (zh) * | 2020-12-01 | 2021-04-06 | 吉林大学 | 基于汽车大数据的初始碳排放权分配方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Thibault et al. | A unified approach for electric vehicles range maximization via eco-routing, eco-driving, and energy consumption prediction | |
| Zhang et al. | Multitype recharge facility location for electric vehicles | |
| Wen et al. | System dynamics modeling and policy simulation for urban traffic: a case study in Beijing | |
| CN110472354A (zh) | 一种新能源汽车渗透的环境影响评估方法 | |
| CN105760960A (zh) | 基于轨道交通的停车换乘设施最优选址和容量确定方法 | |
| CN110472353B (zh) | 基于用户效用最大化的交通网络设计方法 | |
| CN101295326A (zh) | 基于gps数据生成od矩阵的方法及其交通仿真的方法 | |
| CN107798867A (zh) | 一种基于电动汽车和内燃机车混合车流的拥堵交通网络均衡方法 | |
| CN108804801B (zh) | 基于目标配流的城市离散交通网络设计r语言实现方法 | |
| CN110489871A (zh) | 考虑新能源汽车渗透的环境影响评估软件 | |
| Wang et al. | Optimal recharging strategies for electric vehicle fleets with duration constraints | |
| Lovrić et al. | A conceptual design for a national transport model with cross-sectoral interdependencies | |
| Boriboonsomsin et al. | Examination of attributes and value of ecologically friendly route choices | |
| Scarinci et al. | Electrification of urban mobility: The case of catenary-free buses | |
| Madadi et al. | Multi-stage optimal design of road networks for automated vehicles with elastic multi-class demand | |
| Huang et al. | Eco-mobility-on-demand fleet control with ride-sharing | |
| Zakharov et al. | Transit network design for green vehicles routing | |
| Sun et al. | A prediction-evaluation method for road network energy consumption: Fusion of vehicle energy flow principle and Two-Fluid theory | |
| CN110457012A (zh) | 用于可持续交通网络设计的多属性决策软件 | |
| CN108647835B (zh) | 基于设计速度的城市离散交通网络设计r语言实现方法 | |
| CN116358593A (zh) | 考虑非线性能耗的电动车辆路径规划方法、装置和设备 | |
| Amirjamshidi | Assessment of commercial vehicle emissions and vehicle routing of fleets using simulated driving cycles | |
| Valencia et al. | Is switching propulsion technologies the path to sustainable land transport? decarbonizing Bogotá | |
| Batista et al. | On the characterization of eco-friendly paths for regional networks | |
| Xin et al. | A battery electric vehicle transportation network design model with bounded rational travelers |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20191122 |
|
| WD01 | Invention patent application deemed withdrawn after publication |