由于电动汽车拥有减少温室气体排放、促进碳达峰、缓解化石燃料短缺危机等优点，世界各国已经颁布多项法令来加速电动汽车的普及。电动汽车的渗透速度在中国尤其快，其销量和保有量目前均位列世界第一^[1]。在此背景下，电动汽车在路网中的滲透率不断提高，使得路网的电气化程度日益增加。电动汽车电池容量有限、需要一定充电时间、充电站和充电桩数量有限等特征为路网的应急管理带来了挑战。路网是保障社会正常运转的关键基础设施，一旦遭受破坏，会给社会和经济造成巨大的损失。特别是，近年来极端天气在全球范围内频繁发生，给现有关键基础设施的平稳运行造成了巨大的威胁。比如，2021年7月，河南省部分地区发生极端暴雨和特大洪灾，导致交通系统遭受到严重损毁进而引发大规模交通瘫痪，造成了大规模死伤和严重的经济损失。这类例子都在提示，加强路网的脆弱性分析，并将电动汽车考虑在内，已经成为保障路网正常运转的迫切需要。

在过去的几十年里，国内外学者持续对交通系统的脆弱性分析问题展开了大量研究^[2-5]，然而电气化交通路网的脆弱性研究几乎无人涉及。要分析电气化路网的脆弱性，无法避免的一个问题是如何对要研究的系统建立数学模型。然而，目前尚无成熟的方法建立电气化路网系统模型，而传统的交通系统建模方法^[6-7]大致分为2类：

第1类是基于交通网络的拓扑结构进行静态的脆弱性分析。在此类工作中，交通网络被简化为抽象的节点、边的拓扑结构。大多数工作建立在复杂网络理论与静态交通流分配理论的基础上^[8-9]，通过研究路网的拓扑结构、比较事故发生前后路网的平衡运行状态，来建立脆弱性评价指标以评估路网的关键环节。例如，文[10]分析发现，道路车道数是影响链路关键性的最重要参数，其次是道路长度和距城市中心的距离。文[5]探讨了不同道路基础设施属性对网络脆弱性的影响。文[11]采用静态脆弱性评估原理，在不考虑连锁失效的前提下计算结构脆弱性指数；同时，还采用了动态脆弱性评估原理，在考虑交通拥堵传播效应的前提下计算运行状态脆弱性指数；结果表明，路网的脆弱性受网络拓扑结构和网络运行状态2个因素的影响。文[12]基于复杂网络理论的拓扑规则，对成渝地区的城际铁路网络进行建模，通过模拟随机和蓄意攻击，利用PageRank算法识别网络中的关键节点；实验结果表明，重庆、成都、自贡等为关键节点，网络在蓄意攻击下更加脆弱。这类方法的优点是需要的数据较少，通用性较好，并且有着严格的数学推导作为支撑；其弱点是由于抽象后的网络过于简单，可能使得该类方法无法刻画系统在破坏后的某一细化的动态响应。比如在交通系统中，可能无法描述破坏的持续时间、替代路线上流量的增加等。

第2类是基于特定的系统模型进行脆弱性分析。此类模型通常考虑了该特定系统的物理约束和时间因素，因此可以更加合理地描述系统在破坏后的响应。比如，文[13]提出了一种铁路网络脆弱性模型，通过找到对乘客和火车造成最不利后果的铁路链接组合来评估系统的脆弱性。该方法可以提供链接、相应的客流、火车路线和时刻表的关键组合。该实验结果表明，关键铁路链接与交通需求高度相关，与网络拓扑的静态特征关系不大。文[14]提出了一个双层规划模型来设计城市道路网络。其上层模型在保证可靠性的约束条件下，优化路网脆弱性指标；其下层模型为基于后悔理论和效用理论的随机用户均衡模型。实验结果表明，设计路网时需要协同考虑其脆弱性和可靠性。并且，考虑用户出行决策差异可以提高路网设计方案的准确性。文[15]开发了一种日常动态网络脆弱性分析方法，该方法允许在考虑客观旅行时间不确定性和主观感知时间误差的动态模型的基础上，考虑网络脆弱性。文[16]将脆弱性评价与交通流演化过程相结合，从时间维度研究路网脆弱性。基于动态用户最优原则，建立了基于动态路径选择的改进动态流量分配模型，提出了脆弱性评估的累计时变指标。研究结果表明，更准确地识别和评估脆弱节点和脆弱环节，有助于增强路网抗干扰能力。文[17]以广义出行成本为基础，构建城市路网脆弱性评价指标，从用户广义出行成本的角度对城市路网进行脆弱性评价。此类方法的优点和缺点与第1类方法互补。由于第2类方法需要更多数据来对具体的系统特征进行建模，如交通系统中的用户行为，因此一方面可以相对合理地模拟系统破坏后的具体响应，另一方面由于刻画这些细节需要花费更多的计算资源，导致在大规模网络中的应用受到限制。

近年来关键基础设施之间的相互依赖受到了越来越多的关注^[18-19]。与过去只考虑对交通系统建模不同，同时更加全面地考虑交通与电力^[20]、燃气^[21]、供水^[22]网络等其他关键基础设施之间的相互作用，已经逐渐成为未来的热门研究趋势。

电动汽车的快速普及和充电设施的大规模部署使得交通路网的电气化程度不断提高。然而，路网很容易受到几个关键位置发生的故障的影响，进而可能会引发一系列级联故障，并最终导致大范围交通瘫痪。在混合电动汽车和非电动汽车流的背景下，由于电动汽车里程有限、需要一定充电时间等特征，负面影响可能进一步扩散和叠加，从而对交通路网的正常运行造成更加严重的后果。对脆弱道路的保护或备份可以降低电气化路网面对自然灾害和恶意攻击等的脆弱性，从而降低级联失效和破坏事件进一步扩散的风险。因此，对电气化路网进行脆弱性评估并识别系统中的关键道路非常重要。

针对现有文献中尚未有对电气化路网的脆弱性进行研究的现状，本文提出一个双层攻击者-防御者模型来研究电气化路网的脆弱性，并把该双层混合整数线性规划问题转化为一个单层的混合整数线性规划问题，使其求解效率更高。攻击者的目标是使系统性能下降最大化，防御者的目标是使系统的行驶时间最小化。系统最优动态混合交通流模型被用来模拟电气化路网中混合电动汽车和非电动汽车的交通流分布。通过求解该模型，可得到系统中的最脆弱道路集合以及其被破坏后对应情形下系统性能的变化。本文验证了在评估路网脆弱性时，把电动汽车考虑在内的必要性，可以为改善电气化路网的脆弱性提供决策依据和理论支撑。

1 模型建立

本章建立一种考虑电动汽车的链路传输模型、一种双层攻击者-防御者模型, 并提出2个用于评估电气化路网道路脆弱性的动态指标。

1.1 考虑电动汽车的链路传输模型

本节提出了一种考虑电动汽车和快速充电站重要特征的链路传输模型。该模型基于链路传输模型(link transmission model，LTM)。然后，基于该考虑电动汽车的链路传输模型，提出双层攻击者-防御者模型。

本文基于3个基本假设：1) 电动汽车的电量消耗与电动汽车行驶距离呈线性相关。2) 电动汽车在充电站所充的电量与其充电时间呈线性相关。3) 电动汽车上譬如电灯、空调等设备消耗的电量被忽略。

电气化路网被记作$G(\mathcal{N}, \mathcal{A})$, 其中$\mathcal{N}$和$\mathcal{A}$分别为节点与链路的集合。在电气化路网中的链路可被分为源链路$\mathcal{A}_{\mathrm{R}}$、终点链路$\mathcal{A}_{\mathrm{S}}$、普通链路$\mathcal{A}_{\mathrm{G}}$以及充电链路$\mathcal{A}_{\mathrm{C}}$。虚拟的充电链路用于描述实际电气化路网的快速充电站, 其首尾相连。如图 1所示, 具有差异化充电速度的充电桩可以通过多条充电链路表示。节点分为2种类型: 源和终点节点$\mathcal{N}_{\mathrm{SR}}$以及普通节点$\mathcal{N}_{\mathrm{G}}$。在电气化路网中，每个源点节点和终点节点仅连接一个源链路和终点链路。所有充电、源和终点链路都是长度为0的虚拟链路，因此在这些虚拟链路上不计算行驶时间。所有源和终点链路都具有无限的流出、流入和承载车辆的容量，以保证它们不会成为电气化路网模型中流量的瓶颈。对于系统最优动态交通流分配(system optimal dynamic traffic assignment，SODTA)问题，假定所有终点链路的流出容量为0，类似于文[23-25]，即所有车辆在抵达时都会被收集。时间被离散化为一组有限的周期的集合$\mathcal{T}=\{1, 2$, $\left.\cdots, t, \cdots, T_{\text {max }}\right\}$。$T_{\text {max }}=H / \delta$, 其中$\delta$是周期长度, $H$是对系流观测的总时长。周期长度$\delta$应等于或小于最小的链路行驶时间，以使车辆至少需要一个时间单位来穿行链路^[26]。

为了近似道路的宏观特性，LTM中定义了三角形的基本型^[26]。该图由3个参数定义: 堵塞密度$\left(k_{\mathrm{jam}}\right)$、最大流量$\left(q_{\text {max }}\right)$和固定自由流速$\left(v_{\mathrm{f}}\right)$。反向激波速度$w=q_{\text {max }} \cdot v_{\mathrm{f}} /\left(q_{\text {max }}-k_{\text {jam }} \cdot v_{\mathrm{f}}\right)$。

假设存在一辆类型为$c$的电动汽车, 其电池容量为$B_{c}(\mathrm{~kW} \cdot \mathrm{h})$, 能量消耗速率为$\eta(\mathrm{kW} \cdot \mathrm{h} /\mathrm{mile}$, 1 mile $=1.609344 \mathrm{~km})$, 最大行驶里程为$L_{c}^{\max }$ (mile), $L_{c}^{\text {max }}=B_{c} / \eta$。本文将其行驶里程离散化为整数倍的单位能级(energy levels, ELs)。当电动汽车的电量为满时, 它具有最大能级$E_{c}=L_{c}^{\max } /$ $\left(\delta \cdot v_{\mathrm{f}}\right)$。电动汽车行驶$\delta \cdot v_{\mathrm{f}}$需消耗$1 \mathrm{Els}$。本文假设存在着$C$个电动汽车类型, 电动汽车类型的集合记为$\mathcal{C}=\{1, 2, \cdots, c, \cdots, C\}$。电动汽车等级的集合$\mathit{\Xi }=\left\{\mathcal{E}_{1}, \mathcal{E}_{2}, \cdots, \mathcal{E}_{c}\right\}$。每个集合$\mathit{\Xi }$的元素$\mathcal{E}_{i}$本身为一个集合, 其中包含着$c$类电动汽车可能具有的能级, 被记作$\mathcal{E}_{c}=\left\{1, 2, \cdots, e, \cdots, E_{c}\right\}$。在链路传输模型中, 交通流的动态演变是通过计算每个周期$t$各链路累计进入和离开的车辆数目获得的。

在本文提出的电气化交通路网模型中，目标是最小化所有车辆的总行驶时间。总行驶时间包括整个时间范围内所有链路上所有车辆的总停留时间和所有电动汽车的总充电时间。在正常情况下，电气化交通路网的混合交通流分配模型的目标函数为

$ \begin{aligned} & \min \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{t \in \mathcal{T}} \sum\limits_{a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}, \mathcal{A}_{\mathrm{S}}\right\}} \delta\left[\mathrm{UG}_{a}^{s}(t)-\mathrm{VG}_{a}^{s}(t)\right]+ \\ & \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{t \in \mathcal{T}} \sum\limits_{a \in \mathcal{A} \backslash \mathcal{A}_{\mathrm{S}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}} \delta\left[\mathrm{UE}_{a, c}^{s, e}(t)-\mathrm{VE}_{a, c}^{s, e}(t)\right] . \end{aligned} $

(1)

其中: $\mathrm{UE}_{a, c}^{s, e}(t) 、\mathrm{VE}_{a, c}^{s, e}(t)$代表终点节点为$s$的$c$类型具有能级$e$的电动汽车在周期$t$结束前进入、离开链路$a$的累计数量; $\mathrm{UG}_{a}^{s}(t) 、\mathrm{VG}_{a}^{s}(t)$代表终点节点为$s$的燃油汽车$t$结束前进入、离开链路$a$的累计数量。

Newell的简化运动波理论^[27-28]在链路传输模型中被用于计算链路$a$的流出量$S_{a}(t)$和流入量$R_{a}(t)$ :

$ S_{a}(t)=\min \left\{U_{a}\left(t-v_{a}\right)-V_{a}(t-1), f_{a}^{\mathrm{O}}(t)\right\} . $

(2a)

$ \begin{gathered} R_{a}(t)=\min \left\{V_{a}\left(t-\beta_{a}\right)+\right. \\ \left.L_{a} \cdot k_{\text {jam }}-U_{a}(t-1), f_{a}^{\mathrm{I}}(t)\right\} . \end{gathered} $

(2b)

其中: $U_{a}(t) 、V_{a}(t)$代表着在周期$t$结束前进入、离开链路$a$的累计车辆数目; $f_{a}^{\mathrm{l}}(t)$和$f_{a}^{\mathrm{O}}(t)$分别代表着在周期$t$, 链路$a$最大允许的流入量和流出量, 可以通过计算每个周期链路$a$的$\delta \cdot q_{\max }$获得。$v_{a}$和$\beta_{a}$分别为以自由流速度和反向激波速度穿行链路$a$所消耗的时间, 可以通过$v_{a}=L_{a} /\left(\delta \cdot v_{\mathrm{f}}\right)$和$\beta_{a}=L_{a} /(\delta \cdot w)$来计算。$L_{a}$是链路$a$的长度。

链路a在周期t的流入量和流出量被链路对应的最大允许流入和流出量分别约束：

$ \begin{gathered} U_{a}(t)-U_{a}(t-1) \leqslant R_{a}(t), \\ \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t . \end{gathered} $

(3a)

$ \begin{gathered} V_{a}(t)-V_{a}(t-1) \leqslant S_{a}(t), \\ \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t . \end{gathered} $

(3b)

将式(2a)和(2b)代入到不等式组(3)中, 获得了线性的基于链路传输模型的交通流约束组:

$ V_{a}(t) \leqslant U_{a}\left(t-v_{a}\right), \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t. $

(4)

$ \begin{gathered} V_{a}(t)-V_{a}(t-1) \leqslant f_{a}^{\mathrm{O}}(t), \\ \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t . \end{gathered} $

(5)

$ \begin{gathered} U_{a}(t)-U_{a}(t-1) \leqslant f_{a}^{\mathrm{I}}(t), \\ \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t . \end{gathered} $

(6)

$ \begin{gathered} U_{a}(t)-V_{a}\left(t-\beta_{a}\right) \leqslant L_{a} k_{\mathrm{jam}}, \\ \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t . \end{gathered} $

(7)

在本文提出的模型中，考虑电动汽车和燃油汽车的入流量和出流量分别为：

$ \begin{gathered} U_{a}(t)=\sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \mathrm{UG}_{a}^{s}(t)+ \\ \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}} \operatorname{UE}_{a, c}^{s, e}(t), \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t . \end{gathered} $

(8a)

$ \begin{gathered} V_{a}(t)=\sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \mathrm{VG}_{a}^{s}(t)+ \\ \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}} \operatorname{VE}_{a, c}^{s, e}(t), \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t . \end{gathered} $

(8b)

将式(7)代入式(3)—(6)，得到了对电动汽车和燃油汽车的混合流的约束组：

$ \begin{gathered} \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}}\left[\mathrm{VG}_{a}^{s}(t)-\mathrm{VG}_{a}^{s}(t-1)\right]+ \\ \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}}\left[\mathrm{VE}_{a, c}^{s, e}(t)-\mathrm{VE}_{a, c}^{s, e}(t-1)\right] \leqslant \\ f_{a}^{\mathrm{O}}(t), \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t, s . \end{gathered} $

(9)

$ \begin{gathered} \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}}\left[\mathrm{UG}_{a}^{s}(t)-\mathrm{UG}_{a}^{s}(t-1)\right]+ \\ \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}}\left[\mathrm{UE}_{a}^{s}(t)-\mathrm{UE}_{a}^{s}(t-1)\right] \leqslant \\ f_{a}^{\mathrm{I}}(t), \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t, s . \\ \end{gathered} $

(10)

$ \begin{gathered} \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}}\left[\mathrm{UE}_{a, c}^{s, e}(t)-\mathrm{VE}_{a, c}^{s, e}\left(t-\beta_{a}\right)\right]+ \\ \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}}\left[\mathrm{UG}_{a}^{s}(t)-\mathrm{VG}_{a}^{s}\left(t-\beta_{a}\right)\right] \leqslant \\ L_{a} k_{\mathrm{jam}}, \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t, s . \end{gathered} $

(11)

对于燃油汽车，具有不同终点的燃油车流分流的流出量应被流入处的边界条件所约束；对于不同等级、不同能级、不同终点的电动汽车流分流的流出量，也应被流入处的边界条件所约束。因而得到：

$ \begin{gathered} \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \mathrm{VG}_{a}^{s}(t) \leqslant \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \mathrm{UG}_{a}^{s}\left(t-v_{a}\right), \\ \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t . \end{gathered} $

(12)

$ \begin{gathered} \mathrm{VE}_{a, c}^{s, e}(t) \leqslant \mathrm{UE}_{a, c}^{s, e+\rho_{a}}\left(t-v_{a}\right), \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \\ e \in \mathcal{E}_{c} \cap\left\{e \leqslant E_{c}-\rho_{a}\right\}, \forall s, c, t . \end{gathered} $

(13a)

$ \begin{gathered} \operatorname{VE}_{a, c}^{s, e}(t)=0, \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \\ e \in \mathcal{E}_{c} \cap\left\{e>E_{c}-\rho_{a}\right\}, \forall s, c, t . \end{gathered} $

(13b)

其中$\rho_{a}$是穿行链路$a$所消耗的能级, $\rho_{a}=L_{a} /(\delta \cdot$$\left.v_{\mathrm{f}}\right)$。式(13a) 保证了电动汽车的流出量应当小于等于在$v_{a}$前的流入量, 同时也保证了流出车流的能级得到更新。式(13b) 确保了所有电动汽车的能级都小于它们的最大能级。

本文通过让源链路累计流入车辆等于累计需求来满足交通需求：

$ \mathrm{UG}_{a}^{s}(t)=\mathrm{DG}_{a}^{s}(t), \forall a \in \mathcal{A}_{\mathrm{R}}, \forall s, t . $

(14a)

$ \begin{gathered} \quad \mathrm{UE}_{a, c}^{s, e}(t)=\mathrm{DE}_{a, c}^{s, e}(t), \\ \forall a \in \mathcal{A}_{\mathrm{R}}, e \in \mathcal{E}_{c}, \forall s, c, t . \end{gathered} $

(14b)

其中: $\mathrm{DG}_{a}^{s}(t) 、\mathrm{DE}_{a, c}^{s, e}(t)$代表在第$t$周期末尾前, 在源链路$a$入口的燃油汽车、电动汽车的累计交通需求。

为了使得交通流守恒，本文使普通节点处流出量等于其流入量：

$ \begin{gathered} \sum\limits_{a \in B(i)} \mathrm{VG}_{a}^{s}(t)=\sum\limits_{b \in A(i)} \mathrm{UG}_{a}^{s}(t), \\ \forall i \in \mathcal{N} \backslash\left\{\mathcal{N}_{\mathrm{SR}}\right\}, \forall s, t . \end{gathered} $

(15a)

$ \begin{gathered} \sum\limits_{a \in B(i)} \operatorname{VE}_{a, c}^{s, e}(t)=\sum\limits_{b \in A(i)} \mathrm{UE}_{a, c}^{s, e}(t), \\ \forall i \in \mathcal{N} \backslash\left\{\mathcal{N}_{\mathrm{SR}}\right\}, \forall e \in \mathcal{E}_{c}, \forall s, c, t . \end{gathered} $

(15b)

其中: $A(i) 、B(i)$代表尾、首节点为$i$的链路的集合。

对于电动汽车，充电链路a上的现有充电桩占有量应当小于等于该链路上的充电桩数目。

$ \begin{gathered} \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}}\left[\mathrm{UE}_{a, c}^{s, e}(t)-\mathrm{VE}_{a, c}^{s, e}(t)\right] \leqslant \\ \operatorname{NC}_{a}(t), \forall a \in \mathcal{A}_{\mathrm{C}}, \forall t . \end{gathered} $

(16)

其中$\mathrm{NC}_{a}(t)$是充电链路$a$上在周期$t$的实际充电桩数目。

式(17)用于更新充电链路上现有充电桩占用量以及电动汽车的能级，

$ \begin{gathered} \hat{x}_{a, c}^{s, e}(t)=x_{a, s}^{s, e}(t-1)+\left[\mathrm{UE}_{a, c}^{s, e}(t-1)-\right. \\ \left.\mathrm{UE}_{a, c}^{s, e}(t-2)\right]-\left[\operatorname{VE}_{a, c}^{s, e}(t-1)-\mathrm{VE}_{a, c}^{s, e}(t-2)\right], \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall e \in \mathcal{E}_{c}, \forall s, c, t . \end{gathered} $

(17)

其中: $\hat{x}_{a, c}^{s, e}(t)$和$x_{a, s}^{s, e}(t)$分别为在充电链路$a$更新能级之前和之后的电动汽车的数目。

基于获得的充电桩占有量，式(18)通过线性地更新电动汽车能级描述了电动汽车的充电过程：

$ x_{a, c}^{s, E_{c}}(t)=\sum\limits_{l=0}^{\alpha_{a}^{t}} \hat{x}_{a, c}^{s, E_{c}-l}(t), \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s, c, t . $

(18a)

$ \begin{gathered} x_{a, c}^{s, e}(t)=\hat{x}_{a, c}^{s, e-\alpha_{a}^{t}}(t), \forall a \in \mathcal{A}_{\mathrm{C}}, \\ \forall e \in\left\{\alpha_{a}^{t} \leqslant e <E_{c}\right\}, \forall s, c, t . \end{gathered} $

(18b)

$ \begin{gathered} x_{a, c}^{s, e}(t)=0, \forall a \in \mathcal{A}_{\mathrm{C}}, \\ \forall e \in\left\{e <\alpha_{a}^{t}\right\}, \forall s, c, t . \end{gathered} $

(18c)

其中$\alpha_{a}^{t}$代表充电链路$a$在周期$t$的平均充电时间, 即在时间$\delta$内充电链路$a$的充电桩能给电动汽车供应的能级数量。假设充电链路$a$的充电功率$(\mathrm{kW})$是$P_{a}^{\mathrm{ev}}$, 可以计算出$\alpha_{a}^{t}=\frac{P_{a}^{\mathrm{ev}} \cdot \delta}{\eta \cdot \delta \cdot v_{\mathrm{f}}}=\frac{P_{a}^{\mathrm{ev}}}{\eta \cdot v_{\mathrm{f}}}$。式(18a) 和(18c) 分别约束了更新后能级的上下界。式(18b) 线性地描述了电动汽车能级更新过程。并且, 充电链路$a$上各个能级分流的流出量应当小于链路上电动汽车的充电桩占有量,

$ \begin{gathered} \mathrm{VE}_{a, c}^{s, e}(t)-\mathrm{VE}_{a, c}^{s, e}(t-1) \leqslant x_{a, c}^{s, e}(t), \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall e \in \mathcal{E}_{c}, \forall s, c, t . \end{gathered} $

(19)

充电链路的占有量应当是非负的，

$ \begin{gathered} x_{a, c}^{s, e}(t) \geqslant 0, \hat{x}_{a, c}^{s, e}(t) \geqslant 0, \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall e \in \mathcal{E}_{c}, \forall s, c, t . \end{gathered} $

(20)

累计车流量应当非负且非减：

$ \mathrm{VG}_{a}^{s}(t)-\mathrm{VG}_{a}^{s}(t-1) \geqslant 0, \forall a \in \mathcal{A}, \forall s, t. $

(21a)

$ \begin{aligned} & \mathrm{VE}_{a, c}^{s, e}(t)-\mathrm{VE}_{a, c}^{s, e}(t-1) \geqslant 0, \\ & \forall a \in \mathcal{A}, \forall e \in \mathcal{E}_{c}, \forall s, c, t . \end{aligned} $

(21b)

$ \mathrm{UG}_{a}^{s}(t)-\mathrm{UG}_{a}^{s}(t-1) \geqslant 0, \forall a \in \mathcal{A}, \forall s, t. $

(22a)

$ \begin{aligned} & \mathrm{UE}_{a, c}^{s, e}(t)-\mathrm{UE}_{a, c}^{s, e}(t-1) \geqslant 0, \\ & \forall a \in \mathcal{A}, \forall e \in \mathcal{E}_{c}, \forall s, c, t . \end{aligned} $

(22b)

下面的约束使得初始累计车辆数目为0：

$ \mathrm{UG}_{a}^{s}(0)=\mathrm{VG}_{a}^{s}(0)=0, \forall a \in \mathcal{A}, \forall s . $

(23a)

$ \begin{gathered} \mathrm{UE}_{a, c}^{s, e}(0)=\mathrm{VE}_{a, c}^{s, e}(0)=0, \\ \forall a \in \mathcal{A}, \forall e \in \mathcal{E}_{c}, \forall s, c . \end{gathered} $

(23b)

电气化交通路网模型的目标是最小化所有车辆的总行驶时间。总行驶时间由整个时间范围内所有链路上所有车辆的总停留时间和所有电动汽车的总充电时间计算得到。在正常情况下，电气化交通路网的混合交通流分配模型如下：

$ \left\{\begin{array}{l} \text { 目标函数: 式 (1). } \\ \text { 约束条件: 式 (8)-(23). } \end{array}\right. $

(24)

1.2 双层攻击者-防御者模型

本节提出一种双层攻击者-防御者模型用于寻找出电气化交通路网中脆弱链路，这些链路的损坏将对系统性能造成很大的影响。假设存在一个防御者，即路网的管理者，其目标是最优化道路性能，即在系统最优交通流分配中，最小化系统总行驶时间。假设存在一个攻击者，在实际情况下，可能是自然灾害等，对电气化交通路网的链路进行破坏干涉，使系统性能最大程度地降低。由此，本文可以建立最大化-最小化的双层规划模型来描述上述攻击-防御过程，以此来找出网络中的脆弱链路。其目标函数为：

$ \begin{gathered} \max _{\mu \in \mathcal{M}} \min \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{t \in \mathcal{T}} \sum\limits_{a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}, \mathcal{A}_{\mathrm{S}}\right\}} \delta\left[\mathrm{UG}_{a}^{s}(t)-\mathrm{VG}_{a}^{s}(t)\right]+ \\ \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{t \in \mathcal{T}} \sum\limits_{a \in \mathcal{A} \backslash \mathcal{A}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}} \delta\left[\mathrm{UE}_{a, c}^{s, e}(t)-\mathrm{VE}_{a, c}^{s, e}(t)\right] . \end{gathered} $

(25)

本文仅考虑普通链路的损坏，而不考虑终点链路或是源链路的损坏(由于终点链路与源链路均为虚拟链路)。本文采用0-1变量$\mu_{a}, \forall a \in \mathcal{A}_{\mathrm{G}}$, 来表示网络中链路好坏情况, $\mu_{a}=0$代表链路损坏, 此时道路最大承载数量可认为为0, 因此本文可以将式(11) 替换为

$ \begin{gathered} \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}}\left[\mathrm{UE}_{a, c}^{s, e}(t)-\mathrm{VE}_{a, c}^{s, e}\left(t-\beta_{a}\right)\right]+ \\ \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}}\left[\mathrm{UG}_{a}^{s}(t)-\mathrm{VG}_{a}^{s}\left(t-\beta_{a}\right)\right] \leqslant \mu_{a} \cdot L_{a} k_{\mathrm{jam}}, \\ \forall a \in \mathcal{A} \backslash\left\{\mathcal{A}_{\mathrm{C}}\right\}, \forall t, s . \end{gathered} $

(26)

本文可以得到完整的双层攻击-防御者模型如下：

$ \left\{\begin{array}{l} \text { 目标函数: 式 (25). } \\ \text { 约束条件: 式 (8)-(10), (12)-(23), (26). } \end{array}\right. $

(27)

其中: $\mathcal{M}=\left\{\mu: \sum\limits_{a \in \mathcal{A}_{\mathrm{G}}}\left(1-\mu_{a}\right) \geqslant N_{b}, \mu \in\{0\right.$, $\left.1\}^{\left|\mathcal{A}_{\mathrm{G}}\right|}\right\}, N_{b}$为网络中损坏链路的数目。通过求解问题(27)可以获知发生损坏对系统性能影响最大的几个链路，并且可以得知在这几个链路损坏的情况下，电气化交通路网的系统总行驶时间与交通流分配情况。该双层规划模型可以通过YALMIP等工具建模，用CPLEX等商用求解器进行求解，或者可以直接求解与该双层模型等价的单层混合整数线性规划模型。

1.3 脆弱性评估指标

本文使用2个指标来评估电气化路网的脆弱性。一个是给定时间范围内系统中所有用户总行驶时间，即式(1)。总行驶时间越小，代表在给定时间内和给定交通需求的情况下，所有用户总行驶时间越少，系统的通行时间效率越高。

第2个是车辆随时间的到达率p(t)，

$ p(t)=\frac{\sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{a \in \mathcal{A}_{\mathrm{S}}}\left[\mathrm{UG}_a^s(t)+\sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}} \mathrm{UE}_{a, c}^{s, e}(t)\right]}{\sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{a \in \mathcal{A}_{\mathrm{S}}}\left[\mathrm{DG}_a^s(t)+\sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}} \mathrm{DE}_{a, c}^{s, e}(t)\right]} . $

(28)

在给定时间内和交通需求的情形下，到达率越高，则说明到达终点的用户越多，即系统吞吐性能越高。

总行驶时间和到达率都可以反映系统的性能。前者是一个静态指标，后者是一个动态指标，后者可以反映出系统的更多信息。

2 双层攻击者-防御者模型的求解方法

本文提出的双层攻击者-防御者模型的求解方法首先对内层问题取对偶并与外层问题结合，得到一个二次规划，该二次规划的目标函数为：

$ \begin{gathered} \max \sum\limits_{t \in \mathcal{T}}\left\{\sum\limits_{a \in \mathcal{A}_{\mathrm{R}}} \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}}\left[\mathrm{g} \lambda_{a}^{s}(t) \cdot \mathrm{DG}_{a}^{s}(t)\right]+\right. \\ \sum\limits_{a \in \mathcal{A}_{\mathrm{R}} s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}}\left[\mathrm{e} \lambda_{a, c}^{s, e}(t) \cdot \mathrm{DE}_{a, c}^{s, e}(t)\right]+ \\ \sum\limits_{a \in \mathcal{A} /\left\{\mathcal{A}_{\mathrm{C}}\right\}}\left[\pi_{a}^{1}(t) \cdot f_{a}^{\mathrm{O}}(t)+\pi_{a}^{2}(t) \cdot f_{a}^{\mathrm{I}}(t)-\right. \\ \left.\left.\mu_{a} \cdot \pi_{a}^{3}(t) \cdot L_{a} k_{\mathrm{jam}}\right]+\sum\limits_{a \in \mathcal{A}_{\mathrm{C}}}\left[\pi_{a}^{4} \cdot \mathrm{NC}_{a}(t)\right]\right\} . \end{gathered} $

(29)

其中: $g \lambda_{a}^{s}(t)$为燃油汽车需求满足约束的对偶变量, $\mathrm{e} \lambda_{a, c}^{s, e}(t)$为电动汽车交通需求满足约束的对偶变量。$\pi_{a}^{1}(t) 、\pi_{a}^{2}(t) 、\pi_{a}^{3}(t) 、\pi_{a}^{4}(t)$为道路通行能力约束组的对偶变量, 分别代表着交通流流入约束、交通流流出约束、道路容量约束、充电桩数目约束的对偶变量。对偶二次规划所有变量及其取值如下:

$ \begin{gathered} \mathrm{g} \lambda_{a}^{s}(t), \mathrm{e} \lambda_{a, c}^{s, e}(t), \mathrm{g} \varphi_{i}^{s}(t), \mathrm{e} \varphi_{i, c}^{s, e}(t), \eta_{a, c}^{s, e}(t), \\ v_{a, c}^{s, e}(t), \mathrm{g} \kappa_{a}^{s, 1}(t), \mathrm{g} \kappa_{a}^{s, 2}(t), \mathrm{e} \kappa_{a, c}^{s, e, 1}(t), \\ \mathrm{e} \kappa_{a, c}^{s, e, 2}(t), \forall a, i, s, c, e, t . \end{gathered} $

(30a)

$ \begin{gathered} \pi_{a}^{1}(t), \pi_{a}^{2}(t), \pi_{a}^{4}(t), \mathrm{g} \psi_{a}^{s}(t), \mathrm{e} \psi_{a, c}^{s, e}(t), \\ \zeta_{a, c}^{s, e}(t) \leqslant 0, \forall a, c, e, s, t . \end{gathered} $

(30b)

$ \begin{gathered} \pi_{a}^{3}(t), \mathrm{g} \chi_{a}^{s, 1}(t), \mathrm{g} \chi_{a}^{s, 2}(t), \mathrm{e} \chi_{a, c}^{s, e, 1}(t), \\ \mathrm{e} \chi_{a, c}^{s, e, 2}(t) \geqslant 0, \forall a, c, e, s, t . \end{gathered} $

(30c)

其中: $\mathrm{g} \varphi_{i}^{s}(t) 、\mathrm{e} \varphi_{i, c}^{s, e}(t)$分别为燃油车和电动汽车交通流守恒定律约束的对偶变量, 对应式(15)。$\eta_{a, c}^{s, e}(t) 、v_{a, c}^{s, e}(t)$为充电站占用量更新约束和电动汽车充电过程建模约束的对偶变量, 分别对应式(17)和(18)。$\mathrm{g} \kappa_{a}^{s, 1}(t) 、\mathrm{~g} \kappa_{a}^{s, 2}(t) 、\mathrm{e} \kappa_{a, c}^{s, e, 1}(t)$、$\mathrm{e} \kappa_{a, c}^{s, e, 2}(t)$分别为$\mathrm{UG}_{a}^{s}(t) 、\mathrm{VG}_{a}^{s}(t) 、\mathrm{UE}_{a, c}^{s, e}(t)$、$\mathrm{VE}_{a, c}^{s, e}(t)$的置零约束的对偶变量, 对应式(13b) 和(23) 的将变量置零的约束。由于上述对偶变量对应的约束取等号, 因而对偶变量取值自由。

$\mathrm{g} \psi_{a}^{s}(t) 、\mathrm{e} \psi_{a, c}^{s, e}(t)$为约束(12) 的对偶变量, $\zeta_{a, c}^{s, e}(t)$为约束$(19)$的对偶变量。$\mathrm{g} \chi_{a}^{s, 1}(t) 、\mathrm{~g} \chi_{a}^{s, 2}(t)$、$\mathrm{e} \chi_{a, c}^{s, e, 1}(t) 、\mathrm{e} \chi_{a, c}^{s, e, 2}(t)$分别为$\mathrm{UG}_{a}^{s}(t) 、\mathrm{VG}_{a}^{s}(t)$、$\mathrm{UE}_{a, c}^{s, e}(t) 、\mathrm{VE}_{a, c}^{s, e}(t)$的非减约束的对偶变量, 对应式(21)和(22)。

原问题的每一个变量对应着对偶问题的一条约束。为了简洁, 下面只给出$\mathrm{UE}_{a, c}^{s, e}(t) 、\mathrm{VE}_{a, c}^{s, e}(t)$ $\left(\forall a \in \mathcal{A}_{\mathrm{G}} \cup \mathcal{A}_{C}, \forall s, c, e, t\right)$以及中间变量$x_{a, c}^{s, e}(t)$、$\hat{x}_{a, c}^{s, e}(t)\left(\forall a \in \mathcal{A}_{C}, \forall s, c, e, t\right)$所对应的约束(其余原问题变量对应的约束可以通过下面给出的约束的简单变形得到)。

对于变量$\mathrm{UE}_{a, c}^{s, e}(t), \forall a \in \mathcal{A}_{\mathrm{G}}$, 其约束为:

$ \begin{gathered} \mathrm{e} \varphi_{i, c}^{s, e}(t)-\pi_{a}^{2}(t+1)-\pi_{a}^{3}(t)-\mathrm{e} \psi_{a, c}^{s, e-\rho_{a}}\left(t+\nu_{a}\right)- \\ \mathrm{e} \chi_{a, c}^{s, e, 1}(t+1)+\mathrm{e} \kappa_{a, c}^{s, e, 1}(t) \leqslant 1, \\ \forall a \in \mathcal{A}_{\mathrm{G}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \forall e \in \mathcal{E}_{c}, \forall c \in \mathcal{C}, \\ i=A^{-1}(a), t=0 . \end{gathered} $

(31a)

$ \begin{gathered} \mathrm{e} \varphi_{i, c}^{s, e}(t)+\pi_{a}^{2}(t)-\pi_{a}^{2}(t+1)- \\ \pi_{a}^{3}(t)\left(-\mathrm{e} \psi_{a, c}^{s, e-\rho_{a}}\left(t+\nu_{a}\right)\right)+ \\ \mathrm{e} \chi_{a, c}^{s, e, 1}(t)-\mathrm{e} \chi_{a, c}^{s, e, 1}(t+1) \leqslant 1, \\ \forall a \in \mathcal{A}_{\mathrm{G}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, e <\rho_{a}\left(e \geqslant \rho_{a}\right), \\ \forall c \in \mathcal{C}, i=A^{-1}(a), \forall t \in\left\{1, \cdots, T_{\max }-\nu_{a}\right\}. \end{gathered} $

(31b)

$ \begin{gathered} \mathrm{e} \varphi_{i, c}^{s, e}(t)+\pi_{a}^{2}(t)-\pi_{a}^{2}(t+1)-\pi_{a}^{3}(t)+ \\ \mathrm{e} \chi_{a, c}^{s, e, 1}(t)-\mathrm{e} \chi_{a, c}^{s, e, 1}(t+1) \leqslant 1, \\ \forall a \in \mathcal{A}_{\mathrm{G}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \forall e \in \mathcal{E}_{c}, \\ \forall c \in \mathcal{C}, i=A^{-1}(a), \\ \forall t \in\left\{T_{\max }-\nu_a+1, \cdots, T_{\max }-1\right\} . \end{gathered} $

(31c)

$ \begin{gathered} \mathrm{e} \varphi_{i, c}^{s, e}(t)+\pi_{a}^{2}(t)-\pi_{a}^{3}(t)+\mathrm{e} \chi_{a, c}^{s, e, 1}(t) \leqslant 1, \\ \forall a \in \mathcal{A}_{\mathrm{G}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \forall e \in \mathcal{E}_c, \\ \forall c \in \mathcal{C}, i=A^{-1}(a), t=T_{\max } . \end{gathered} $

(31d)

其中$A^{-1}(a)$为集合$A(i)$的逆映射, 其含义为链路$a$的头节点。当式(31b) 中能级$e$满足$e <\rho_{a}$时, $-\mathrm{e} \psi_{a, c}^{s, e-\rho_{a}} \left(t+\nu_{a}\right)$不存在; 当能级$e$满足$e \geqslant \rho_{a}$时, $-\mathrm{e} \psi_{a, c}^{s, e-\rho_{a}}\left(t+\nu_{a}\right)$存在, 应将括号内的$-\mathrm{e} \psi_{a, c}^{s, e-\rho_{a}}(t+$ $\left.\nu_{a}\right)$添加进约束中。

对于变量$\mathrm{VE}_{a, c}^{s, e}(t), \forall a \in \mathcal{A}_{\mathrm{G}}$, 其约束为:

$ \begin{gathered} -\mathrm{e} \varphi_{i, c}^{s, e}(t)-\pi_{a}^{1}(t+1)+\pi_{a}^{3}\left(t+\beta_{a}\right)+ \\ \mathrm{e} \psi_{a, c}^{s, e}(t)\left(\mathrm{e} \kappa_{a, c}^{s, e, 2}(t)\right)-\mathrm{e}\chi_{a, c}^{s, e, 2}(t+1)+ \\ \mathrm{e} \kappa_{a, c}^{s, e, 2}(t) \leqslant-1, \quad \forall a \in \mathcal{A}_{\mathrm{G}}, \\ \forall s \in \mathcal{N}_{\mathrm{SR}}, e \leqslant E_{c}-\rho_{a}\left(e>E_{c}-\rho_{a}\right), \\ \forall c \in \mathcal{C}, i=B^{-1}(a), t=0 . \end{gathered} $

(32a)

$ \begin{gathered} -\mathrm{e} \varphi_{i, c}^{s, e}(t)+\pi_{a}^{1}(t)-\pi_{a}^{1}(t+1)+\pi_{a}^{3}\left(t+\beta_{a}\right)+ \\ \mathrm{e} \psi_{a, c}^{s, e}(t)\left(\mathrm{e} \kappa_{a, c}^{s, e, 2}(t)\right)+\mathrm{e} \chi_{a, c}^{s, e, 2}(t)- \\ \mathrm{e} \chi_{a, c}^{s, e, 2}(t+1) \leqslant-1, \quad \forall a \in \mathcal{A}_{\mathrm{G}}, \\ \forall s \in N_{\mathrm{SR}}, e \leqslant E_c-\rho_a\left(e>E_c-\rho_a\right), \\ \forall c \in \mathcal{C}, i=B^{-1}(a), \forall t \in\left\{1, \cdots, T_{\max }-\beta_{a}\right\}. \end{gathered} $

(32b)

$ \begin{gathered} -\mathrm{e} \varphi_{i, c}^{s, e}(t)+\pi_{a}^{1}(t)-\pi_{a}^{1}(t+1)+ \\ \mathrm{e} \psi_{a, c}^{s, e}(t)\left(\mathrm{e} \kappa_{a, c}^{s, e, 2}(t)\right)+\mathrm{e} \chi_{a, c}^{s, e, 2}(t)- \\ \mathrm{e} \chi_{a, c}^{s, e, 2}(t+1) \leqslant-1, \quad \forall a \in \mathcal{A}_{\mathrm{G}}, \\ \forall s \in \mathcal{N}_{\mathrm{SR}}, e \leqslant E_{c}-\rho_{a}\left(e>E_{c}-\rho_{a}\right), \forall c \in \mathcal{C}, \\ i=B^{-1}(a), \forall t \in\left\{T_{\max }-\beta_{a}+1, \cdots, T_{\max }-1\right\}. \end{gathered} $

(32c)

$ \begin{gathered} -\mathrm{e} \varphi_{i, c}^{s, e}(t)+\pi_{a}^{1}(t)+\mathrm{e} \psi_{a, c}^{s, e}(t)\left(\mathrm{e} \kappa_{a, c}^{s, e, 2}(t)\right)+ \\ \mathrm{e} \chi_{a, c}^{s, e, 2}(t) \leqslant-1, \forall a \in \mathcal{A}_{\mathrm{G}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \\ e \leqslant E_{c}-\rho_{a}\left(e>E_{c}-\rho_{a}\right), \forall c \in \mathcal{C}, \\ i=B^{-1}(a), t=T_{\max } . \end{gathered} $

(32d)

其中$B^{-1}(a)$为集合$B(i)$的逆映射, 其含义为链路$a$的尾节点。当式(32) 中能级$e$满足$e \leqslant E_{c}-\rho_{a}$时, 说明电动汽车的能级没有超出最大能级, 电动汽车可以流出链路$a$, 对应式(13a), 式(13a) 的对偶变量是$\mathrm{e} \psi_{a, c}^{s, e}(t)$; 当式(32) 中能级$e$满足$e>$ $E_{c}-\rho_{a}$时, 说明电动汽车的能级超出了最大能级, 电动汽车不可流出链路$a$, 对应$\mathrm{VE}_{a, c}^{s, e}(t)$的置零约束式(13b), 此时应将$\mathrm{e} \psi_{a, c}^{s, e}(t)$替换为式(13b) 所对应的对偶变量, 即$\mathrm{e} \psi_{a, c}^{s, e}(t)$后括号内的$\mathrm{e} \kappa_{a, c}^{s, e, 2}(t)$。

对于变量$\mathrm{UE}_{a, c}^{s, e}(t), \forall a \in \mathcal{A}_{\mathrm{C}}$, 其约束为:

$ \begin{gathered} \mathrm{e} \varphi_{i, c}^{s, e}(t)+\eta_{a, c}^{s, e}(t+2)+\pi_{a}^{3}(t)-\mathrm{e} \chi_{a, c}^{s, e, 1}(t+1)+ \\ \mathrm{e}\kappa_{a, c}^{s, e, 1}(t) \leqslant 1, \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \\ \forall e \in \mathcal{E}_{c}, \forall c \in \mathcal{C}, t=0 . \end{gathered} $

(33a)

$ \begin{gathered} \mathrm{e} \varphi_{i, c}^{s, e}(t)+\eta_{a, c}^{s, e}(t+2)-\eta_{a, c}^{s, e}(t+1)+\pi_{a}^{4}(t)+ \\ \mathrm{e} \chi_{a, c}^{s, e, 1}(t)-\mathrm{e} \chi_{a, c}^{s, e, 1}(t+1) \leqslant 1, \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, e <\rho_{a}\left(e \geqslant \rho_{a}\right), \\ \forall c \in \mathcal{C}, \forall t \in\left\{1, \cdots, T_{\max }-2\right\}. \end{gathered} $

(33b)

$ \begin{gathered} \mathrm{e} \varphi_{i, c}^{s, e}(t)-\eta_{a, c}^{s, e}(t+1)+\pi_{a}^{4}(t)+\mathrm{e} \chi_{a, c}^{s, e, 1}(t)- \\ \mathrm{e} \chi_{a, c}^{s, e, 1}(t+1) \leqslant 1, \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \\ \forall e \in \mathcal{E}_{c}, \forall c \in \mathcal{C}, t=T_{\text {max }}-1. \end{gathered} $

(33c)

$ \begin{gathered} \mathrm{e} \varphi_{i, c}^{s, e}(t)+\pi_{a}^{4}(t)+\mathrm{e} \chi_{a, c}^{s, e, 1}(t) \leqslant 1, \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \forall e \in \mathcal{E}_{c}, \\ \forall c \in \mathcal{C}, t=T_{\max } . \end{gathered} $

(33d)

对于变量$\mathrm{VE}_{a, c}^{s, e}(t), \forall a \in \mathcal{A}_{\mathrm{C}}$, 其约束为:

$ \begin{gathered} -\mathrm{e} \varphi_{i, c}^{s, e}(t)-\eta_{a, c}^{s, e}(t+2)-\pi_{a}^{4}(t)- \\ \mathrm{e} \chi_{a, c}^{s, e, 2}(t+1)-v_{a, c}^{s, e}(t)+ \\ \mathrm{e} \kappa_{a, c}^{s, e, 2}(t+1) \leqslant-1, \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \\ \forall e \in \mathcal{E}_{c}, \forall c \in \mathcal{C}, i=B^{-1}(a), t=0 . \end{gathered} $

(34a)

$ \begin{gathered} -\mathrm{e} \varphi_{i, c}^{s, e}(t)-\eta_{a, c}^{s, e}(t+2)+\eta_{a, c}^{s, e}(t+1)-\pi_{a}^{4}(t)+ \\ \mathrm{e} \chi_{a, c}^{s, e, 2}(t)-\mathrm{e} \chi_{a, c}^{s, e, 2}(t+1)+v_{a, c}^{s, e}(t)- \\ v_{a, c}^{s, e}(t+1) \leqslant-1, \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \forall e \in \mathcal{E}_{c}, \forall c \in \mathcal{C}, \\ i=B^{-1}(a), \forall t \in\left\{1, \cdots, T_{\max }-2\right\} . \end{gathered} $

(34b)

$ \begin{gathered} -\mathrm{e} \varphi_{i, c}^{s, e}(t)+\eta_{a, c}^{s, e}(t+1)-\pi_{a}^{4}(t)+\mathrm{e} \chi_{a, c}^{s, e, 2}(t)- \\ \mathrm{e} \chi_{a, c}^{s, e, 2}(t+1)+v_{a, c}^{s, e}(t)-v_{a, c}^{s, e}(t+1) \leqslant-1, \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \forall e \in \mathcal{E}_{c}, \forall c \in \mathcal{C}, \\ i=B^{-1}(a), t=T_{\max }-1 \text {. } \end{gathered} $

(34c)

$ \begin{gathered} -\mathrm{e} \varphi_{i, c}^{s, e}(t)-\pi_{a}^{4}(t)+\mathrm{e} \chi_{a, c}^{s, e, 2}(t)+v_{a, c}^{s, e}(t) \leqslant-1, \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \forall e \in \mathcal{E}_{c}, \\ \forall c \in \mathcal{C}, i=B^{-1}(a), t=T_{\text {max }} . \end{gathered} $

(34d)

对于变量$x_{a, c}^{s, e}(t)$, 其约束为:

$ \begin{gathered} \xi_{a, c}^{s, e}(t)-\eta_{a, c}^{s, e}(t+1)-v_{a, c}^{s, e}(t) \leqslant 0, \\ \forall a \in \mathcal{A}_{\mathrm{C}}, \forall s \in \mathcal{N}_{\mathrm{SR}}, \\ \forall e \in \mathcal{E}_{c}, \forall c \in \mathcal{C}, i=B^{-1}(a), \\ \forall t \in\left\{1, \cdots, T_{\max }-1\right\} . \end{gathered} $

(35a)

$ \begin{gathered} \xi_{a, c}^{s, e}(t)-v_{a, c}^{s, e}(t) \leqslant 0, \forall a \in \mathcal{A}_{\mathrm{C}}, \\ \forall s \in \mathcal{N}_{\mathrm{SR}}, \forall e \in \mathcal{E}_{c}, \\ \forall c \in \mathcal{C}, i=B^{-1}(a), t=T_{\max } . \end{gathered} $

(35b)

对于变量$\hat{x}_{a, c}^{s, e}(t)$, 其约束为:

$ \begin{gathered} -\xi_{a, c}^{s, E_{c}}(t)+\eta_{a, c}^{s, e}(t) \leqslant 0, \forall a \in \mathcal{A}_{\mathrm{C}}, \\ e \geqslant E_{c}-\alpha_{a}^{t}, \forall s, c, t . \end{gathered} $

(36a)

$ \begin{gathered} -\xi_{a, c}^{s, e+\alpha_{a}^{t}}(t)+\eta_{a, c}^{s, e}(t) \leqslant 0, \forall a \in \mathcal{A}_{\mathrm{C}}, \\ e <E_{c}-\alpha_{a}^{t}, \forall s, c, t . \end{gathered} $

(36b)

式(29)-(36) 即本文中双层模型的等价混合整数二次规划模型, 目标函数的二次项可被进一步线性化为一次项, 使得混合整数二次规划模型转化为混合整数线性规划模型。引入一变量$\gamma_{a}(t)=\mu_{a}\cdot$$\pi_{a}^{3}(t)$, 其中$\pi_{a}^{3}(t)$为大于0的决策变量, 其含义为道路$a$的容量变化1单位对目标函数的影响, 因而可以认为$\pi_{a}^{3}(t) \in[0, M]$, 其中$M$为预设的足够大的正整数, 一般情况下可以将$M$设为1000。对于新引入的变量$\gamma_{a}(t)$, 基于大$M$法, 本文有着以下的线性化约束组:

$ \gamma_{a}(t) \leqslant \pi_{a}^{3}(t), \quad \forall a \in \mathcal{A}_{\mathrm{G}}, \forall t . $

(37a)

$ \gamma_{a}(t) \leqslant M \cdot \mu_{a}, \quad \forall a \in \mathcal{A}_{\mathrm{G}}, \forall t . $

(37b)

$ \gamma_{a}(t) \geqslant \pi_{a}^{3}(t)+M \cdot\left(\mu_{a}-1\right), \forall a \in A_{\mathrm{G}}, \forall t . $

(37c)

当$\mu_{a}=0$时, 式(37b) 使得$\gamma_{a}(t) \leqslant 0$, 式(37c) 使得$\gamma_{a}(t) \geqslant \pi_{a}^{3}(t)-M \geqslant 0$, 因而可得$\gamma_{a}(t)=0$。当$\mu_{a}=1$时, 式(37a) 使得$\gamma_{a}(t) \leqslant \pi_{a}^{3}(t)$, 式(37c) 使得$\gamma_{a}(t) \geqslant \pi_{a}^{3}(t)$, 因而$\gamma_{a}(t)=\pi_{a}^{3}(t)$。原问题被线性化为

$ \begin{aligned} & \max _{\mu \in \mathcal{M}} \sum\limits_{t \in \mathcal{T}} \sum\limits_{a \in \mathcal{A}_{\mathrm{R}}}\left\{\sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}}\left[\mathrm{g} \lambda_{a}^{s}(t) \cdot \mathrm{DG}_{a}^{s}(t)\right]+\right. \\ & \sum\limits_{a \in \mathcal{A}_{\mathrm{R}}} \sum\limits_{s \in \mathcal{N}_{\mathrm{SR}}} \sum\limits_{c \in \mathcal{C}} \sum\limits_{e \in \mathcal{E}_{c}}\left[\mathrm{e} \lambda_{a, c}^{s, e}(t) \cdot \mathrm{DE}_{a, c}^{s, e}(t)\right]+ \\ & \sum\limits_{a \in \mathcal{A}_{\mathrm{C}}}\left[\pi_{a}^{1}(t) \cdot f_{a}^{\mathrm{O}}(t)+\pi_{a}^{2}(t) \cdot f_{a}^{\mathrm{I}}(t)-\right. \\ & \left.\left.\gamma_{a}(t) \cdot L_{a} k_{\mathrm{jam}}\right]+\sum\limits_{a \in \mathcal{A}_{\mathrm{C}}}\left[\pi_{a}^{4} \cdot \mathrm{NC}_{a}(t)\right]\right\} . \end{aligned} $

(38)

$ \begin{array}{l} \text { s.t. } \\ & \text { 式(31)-(37). } \end{array} $

该单层混合整数线性规划问题与原双层问题等效，并且其求解速率更高。

3 实验

本文改进了文[29]中的美国北卡罗来纳州的部分高速路网的例子来验证提出的方法。图 2展示了这个区域的高速路网^[30]被抽象化后的拓扑结构。该电气化路网的节点ID对应城镇的名字和人口数据^[29]列在表 1中，高速路网参数和模型所用的参数分别列在表 2和3中。考虑到这些城镇之间的地理距离和人口，重力模型被用来生成每天的交通需求。重力模型的一般形式^[31]被写作：$f_{o d}=P_{o}^{\alpha} P_{d}^{\beta} /$ $D_{o d}^{\gamma}$。其中: $P_{o}$和$P_{d}$分别是起点$o$和终点$d$的人口规模, $D_{o d}$是它们之间的最短路径, $\alpha 、\beta$和$\gamma$是拟合参数。本文设置$\alpha=\beta=0.92, \gamma=1$。为了考虑最坏的情况，本文采用17：00到20：00时的交通需求，这个时间段的交通需求占到了全天交通需求的大约15.3%^[32]。不同方向的交通流量通常具有差异性，获得具体的数据是非常困难的^[32]。因而，为了简化，本文仅随机考虑每个od对的一个方向，而另一个方向的交通流量被忽略。获得的交通出行需求如表 4所示。

所有的数值实验在处理器为Intel Core i7-11700、2.50 GHz，内存为32 GB的个人计算机上运行。所有的问题都使用商业软件IBM ILOG CPLEX(版本12.10.0)求解。

3.1 不同电动汽车滲透率对比

本节假设攻击者最多只能攻击路网中的一条链路，以此来研究在不同电动汽车滲透率情况下的电气化路网的脆弱性。

图 3展示了在不同电动汽车滲透率情况下车辆到达率随着时间的变化。在时间等于0的时刻，系统中没有车辆到达终点。随着时间的流逝，到达终点的车辆数量逐渐增加。表 5展示了在不同电动汽车滲透率下电气化路网中最脆弱的链路。如图 3和表 5所示，随着电动汽车滲透率的升高，系统的整体性能(总行驶时间和车辆到达率)下降。当电动汽车滲透率从0%上升到100%，系统总行驶时间从189 492 s上升到473 734 s，车辆到达率从0.81急剧下降到0.13。造成这一现象的原因主要有2个：1) 相较于传统的非电动汽车，电动汽车需要额外的充电时间；2) 系统中的充电桩数量是有限的，因此随着电动汽车滲透率的上升，系统中的充电站出现拥挤导致较长的额外排队时间。充电桩的数量限制了该系统的整体性能，使得系统性能随着电动汽车滲透率的上升出现了急剧下降。当电动汽车滲透率大于等于25%时，系统中最脆弱的链路稳定在105号。当路网中没有电动汽车时，最脆弱的链路为110号。考虑电动汽车与不考虑电动汽车时，系统中最脆弱的链路完全不同。由此可见，在评估电气化路网中最脆弱的链路时，仅考虑非电动汽车流是不够全面的，把电动汽车流考虑在内是有必要的。

3.2 不同攻击资源等级对比

为了不丢失一般性，本节假设电动汽车滲透率为50%，以此研究攻击者的资源等级破坏0到4条链路的情况。如图 4和表 6所示，随着攻击资源等级从0上升到4，系统总行驶时间从285 196 s增加到373 924 s，车辆到达率从0.54下降到0.35。表 6展示了低攻击资源等级的最脆弱链路的集合并不一定是高攻击资源等级的最脆弱链路集合的子集。比如，当攻击资源等级为2时，最脆弱链路的集合是3号和105号；而当攻击资源等级为3时，最脆弱链路的集合为9号、16号和117号。如图 4所示，当攻击资源等级小于等于2时，系统中车辆到达率的下降并不明显；当攻击资源等级从2变为3时，车辆到达率急剧下降；当攻击资源等级从3变为4时，车辆到达率的下降程度又变得不很明显。这样的现象印证了复杂系统中存在的临界点和相变现象。在复杂系统中，初始的攻击可能会造成系统以级联失效的方式崩溃，这种崩溃方式往往表现为相变的形式。在本例中，当攻击资源达到到某一等级的临界点时，系统中可能出现了一定范围内的级联失效，导致系统性能的急剧下降。也就是说，当路网中某些关键道路的集合失去通行能力时，会导致大范围内车辆的出行受到阻碍，很大程度地影响了整个路网的运行性能。对于改善电气化路网的脆弱性问题，识别出系统中这些临界点参数和关键道路集合尤为重要。

4 结论

本文提出一种双层攻击者-防御者混合整数规划模型来研究电气化路网的脆弱性。外层模型通过攻击系统中的道路，使其失去通行能力，从而最大程度降低系统的性能，即最大化系统总行驶时间。内层模型作为防御者，通过动态地最优化分配包含电动汽车和非电动汽车的交通流量，使系统总行驶时间最小化。此外，本文还给出了所提出的模型的详细求解方法和理论分析。通过对内层问题取对偶并与外层问题结合，得到一个混合整数二次规划问题，并进一步通过大M法，将该问题转化成混合整数线性规划问题。最后，将提出的模型应用于美国北卡罗莱纳州的部分高速路网，实验结果表明：1) 考虑电动汽车与不考虑电动汽车时，得到的关键道路是完全不同的。因此，在分析路网的脆弱性时，有必要把电动汽车考虑在内。2) 在不同攻击资源等级情况下，关键道路的集合是不同的。3) 实验结果验证了攻击资源等级参数存在临界点，达到该临界点，系统性能出现相变现象，即系统性能出现大幅下降。以上研究结果，可以为改善电气化路网的脆弱性提供理论支撑。

本文的工作还可以从以下几个方向进行拓展：1) 虽然目前该工作已经把原始模型的双层混合整数规划问题转换为混合整数线性规划问题，且相较于原模型，其求解效率已有较大提高，但考虑到内层问题是一个复杂的动态混合交通流分配问题，其求解效率还可进一步提高。在接下来的工作中，可以考虑使用Benders解耦算法或列生成算法，进一步优化求解速度。2) 电气化路网与电网是紧密相连的。电气化路网的扰动可以通过电动汽车的充电行为传播到电网，反之亦然。如果将电网和路网视为一个整体，全面地分析其整体的脆弱性，识别整个耦合系统中最脆弱的组件(道路、充电站、电线等)，将是一个更具挑战性但也更有指导意义的工作。3) 本文建立的是一个确定性的模型，在未来的工作中，可以将交通需求的不确定性考虑在内，使得模型能够考虑更加复杂的现实情况。4) 改善电气化交通路网的脆弱性是十分重要且具有实际意义的工作。本研究是改善交通路网脆弱性的一项必不可少的前置工作，具体的改善交通路网脆弱性的优化方法，还需要结合实际的资源约束(如时间、人力和花费等)和物理系统约束(如道路的宽度、位置和长度)等进一步研究，这也是未来非常值得探索的研究方向。