数据-模型协同驱动的Robotaxi车队充电设施选址定容方法

摘要

自动驾驶电动出租车Robotaxi车队的规模化商业运营需要完善的充电基础设施作为前提，但目前充电设施资源仍然存在容量不足、利用率低、布局不合理等问题。针对Robotaxi车队充电站选址定容问题，首先，以满足乘客出行需求为约束条件，提出基于出行网络行程衔接和Hopcroft-Karp算法的最小车队规模计算方法；在此基础上，通过蒙特卡洛仿真量化获取车队充电需求的时空分布；然后，考虑充电站建设投资、运营、维护的折旧成本以及车队空载行驶、排队充电和订单损失的机会成本，以综合成本最小化为目标，构建充电站选址定容优化模型，并针对模型求解提出一种基于遗传算子和自适应惯性权重的改进粒子群算法；最后，基于中国成都市的用户出行订单和地理数据验证了所提模型及算法的有效性。

关键词

自动驾驶出租车; 充电设施选址; 粒子群优化

新能源汽车的发展承载着我国汽车产业转型升级、降低化石资源依赖、改善生态环境等多重任务。据公安部统计，截至2022年6月底，全国汽车保有量达3.10亿辆^［1］，比2012年增加近2倍。然而，汽车尤其是电动汽车（EV）保有量的快速增长也显露出一系列挑战，特别是充电基础设施建设的滞后。对此，国务院办公厅印发了《新能源汽车产业发展规划（2021—2035年）》，指出要加快推进充换电基础设施的科学布局^［2］。随着电动化、智能化、网联化、共享化等创新技术在城市交通系统的应用于发展，自动驾驶电动出租车Robotaxi车队应运而生。相较于传统出租车，Robotaxi可以在共享模式下以较小的车队规模满足更多的乘客出行需求，显著提高车辆利用率，大幅削减劳动力成本，因而被认为是一种灵活、智能、节能、经济的全新出行方式^［3］。Robotaxi的大规模运营需要完善的充电设施作为保障，然而，目前充电设施的开发大多集中在传统人类驾驶车辆上，少有考虑无人驾驶的运营场景。

基于上述背景，本研究针对Robotaxi车队充电设施选址定容问题开展研究。首先，以满足乘客出行需求为约束条件，提出基于出行行程网络衔接和Hopcroft-Karp算法的最小车队规模计算方法；在此基础上，通过蒙特卡洛仿真获取车队充电需求的时空分布；然后，考虑充电站建设投资、运营维护和网络损耗成本及车队空载行驶、排队充电和订单损失机会成本，以综合成本最小化为目标构建充电站选址定容优化模型，并针对模型求解提出一种基于遗传算子和自适应惯性权重的改进粒子群优化（IPSO）算法；最后，基于中国成都市用户出行订单和地理数据对所提模型和算法的有效性进行了全面验证。

1 文献综述

目前，国内外针对充电站布局规划已有大量研究。考虑到本研究主要涉及车队充电需求预测、选址定容优化模型构建和求解算法设计等层面，因此，针对以上问题对现有研究进展综述如下。

在充电需求预测方面，Zhang等^［4］提出了一种基于准动态交通流的EV充电需求时空分布预测方法，同时考虑了外部环境和交通拥堵对EV出行的影响，较好地平衡了仿真精度和计算效率之间的问题。针对用户随机性和间歇性的出行模式，Zhou等^［5］开发了EV充放电行为预测模型，并在英国电车行为数据集上对模型进行了验证。Guo等^［6］采用成都二环区域的相关数据，基于兴趣点分布间接表达充电需求的静态关联关系，并由用户出行行为衡量动态的实际充电需求，进而采用车辆历史行驶轨迹推断可能的充电需求。可以发现，上述面向电动出租车的充电设施布局研究，均未考虑Robotaxi运营场景，为了获取研究区域内的充电需求，以往研究一般基于车辆充电日志或电网负荷数据。然而，这些方法并不适用于Robotaxi，因为目前难以获取大量真实的Robotaxi运营数据。因此，本研究从用户出行需求角度出发，考虑Robotaxi车队的运营和充电特征，建立蒙特卡洛仿真流程量化获取车队充电需求的时空分布。

在选址定容模型构建方面，Zhao等^［7］以最大化充电站总利润为目标，采用改进的Huff重力模型描述司机的概率充电行为，并据此建立双层优化模型确定新充电站的位置和容量。刘等^［8］利用Huff模型分析用户对不同充电站的选择概率，并基于用户选择概率确定充电站服务范围。任等^［9］以充电成本最小化为优化目标，提出了一种考虑电动出租车运行轨迹、车辆进入停车场时序性和停车场现有负荷水平的电动出租车充电桩规划方法。苏等^［10］提出了一种基于轨迹挖掘的电动出租车充电站规划方法，综合考虑充电站建设运行成本、车辆到站时间成本及等待时间成本。李等^［11］综合考虑EV流量和初始剩余电量的随机特性，引入扩张网络理论构建了多路径条件下潜在的换电站选址数学模型。盛等^［12］从宏观角度对车辆出行-充电行为进行了分析，提出了考虑充电站建设运营成本与用户充电时间成本的充电设施扩展规划架构。以上文献梳理表明，一方面，针对当前充电站利用率低、盈利难、用户充电体验差等问题，有必要综合考虑站点和车队双方利益，统筹充电站建设投资、运营维护和网络损耗成本以及车队空载行驶、排队充电和订单损失机会成本等多优化目标。另一方面，充电设施规划时不仅需要考虑交通特性，还需考虑充电站的负荷约束^［13］。

在模型求解方面，Pahlavanhoseini等^［14］将充电站规划建模为一个混合整数非线性规划问题，并提出了一种采用了混合交叉算子和动态非均匀突变算子的协同进化遗传算法。Pan等人^［15］通过模拟EV驾驶员的充电选择行为建立了一个EV公共充电站选址模型，利用遗传算法进行求解。Zhang等^［16］将用户的偏好和等待时间纳入充电站选址规划模型中，并提出了一种结合了嵌入式水平确定算法和部分枚举算法的混合非主导排序遗传算法II来求解该模型。Gan等^［17］将快速充电站部署问题表述为一个非线性整数问题，并采用一种基于遗传算法的启发式算法进行求解。齐琳等^［18］构建了以充电站服务成本和用户需求成本之和最小为目标的优化模型，采用改进的K中心点云模型混合自适应PSO算法进行求解。Reddy等^［19］提出了一种基于PSO算法寻找EV充电站最优位置和最优容量的方法。综上所述，单独采用某一启发式算法容易存在收敛速度慢或陷入局部最优等缺点，难以满足实际求解需求。因此，相关研究或针对某一特定算法进行针对性改进，或将不同的算法相结合以兼顾收敛速度和求解质量的优化需求。本研究针对PSO算法容易陷入局部最优的缺点，在算法迭代过程中引入选择、交叉、变异等遗传算子，并设计了自适应惯性权重取值机制，以增强种群多样性。案例分析结果表明，上述定制化改进能够有效提高算法全局寻优能力。

2 模型构建

为构建充电站选址布局优化模型，需首先确定车队充电需求的时空分布，其中的关键因素是车队规模。因此，首先通过出行行程网络衔接计算满足乘客出行需求的Robotaxi最小车队规模，进而基于蒙特卡洛仿真获取车队充电需求的时空分布，最终构建以综合成本最小化为目标的充电站选址定容优化模型。模型构建涉及的参数与辅助变量如表1所示。

表1  模型参数与辅助变量

Tab.1  Model parameters and auxiliary variables

符号	参数和辅助变量
$ch$	充电桩充电功率/（kWh·min-1）
$CE$	电力传输效率/%
$d_{ij}$	从位置$i$到位置$j$所需的车辆行驶距离/km
$EF$	二氧化碳排放系数/（t·kWh-1）
$M$	研究区域内充电需求点总数
$n_{\mathrm{i}}$	位置$i$的充电需求总数
$nu{m}_{min}$	单个充电站内配置的充电桩最少数目
$nu{m}_{max}$	单个充电站内配置的充电桩最大数目
$N_{ij}$	从充电需求点$i$前往充电站$j$进行充电的充电需求总数
$p_{invest}$	车辆单位时间投入成本/（元·min-1）
$p_{income}$	车辆单位时间产生收益/（元·min-1）
$p_q$	车辆充电成本/（元·kWh-1）
$Plca{e}_{time}$	车辆当前所在位置
$q$	车辆单位里程能耗/（kWh·km-1）
$SO{C}_0$	车辆初始电池荷电状态/%
$SO{C}_{threshold}$	判断车辆是否产生充电需求的电池电量阈值/%
$SO{C}_{time}$	车辆当前电池荷电状态/%
$Su{m}_{charge}$	研究区域内全天充电需求量
$T$	充电站全天充电服务总时数（$T=24$ h）
$TE$	车辆充电效率/$\mathrm{\%}$
$Tim{e}_{max}$	车辆全天运营时长（$Tim{e}_{max}=1440$ min）
$W_q$	车辆充电平均等待时间/min
$v$	车辆平均行驶速度/（$\mathrm{k}\mathrm{m}·{\mathrm{h}}^{-1})$
$\eta$	碳排放成本系数（$\eta =0.0$1）
$\epsilon$	车辆所在区域接单概率（$\epsilon \in \left[\mathrm{0,1}\right]$）
符号	决策变量
$e_{ab}$	行程$T_a$和行程$T_b$之间是否存在有向连接（二进制变量）
$la{t}_i$	充电站$i$的纬度坐标（基于GCJ-02坐标系）
$ln{g}_i$	充电站$i$的经度坐标（基于GCJ-02坐标系）
$nu{m}_i$	充电站$i$内配置的充电桩总数（正整数变量）
$S$	研究区域内建设的充电站总数（正整数变量）

2.1　最小车队规模核算

计算满足乘客出行需求所需最小车队规模的内在逻辑可分为以下4个步骤：乘客出行状态定义、出行网络行程衔接、有向无环图构建和最小车队规模核算。

定义乘客$i$的单次行程为$T_i$（$t_i^p$，$t_i^d$，$l_i^p$，$l_i^d$），$\mathrm{其中}t$为时间变量，$l$为位置变量，$t_i^p$，$t_i^d$，$l_i^p$，$l_i^d$分别代表乘客$i$的上、下车时间和上、下车位置。定义乘客出行集合为$\Gamma (T_i\in \Gamma )$，对于一个给定的出行集合$\Gamma$，其表示区域某个时间间隔内的全部交通出行需求。为减小车队规模，需要考虑每辆车连续服务多个行程的情形。假设一辆车连续服务两个行程$T_a$和$T_b$，连接行程$T_a$和$T_b$的时间为连接时间$t_{ab}(\mathrm{也即}\mathrm{由位}\mathrm{置}{l}_a^d\mathrm{到位}\mathrm{置}{l}_b^p\mathrm{的车}\mathrm{辆行}\mathrm{程时}\mathrm{间}),t_{ab}=t_b^p-t_a^d$。显然$t_{ab}$过大会导致车辆空载时间太长、车辆利用率过低。因此将最小车队规模问题定义为：找到一个最小的车队规模满足出行集合$\Gamma$，并且满足以下约束：① 在$t_i^p$时刻之前，有一辆可用车辆可以到达位置$l_i^p$为乘客$i$提供服务；② 连接时间$t_{ab}$不超过$\delta$。

由于时间的无环性质，可基于乘客出行网络定义有向无环图$\mathrm{D}\mathrm{A}\mathrm{G}$=（$N$， $E$），其中，$\mathrm{N}$为节点集合，$\mathrm{E}$为有向边集合。节点$n_a,n_b$（$n_a,n_b$∈$N$）分别对应出行需求$T_a,T_b$，当且仅当$t_a^d+t_{ab}\le t_b^p$（即上述约束1）及$t_b^p-t_a^d\le \delta$（即上述约束2）时，两次行程可由一辆车提供，此时$\mathrm{节点}{n}_a$和$n_b$之间存在有向连接，即存在有向边$e_{ab}$∈$E$，行程衔接拓扑关系构建流程如算法1所示（见表2）。

表2  算法1 行程衔接拓扑关系构建

Tab.2  Construction of trip articulation topological relations

步骤	构建说明
输入	乘客出行行程信息$\Gamma$
参数	最大连接时间$\delta$
输出	乘客出行行程拓扑关系$\mathrm{D}\mathrm{A}\mathrm{G}$
1：	For $T_a$ in$\Gamma$： #遍历行程集合$\Gamma$中的每个行程
2：	$n_a$=$T_a$#构建有向无环图节点$n_a$
3：	$\mathrm{D}\mathrm{A}\mathrm{G}.N.\mathrm{a}\mathrm{d}\mathrm{d}(n_a)$#将节点$n_a$添加至有向无环图节点集合$\mathrm{N}$中
4：	End #结束For循环
5：	For $n_a$ in $\mathrm{D}\mathrm{A}\mathrm{G}.N$：#以$N$的每个节点$n_a$为起点
6：	For $n_b$ in $\mathrm{D}\mathrm{A}\mathrm{G}.N$ and b！=a #以$N$中的每个节点$n_b$为终点，且$n_a$！=$n_b$
7：	If $n_b.t_b^p-n_a.t_a^d\le \delta$： #判断节点$n_a$的完成时间和节点$n_b$的开始时间是否超过最大连接时间$\mathrm{\delta }$
8：	$e_{ab}=1$#赋值$e_{ab}\mathrm{为}1$，即$\mathrm{节点}{n}_a\mathrm{至节}\mathrm{点}{n}_b$存在有向边
9：	$\mathrm{D}\mathrm{A}\mathrm{G}.E.\mathrm{a}\mathrm{d}\mathrm{d}({\mathrm{e}}_{\mathrm{a}\mathrm{b}})$#将$e_{ab}$添加至有向无环图边集合$\mathrm{E}$中
10：	Else：
11：	$e_{ab}=0$#赋值$e_{ab}\mathrm{为}0$，即$\mathrm{节点}{n}_a\mathrm{到节}\mathrm{点}{n}_b$不存在有向边
12：	$\mathrm{D}\mathrm{A}\mathrm{G}.E.\mathrm{a}\mathrm{d}\mathrm{d}(e_{ab})$#将$e_{ab}$添加至有向无环图边集合$E$中
13：	End if #结束If循环
14：	End #结束For循环
15：	End #结束For循环
16：	Return #返回$\mathrm{D}\mathrm{A}\mathrm{G}$

将DAG中的节点$N$拆分为$N_x$和$N_y$，分别保留节点$N$有向边的指出方向和指入方向。对于节点$A\to B$的有向边，可以分别将$A$和$B$拆分为$A_x$，$A_y$，$B_x$，$B_y$，将该有向边表示为$A_x\to B_y$。对于一个不对外指出有向边的节点$C$，也需要对它建立一个有向连接以确保输入的完整性，把它拆分为$C_x$和$C_y$，并将有向连接表示为$C_x\to C_y$。通过上述方法，可以将DAG中的所有节点建立连接，也就将DAG最小路径覆盖问题转化为二分图最大匹配问题，该问题可以通过Hopcroft-Karp算法进行快速有效的求解^［20］，DAG中最小不相交路径条数就等于节点总数减去二分图最大匹配数。

2.2　车队充电需求仿真

本部分首先定义出行特征矩阵描述乘客出行事件，进而构建基于蒙特卡洛模拟的乘客出行和车辆运行仿真流程，在此基础上，通过设定Robotaxi充电需求产生的判定条件，最终得到车队充电需求的时空分布。

单次乘客出行事件可以由起止时间、起止地点和出行距离等参数进行完整描述，这些参数可以通过起始概率矩阵$\mathit{S}$、转移概率矩阵${\mathit{P}}^t$和OD距离矩阵${\mathit{D}}^t$来捕捉。定义乘客出行概率矩阵${\mathit{P}}^t$如式（1）所示，矩阵元素$p_{ij}^t$表示时刻$t$以位置$i$为起点的出行需求以位置$j\mathrm{为终}\mathrm{点}$的概率，显然矩阵元素$p_{ij}^t$满足约束（2）。

$$\begin{array}{c}{\mathit{P}}^t=\left[\begin{array}{cccc}& & & \\ p_{\mathrm{11}}^t& p_{\mathrm{12}}^t& \cdots & p_{\mathrm{1}n}^t\\ p_{\mathrm{21}}^t& p_{\mathrm{22}}^t& \cdots & p_{\mathrm{2}n}^t\\ ⋮& ⋮& p_{ij}^t& ⋮\\ p_{n\mathrm{1}}^t& p_{n\mathrm{2}}^t& \cdots & p_{nn}^t\end{array}\right]i,j\in n\end{array}$$

（1）

$${\displaystyle \sum _i^n}{p}_{ij}^t={\displaystyle \sum _j^n}{p}_{ij}^t=\mathrm{1}$$

（2）

定义起点概率矩阵S描述车辆一天首次接单的时间和位置，如式（3）所示。在本研究中，S被设定为24×$n$的二阶矩阵，矩阵元素$s_{ij}$表示车辆在时刻$i$位置$j$接到全天第一个订单的概率。

$$\begin{array}{c}\mathit{S}=\left[\begin{array}{cccc}{s}_{\mathrm{1,1}}& s_{\mathrm{1,2}}& \cdots & s_{\mathrm{1},n}\\ s_{\mathrm{2,1}}& s_{\mathrm{2,2}}& \cdots & s_{\mathrm{2},n}\\ ⋮& ⋮& s_{i,j}& ⋮\\ s_{\mathrm{24,1}}& s_{\mathrm{24,2}}& \cdots & s_{\mathrm{24},n}\end{array}\right]i\in \left[\mathrm{1,24}\right],j\in n\end{array}$$

（3）

OD距离矩阵${\mathit{D}}^t$为n阶方阵，矩阵元素$d_{ij}$表示由起点位置$i$到达目的地位置$j$所需的车辆行驶里程，如式（4）所示。在本研究中，${\mathit{D}}^t$的元素取值通过调用高德开放平台API实现（https：//lbs.amap.com/api/webservice/summary），以考虑交通状况对车辆路线选择的影响。

$$\begin{array}{c}{\mathit{D}}^t=\left[\begin{array}{cccc}& & & \\ d_{\mathrm{11}}^t& d_{\mathrm{12}}^t& \cdots & d_{\mathrm{1}n}^t\\ d_{\mathrm{21}}^t& d_{\mathrm{22}}^t& \cdots & d_{\mathrm{2}n}^t\\ ⋮& ⋮& d_{ij}& ⋮\\ d_{n\mathrm{1}}^t& d_{n\mathrm{2}}^t& \cdots & d_{nn}^t\end{array}\right]i,j\in n\end{array}$$

（4）

蒙特卡洛仿真通过生成许多与其发生概率相匹配的事件来模拟现实情况。基于所构建的乘客出行特征矩阵，Robotaxi车队充电需求仿真流程如图1所示。该过程包括乘客出行和Robotaxi运行两个事件。在模拟过程中，当出现订单需求时，乘客出行事件产生。Robotaxi在接单时需提前判断电池状态，如果当前电池剩余电量高于阈值$SO{C}_{threshold}$，且能够满足该订单所需的车辆行驶里程，则由该车完成该次订单。否则，该出行订单无法被该辆车完成，同时在当前时间和当前位置产生车辆充电需求并记录。按照时间递增顺序推进这两个事件，每次迭代结束可以得到一辆Robotaxi的全天充电需求分布。重复仿真模拟多辆车在同一天的运行情况，得到Robotaxi车队充电需求的时空分布。

图1  Robotaxi车队充电需求仿真流程

Fig.1  Robotaxi Fleet Charging Demand Simulation Process

2.3　选址定容模型构建

本节构建以综合成本最小化为目标的充电设施选址定容优化模型，其中，综合成本包含充电站建设投资成本、运营维护成本和网络损耗成本，以及Robotaxi车队充电过程产生的空载行驶成本、订单损失机会成本和排队充电成本，具体如下。

2.3.1　建设投资成本

充电站建设投资成本指充电站建设的土地成本和充电桩等固定资产投资成本。根据国家标准《电动汽车充电站设计规范GB 50966-2014》，本研究设定Robotaxi充电站平面布局如图2所示。

图2  Robotaxi充电站平面布局

Fig.2  Robotaxi Charging Station Layout

设充电站$i$配置的充电桩总数为$nu{m}_i$，综合考虑充电站消防设施、变压装置等占地面积，则充电站$i$占地面积$S_i$可通过式（5）得到。

$$S_i=\mathrm{133}+\mathrm{23}nu{m}_i$$

（5）

设土地价格为$C{(i)}_{land}$，则充电站$i$的土地成本$C{(i)}_{allland}$可通过式（6）得到。

$$C{\left(i\right)}_{allland}=C{\left(i\right)}_{land}\times S_i$$

（6）

设单个充电桩价格为$C_{pile}$，则充电站$i$的充电桩采购成本$C{(i)}_{allpile}$可通过式（7）得到。

$$C{\left(i\right)}_{allpile}=C_{pile}\times nu{m}_i$$

（7）

对于充电桩安装、供配电系统、监控及通信系统、消防灭火设施等其它费用，假设其与充电站配置的充电桩数量呈线性相关关系，如公式（8）所示。式中折算系数$\alpha$设为$\mathrm{0.03}\times {\mathrm{10}}^{\mathrm{6}}$。

$$C{\left(i\right)}_{other}=\alpha \times nu{m_i}^{\mathrm{2}}$$

（8）

最后，将充电站建设投资成本分摊至最大运营年限year内，同时考虑贴现率$\beta$，则充电站年均建设投资成本$C_{con}$如公式（9）所示。

$$\begin{array}{l}{C}_{con}={\displaystyle \sum _i^S}\left[C{\left(i\right)}_{allland}+C{\left(i\right)}_{allpile}+C{\left(i\right)}_{other}\right]\times \\ \frac{\beta {\left(\mathrm{1}+\beta \right)}^{year}}{{\left(\mathrm{1}+\beta \right)}^{year}-\mathrm{1}}\end{array}$$

（9）

2.3.2　运营维护成本

充电站运营维护成本是指进行日常的设施维修养护、管理和清洁等产生的费用，假设其与充电站建设投资成本呈线性相关关系，如公式（10）所示。式中相关系数$\gamma$一般取0.01。

$$C_{ope}=\gamma \times C_{construction}$$

（10）

2.3.3　网络损耗成本

网络损耗成本是指电网电力运输过程产生的网络损耗费用，电网损耗主要来源于线损，可通过焦耳定律计算配电系统产生的电能损耗成本^［21］，如公式（11）所示。

$$C_{loss}=p_q{\displaystyle \sum _{t=\mathrm{1}}^T}{\displaystyle \sum _{i=\mathrm{1}}^l}{I}_{i,t}^{\mathrm{2}}{R}_i$$

（11）

式中：$l$为配电系统中的支路数； $I_{i,t}$为支路电流； $R_i$为支路电阻。

2.3.4　空载行驶成本

空载行驶成本是指车辆空载状态下前往充电站充电过程所产生的成本，包含耗电成本、车辆损耗成本和碳排放成本。其中，车辆损耗成本通过收入法^［22］计算得到，假设车辆在位置$i$产生充电需求，前往充电站$j$进行充电，则该过程产生的车辆损耗和耗电成本$C{\left(ij\right)}_{tear}$如公式（12）所示。

$$C{\left(ij\right)}_{tear}=p_{invest}\frac{d_{ij}}{v_{}}+p_{q}q{d}_{ij}$$

（12）

尽管Robotaxi使用阶段是零排放，但包括能源开采、加工、运输、储存等过程的电能生命周期碳排放不可忽略。从燃料生命周期的角度核算单位电能的碳排放系数，并计算车辆空载行驶碳排放成本^［23］，如公式（13）所示。

$$C{\left(ij\right)}_{carbon}=\eta \frac{d_{ij}qEF}{TE\times CE}\times \frac{\mathrm{12}}{\mathrm{44}}$$

（13）

综合公式（12）和（13），得到车队年均空载行驶成本$C{\left(ij\right)}_{emp}$如公式（14）所示。

$$\begin{array}{l}C{\left(ij\right)}_{emp}=\\ \mathrm{365}\left({\displaystyle \sum _i^S}{\displaystyle \sum _j^M}{N}_{ij}\left(\begin{array}{c}{p}_{invest}\frac{d_{ij}}{v_{}}+p_{q}q{d}_{ij}+\\ \eta \frac{d_{ij}qEF}{TE\times CE}\times \frac{\mathrm{12}}{\mathrm{44}}\end{array}\right)\right)\end{array}$$

（14）

2.3.5　订单损失机会成本

Robotaxi充电过程中会损失载客机会，而车辆所在区域的订单分布越密集，接单概率越高，则充电过程产生的机会成本越大。本研究基于车辆所在区域的接单概率$\epsilon$，通过车辆空载行驶时长估算车队的年订单损失机会成本$C_{opp}$，如公式（15）所示。

$$C_{opp}=\mathrm{365}\times p_{income}\epsilon {\displaystyle \sum _i^S}{\displaystyle \sum _j^M}{N}_{ij}\frac{d_{ij}}{v_{}}$$

（15）

2.3.6　车队排队充电成本

车辆到站时间和服务时间的随机性引起的车辆排队充电现象广泛存在，本研究基于M/M/C/∞排队系统考虑Robotaxi车辆排队充电过程。其中，车辆达站时间服从参数为$\lambda$的泊松分布，充电桩服务时间服从参数为$\mu$的负指数分布，车辆排队充电方式为先到先服务制，则车队排队充电产生的年等待成本$C_{que}$如式（16）所示。

$$C_{que}=\mathrm{365}{p}_{invest}{W}_q{\displaystyle \sum _i^S}{\displaystyle \sum _j^M}{N}_{ij}$$

（16）

2.3.7　目标函数

综合考虑充电站建设、运营和网络损耗成本，以及车队充电过程产生的空载行驶、订单损失和排队充电成本，得到综合成本最小化目标函数如公式（17）所示。式中，$\alpha$、$\beta$、$\gamma$、$\theta$、$\delta 、\mu$分别为各项成本的重要性系数，在实际制定选址布局方案时需要结合自身情况和重视程度设定。

$$\begin{array}{l}COST=\alpha C_{con}+\beta C_{ope}+\gamma C_{loss}+\\ \theta C_{emp}+\delta C_{opp}+\mu C_{que}\end{array}$$

（17）

2.3.8　约束条件

约束条件如式（18）—式（24）所示。公式（18）和（19）为研究区域内建设的充电站总数约束，以同时满足车队充电需求和保证充电站利用率；公式（20）为任意两个充电站之间的最小距离约束；公式（21）为每个充电站内配置的充电桩数量约束；约束（22）保证任意充电需求点到最近充电站的行驶距离须在电池阈值电量$SO{C}_{threshold}$所能承受的范围之内；约束（23）为研究区域内充电站服务能力上限；约束（24）为充电站总功率约束^［8］，即，当区域内存在多座充电站时，充电站总额定功率与区域内其他基础负荷的最大功率$P_{total,max}$之和不应超过电网允许接入的最大功率$P_{max}$。

$$N_{min}=\frac{Su{m}_{charge}}{nu{m}_{max}\mu T}$$

（18）

$$N_{max}=\frac{Su{m}_{charge}}{nu{m}_{min}\mu T}$$

（19）

$$d_{int}\ge d_{in{t}_{min}}$$

（20）

$$nu{m}_{min}\le nu{m}_i\le nu{m}_{max}$$

（21）

$$q{d}_{ij}\le SO{C}_{threshold}$$

（22）

$${\displaystyle \sum _i^S}{\displaystyle \sum _j^M}{N}_{ij}\le Serv{e}_{max}$$

（23）

$${\displaystyle \sum _{i=\mathrm{1}}^S}{n}_{i}c{h}_{}+P_{total,max}\le P_{max}$$

（24）

3 算法设计

充电设施选址布局是一个多目标优化问题，且涉及到诸多离散约束条件，可采用PSO算法进行求解。PSO算法是1995年由美国学者Eberhart和Kennedy提出的一种群体智能优化算法。粒子群由$n$维空间中的$m$个粒子组成的，每个粒子有适应度值、速度、位置三个属性。在本研究中，每个粒子均代表一种充电站布局方案，粒子的速度和位置是形如$\left[ln{g}_{\mathrm{1}},ln{g}_{\mathrm{2}},\cdots ,ln{g}_n,la{t}_{\mathrm{1}},la{t}_{\mathrm{2}},\cdots ,la{t}_n,nu{m}_{\mathrm{1}},nu{m}_{\mathrm{2}},\cdots ,nu{m}_n\right]$的向量，适应度用于衡量粒子当前位置的优劣程度。粒子根据适应度值不断更新当前速度和位置，从而逐渐逼近最优位置。粒子的速度和位置按照如公式（25）和（26）进行更新。

$$\begin{array}{l}{V}_{id}^{t+\mathrm{1}}=\omega V_{id}^t+c_{\mathrm{1}}{r}_{\mathrm{1}}\left(G_{best}{}_{id}^t-X_{id}^t\right)+\\ c_{\mathrm{2}}{r}_{\mathrm{2}}\left(P_{best}{}_{id}^t-X_{id}^t\right)\end{array}$$

（25）

$$X_{id}^{t+\mathrm{1}}=X_{id}^t+V_{id}^{t+\mathrm{1}}$$

（26）

式中：$\omega$表示惯性权重；$c_{\mathrm{1}}$和$c_{\mathrm{2}}$表示个体学习因子和社会学习因子；$r_{\mathrm{1}}$和$r_{\mathrm{2}}$表示随机速度扰动；$V_{id}^t$、$X_{id}^t、G_{best}{}_{id}^t、P_{best}{}_{id}^t$分别表示第$t$次迭代时第$i$个粒子的第$d$个维度的速度向量、位置向量、群体最优值和个体最优值。

3.1　改进机制1：基于遗传算子的种群更新机制

PSO算法种群多样性较差，容易陷入局部最优。遗传算法通过选择、交叉、变异等遗传算子实现对种群个体的优胜劣汰，种群更新方式更加多样，因而具有较好的全局寻优能力。但遗传算法的种群没有记忆，种群迭代速度较慢，对于高维度问题的求解收敛速度较慢。因此，本研究将遗传算子引入PSO算法，以提升种群多样性，进而提高算法全局寻优能力。在更新群体位置之前，计算适应度值以评估每个粒子在综合成本最小化目标下的性能。为确保较优个体具有更高的遗传概率，根据轮盘赌法则确定父代粒子。然后，根据概率$p_{\mathrm{1}}$选择一定数量的粒子进行两两交叉，并根据公式（27）—式（30）更新粒子的速度和位置。

$${\widehat{V}}_{id}^t=\frac{V_{id}^t+V_{jd}^t}{\left|\left|V_{id}^t+V_{jd}^t\right|\right|}\left|\left|V_{id}^t\right|\right|$$

（27）

$${\widehat{V}}_{jd}^t=\frac{V_{id}^t+V_{jd}^t}{\left|\left|V_{id}^t+V_{jd}^t\right|\right|}\left|\left|V_{jd}^t\right|\right|$$

（28）

$${\widehat{X}}_{id}^t=P_{\mathrm{c}}{X}_{id}^t+\left(\mathrm{1}-r\right)X_{jd}^t$$

（29）

$${\widehat{X}}_{jd}^t=P_{\mathrm{c}}{X}_{jd}^t+\left(\mathrm{1}-r\right)X_{id}^t$$

（30）

对于经过交叉操作后的每个粒子，根据变异率$P_m$，并结合高斯变异策略对粒子的位置按照式（31）—式（33）进行更新。

$$X_{id}^t=X_{id}^t+N\left(\mu ,{\sigma }^{\mathrm{2}}\right)$$

（31）

$$\mu =\frac{max\left(X_{id}^t\right)+min\left(X_{id}^t\right)}{\mathrm{2}}$$

（32）

$$\sigma =\frac{max\left(X_{id}^t\right)-min\left(X_{id}^t\right)}{\mathrm{6}}$$

（33）

3.2　改进机制2：基于适应度值的惯性权重自适应取值机制

惯性权重是影响PSO算法全局寻优能力的关键参数，较大权重有利于跳出局部最优，较小权重则能够加速全局收敛。因此，本研究基于粒子的速度和位置变化，并结合适应度值对惯性权重进行自适应动态调整，以兼顾算法的收敛速度和全局寻优能力。首先以较大值初始化惯性权重$\omega$，每次迭代后对当前粒子移动效果进行评估，并通过sigmod函数改变$\omega$，使得当粒子搜索到最优点附近时减小$\omega$^［24］，如式（34）所示。

$$\omega =sigmod\left(\alpha \frac{\mathrm{\Delta }h}{V_{id}^t}\right)$$

（34）

式中： $\alpha$为变化系数，$\mathrm{\Delta }h$为粒子当前适应度值与上次适应度值的差值。

综上，基于遗传算子和自适应惯性权重的IPSO算法求解充电站选址布局优化问题的流程如图3所示。

图3  IPSO算法求解流程

Fig.3  IPSO algorithm solution process

4 案例分析

本节首先基于出行轨迹数据提取乘客出行订单信息，然后基于乘客出行需求进行Robotaxi最小车队规模核算，进而根据最小车队规模进行车队充电需求仿真，最后，构建充电设施布局模型并分别采用PSO算法和IPSO算法进行求解，得到面向Robotaxi车队的充电设施布局方案。数据处理和模型编程基于Python语言，测试基于一台Intel（R） Core（TM） i5-8250U CPU @ 1.60GHz 1.80 GHz、内存8GB的笔记本电脑和一台Inter（R） Xeon（R） Gold 6248R CPU @ 3.00GHz 2.99GHz（2处理器）、内存192GB的服务器。

4.1　数据处理

数据集为滴滴出行平台在四川省成都市成华区、青羊区、金牛区2016年11月1日至2016年11月30日产生的车辆GPS轨迹数据，共计超过10亿条信息。本研究所用数据特征共包含五列，分别为司机编号、乘客编号、时间、经度和纬度，其中时间为Unix格式，如图4所示。

图4  原始数据结构示例

Fig.4  Example of original data structure

研究区域内的乘客出行轨迹分布如图5所示，靠近市中心的西南方向的乘客轨迹相较于远离市中心的东北方向更加密集，表明乘客在此区域出现的频率更高。为了进一步定位乘客出行起讫点，根据Asamer等^［25］的建议，本研究采用正六边形蜂窝结构对研究区域进行完整镶嵌划分。综合考虑乘客的步行范围和计算机运算能力，确定六边形边长为500米，最终将研究区域划分为120个六边形网格，如图5所示。基于此，将研究区域内所有乘客出行起讫点聚类至蜂窝网中心，即当乘客的起讫点在某一六边形中时，则用该六边形的中心点代表乘客出行的起讫点。

图5  车辆GPS轨迹分布及研究区域网格划分

Fig.5  Vehicle GPS track distribution and gridding of the study area

4.2　实例验证

4.2.1　最小车队规模核算

基于前文介绍的最小车队规模计算流程，以2016年11月1日研究区域内产生的乘客出行订单为例（共包含由18 863辆车服务的172 651个出行订单），构建有向无环图并计算行程衔接匹配方案，Hopcroft-Karp算法表明可以利用10293条最小不相交路径将所有订单进行连接，因此，为满足乘客出行需求所需要的Robotaxi最小车队规模为10 293辆，这相较于实际车队规模减少了45.43%（见表3）。

表3  最小不相交路径结果示例

Tab.3  Example of minimum disjoint path result

编号	最小不相交路径节点轨迹
1	［9， 51， 291， 501， 947， 1961， 2748， 3971］
2	［40， 76， 163， 539， 614， 931， 1629， 2390， 2836， 4116， 4991， 5127］
3	［32， 89， 350， 464， 933， 1227， 1995， 2554， 2789， 3311， 3647， 3908， 4696， 5316， 5567， 5779， 6197， 6407， 6680］
︙	︙
7164	［4， 164， 349， 943， 1646， 2086， 3002， 5048， 5347， 5816， 6605， 6793， 7050， 7292， 7640， 7991， 8149， 8549， 8758， 9277， 9586， 9819， 10040， 10364， 10635， 10860， 10981， 11803， 12186， 13049， 13620， 15201， 17639， 19085， 21658， 24126， 25768， 28780， 34331， 35704， 37089， 40759， 43716， 45365， 47306， 50842， 51858， 53499， 55015， 56158， 59019， 66061， 69570， 71620， 74068， 77265， 79541， 80325， 88560， 91734， 95946， 99578， 102914， 104612， 105260， 107127， 108925， 111510， 113281， 115358， 116815，118971， 126153， 128939］

4.2.2　车队充电需求分析

当前成都市电动出租车多为吉利帝豪EV450，配置为电池容量52 kWh，续航为300 km，快充1 h。据此，本研究设定Robotaxi车辆相关参数为：车辆电池容量为5 0 kWh，充电阈值$SO{C}_{threshold}$为20%，车辆单位里程耗电$q$为0.2 kWh·km^-1，电池充电速度$ch$为1.2 kWh·min^-1。基于上文计算的最小车队规模，考虑实际运营过程车队需要有一定冗余，本研究假定投放的车队规模为15 000辆，进而按照充电需求仿真流程（见图1）记录研究区域内产生的充电需求，结果如表4所示。以表4第一行数据为例进行说明，该辆车于21时在位置92产生充电需求，所需充电电量为28.676 kWh。

表4  充电需求数据示例

Tab.4  Driver transportation behavior feature variables

位置编号	经度	纬度	产生 时刻/h	容量需求/（kWh）
92	104.115 0	30.719 282	21	28.676
112	104.085 0	30.719 282	20	27.314
114	104.092 5	30.723 612	21	28.956
117	104.100 0	30.650 000	22	31.054
16	104.092 5	30.654 330	22	31.228
35	104.062 5	30.662 990	24	35.176
…	…	…	…	…

车队充电需求在研究区域内的空间分布如图6所示，深色区域表示充电需求较强，浅色区域充电需求较少。可以看出充电需求在研究区域内的分布具有明显的空间异质性：在左中和右下位置形成两个较为明显的充电需求旺盛区，而在远离市中心的右上区域则充电需求相对较低，这与乘客出行轨迹分布具有一定的相似性。考虑到充电需求在空间的分布并不均匀，因此结合充电需求的时空分布进行充电设施的选址定容优化是必要的。

图6  充电需求空间分布热力图

Fig.6  Heat map of spatial distribution of charging demand

如图7所示，车队充电需求随时间呈现明显的波动性：从00：00时起，车队充电需求从中等水平逐渐降低，至凌晨05：00时左右降至全天最低谷，随后快速上升，至14：00时左右达到全天最高峰，之后出现回落，在15：00时至19：00时期间维持一个相对较高的充电需求，之后再次快速上升，在21：00时左右第二次到达全天充电需求的最高峰，之后呈持续下降趋势直至结束。

图7  充电需求的全天分布

Fig.7  All-day distribution of charging demand

4.2.3　站点选址定容方案

为计及土地价格对充电设施选址布局的影响，根据成都市中心城区商服用地土地基准地价（见表3），将研究区域划分成Ⅰ、Ⅱ、Ⅲ和Ⅳ四个级别，如图8所示，每个级别对应的地价分别为38 850元/m²、26 700元/m²、18 765元/m²和12 450元/m²。此外，本研究对目标函数中各项成本同等看待而不偏向于特定主体，因此各项成本系数取值均为1。

图8  区域地价划分结果

Fig.8  Results of regional land value classification

表3  成都市中心城区商服用地基准地价

Tab.3  Benchmark land price for commercial land use in Chengdu city center

级别/基准地价	元/m2	万元/亩
Ⅰ级	38850	2590
Ⅱ级	26700	1780
Ⅲ级	18765	1251
Ⅳ级	12450	830
Ⅴ级	7650	510
Ⅵ级	5550	370
Ⅶ级	4200	280

注：  1亩=666.667 m²。

首先采用PSO算法，计算建设不同数量充电站的综合成本，如图9所示。可以发现：①随着充电站建设数量的增加，综合成本呈现先下降后上升的趋势，并在充电站建设数量为17个时，综合成本取得最小值6 881.865 4万元，其中建设投资成本占比达75%，此时的充电站位置分布和容量配置如图10所示。②充电站建设投资成本整体上随充电站数量的增加呈增加趋势，这意味着建设的充电站数量增加时，前期需要投入更多成本。然而应该指出的是，这种趋势并不严格时时成立，如当充电站的数量从13增加至17时，建设投资成本反而逐渐降低，其主要原因是充电站建设数量增加的过程中，建设的充电桩总数可能得到优化而减少，这使得总建设投资成本得到降低。③其他成本随充电站建设数量的增加而减少，这是因为随着充电站建设数量的增加，运营车队的充电需求能够得到更便捷的满足，因而车队空载行驶成本、排队充电成本和订单损失成本均得到降低，尽管其中的运营维护成本会增加，但由于其占其他成本的比例较小，因而并没有改变整体变化趋势。④建设投资成本占综合成本比例随充电站数量的增加呈略微增加趋势。

图9  基于PSO算法的充电站布局成本优化

Fig.9  Charging station layout cost optimization based on PSO algorithm

注：  其他成本=运营维护成本+车队空载行驶成本+车队排队充电成本+订单那损失机会成本+网络损耗成本。

图10  基于PSO算法的充电站选址定容方案

Fig.10  PSO algorithm based charging station siting and capacity determination scheme

为对比验证本研究所提出的基于遗传算子和自适应惯性权重的IPSO算法的有效性，进一步采用IPSO法计算建设不同数量充电站时的综合成本，如图11所示。由图11可以看出：① 采用改进PSO算法进行求解时，随着建设充电站数量的增加，综合成本呈现先下降，再上升的变化趋势。当充电站建设数量为17时，综合成本取得最优值6 378.751 6万元，这相比PSO算法降低了7.31%。② 对比图9和图11可以发现，建设不同数量的充电站，IPSO算法相比PSO算法的综合成本均有不同程度的降低，这表明遗传算子和自适应惯性权重机制的引入能够有效提升算法全局寻优能力。③ 此外，尽管两种算法所输出的充电站最优数量均为17，但具体的位置分布和容量配置有所不同。

图11  基于IPSO算法的充电站布局成本优化

Fig.11  Cost optimization of charging station layout based on IPSO algorithm

基于IPSO算法的充电设施选址布局方案如图12所示。由图12可以看出，基于IPSO算法的充电站布局方案在区域Ⅰ、Ⅱ、Ⅲ、Ⅳ建设的充电站数量分别为0、2、1和14，而基于PSO算法的充电站布局方案在区域Ⅰ、Ⅱ、Ⅲ、Ⅳ建设的充电站数量分别为1、1、1和14，进一步对比图10和图12可知，IPSO算法通过将原本位于区域Ⅰ的一个充电站转移至区域1与区域2的交界处（位于区域2），进而能够在有效降低充电站建设投资成本的同时，满足各区域的充电需求。此外，基于IPSO算法的充电站设施布局方案的充电桩建设总数为381，这相比PSO算法的最优充电桩数量403降低了5.46%。

图12  基于IPSO算法的充电站选址定容方案

Fig.12  IPSO algorithm-based charging station siting and capacity determination scheme

综上，遗传算子和自适应惯性权重的引入能够有效提升了算法的全局寻优能力，有效降低充电站建设综合成本。在本文案例分析中，IPSO算法相比PSO算法使得充电站综合投入成本得到降低主要表现在两方面，即一方面使得充电站的位置分布得到优化，另一方面使得建设的充电桩总数得到降低。

5 结论

本研究以乘客出行数据为基础，综合考虑充电需求、车队规模和成本因素，构建了以综合成本最小化为目标的充电设施布局优化模型，并针对模型求解提出一种基于遗传算子和自适应惯性权重的IPSO算法，最后，基于成都市乘客出行轨迹数据对所提模型和算法展开了案例验证，最终得到以下主要结论：

（1）充电需求分布、车队规模配置和综合投入成本是充电设施布局优化的关键考量因素。

（2）Robotaxi车队的充电需求具有典型的时空异质性，从时间维度来看，充电需求在凌晨5时达到全天最低谷，并分别在14时和21时两次迎来全天最高峰。

（3）在基于成都市的案例分析中，当建设的充电站数量为17个时，综合成本取得最小值6356万元，其中建设投资成本占比最高，达74%。

（4）本研究所提出的IPSO算法相比PSO算法使得充电站布局方案综合成本降低7.31%，这表明遗传算子和自适应惯性权重的引入可有效改善粒子种群多样性，提高算法寻优能力。

对于未来的研究，一方面可以考虑多种能源补给方式，如将换电模式和充电模式结合起来，制定多样化充电设施布局方案，以满足不同场景的车队补电需求；另一方面，可以扩大研究区域范围，并在不同地区开展案例分析，以进一步验证本文所提出的优化模型和求解算法的普适性。