共享自行车系统(bike sharing system, BSS)在世界各地得到了广泛的应用, 目前分为锁桩式BSS(公共自行车系统)和无锁桩式BSS(共享单车).锁桩式BSS系统通过固定自助服务站提供租还服务, 共享单车通过电子围栏方式来解决无序停放问题.无论是固定服务站还是电子围栏, 都有停放自行车容量的限制, 由于交通出行的潮汐时空分布不均衡特性, 租还车难将成为BSS发展中突出的长期问题.通过调度车辆对BSS进行再平衡, 可以有效地缓解租还车难的现象.

目前, 国内外对共享自行车系统调度问题的研究成果较多, 如董红召等^[1]提出城市公共自行车系统的滚动时域调度算法, 基于大样本历史数据对公共自行车服务点自然租赁需求进行估算及验证^[2].Salaken等^[3-5]对公共自行车站点需求量进行预测, 以降低调度成本.Adham等^[6-8]对BSS调度方案提出改进的调度算法, 优化调度路径, 提高调度效率, 但这些成果缺少调度区域划分的研究.尤其是大型BSS系统, 若不采用分区调度, 则往往算法求解耗时长, 调度实时性差, 成本增加, 因此有必要对BSS调度区域进行科学划分, 实行分区调度.

研究BSS调度分区的文献较少, Come等^[9]根据公共自行车站点的使用情况进行聚类分区, 分析用户出行特征与站点地理位置之间的关系, 但未涉及到BSS调度区域的划分.陈昕昀等^[10]利用K-means聚类划分调度区域, 使得各站点离调度中心的平均距离最短, 但没有考虑各区域调度需求量的平衡.董红召等^[11]基于服务点之间的关联特性进行调度区域的划分, 但划分的每个调度区域比较笼统, 不同区域之间出现不平衡后不能进行协调互补.徐建闽等^[12]的多层次分区调度方法虽然采用了区域划分和分层调度, 但调度小区的划分停留在经验和行政区划阶段, 缺少区域内自行车租/还自平衡的考虑.

本文针对上述问题, 提出基于分形树的自平衡调度区域划分模型^[13], 即根据分形树的自相似性特征, 从叶子节点开始, 将租/还需求互补的节点层层聚类形成多级分形自平衡区域, 使得每层调度区域的需求尽量达到自平衡来减少跨区调度.根据不同时段调度需求量的动态变化, 提出基于周转率杠杆系数的共协矩阵聚类融合算法.对聚类结果进行优化融合^[14], 提高了分区结果的鲁棒性.

1 基于分形树的自平衡调度区域划分模型

BSS自行车处于分区、分级调度模式下, 调度车辆首先负责所属区域内的自行车调度.只有当区域内租还总量差距较大、无法平衡时, 开展跨区域间调度.为了尽量减少区域间调度的次数, 划分的每一个调度区域都需要能够尽量达到租/还自平衡, 因此在开展调度区域划分时, 将在某时间区间内具有自行车调入需求的服务点/区域和调出需求的服务点/区域放入同一个上级区域, 这样调度车辆可以在这两个不同需求的服务点/区域间进行自行车运送, 从而达到该级区域内的自平衡.

BSS自平衡分级调度区域划分模型如图 1所示.先划分一级自平衡区域S¹, 按照服务点的自行车调入和调出需求的互补关系, 将距离较近且有需求互补关系的服务点划在同一区域形成自平衡区域S¹；再将具有调入和调出需求的互补区域S¹合并形成二级自平衡区域S², 以此类推, 最终形成最高级自平衡区域, 即整个BSS系统.定义为：对于一个H级分区调度BSS, 每个级别区域定义为Sⁿ(n=0, 1, 2, …, H), 其中, 当n=H时, 最高级S^H表示整个BSS, BSS已达到自平衡；当n=0时, 最低级S⁰表示服务点, 服务点平衡是调度系统服务的最终目标.

分形理论是非线性科学的一种方法论和认识论, 它可以依靠简单结构与规则的迭代来生成复杂系统, 这种系统具有鲜明的整体与局部结构的自相似特性^[12].显然, 上述划分的BSS自平衡区域各个层级的组织和调度方法具有自相似性, 即鲜明的分形结构特性, 因此, 图 1所示的BSS调度模型可以用一个自平衡分形树表达, 如图 2所示.分形树中, 每个叶子节点表示一个自行车服务点S⁰, 枝节点表示某一级的自平衡区域Sⁿ(0 < n < H), 根节点表示整个BSS系统S^H, 各个层级的调度区域结构和规则根据分形系统的自相似性来生成.图 2中, 所有叶子节点根据地理位置和需求互补关系, 可以形成多个叶子级调度区域；这些叶子级调度区域是形成上一级枝节级调度区域的枝节点, 形成的枝节级调度区域将作为上一级枝节调度区域的枝节点；采用自形似方式层层迭代, 最终可以构成根级调度区域.

BSS自平衡分形树区域的生成与划分采用如下关键参数.

1) 叶子节点i的不平衡度W_i(τ)定义为在时间段τ内, 服务点i借出和还入的自行车数量的差值.

$ {W_i}\left( \tau \right) = Z_i^{{\rm{out}}}\left( \tau \right) - Z_i^{{\rm{in}}}\left( \tau \right). $

(1)

式中：Z_i^out(τ)和Z_iⁱⁿ(τ)分别为服务点i在时间段τ的借车和还车总数.可以看出, 当W_i(τ)>0时表示自行车借出量大于归还量, 若W_i(τ)>0持续较长时间, 则服务点将进入空位状态, 借车困难；当W_i(τ) < 0时表示自行车借出量小于归还量, 若W_i(τ) < 0持续较长时间, 则服务点进入满位状态, 还车困难.

2) 枝节点α的不平衡度W_α(τ)定义为在时间段τ内, 区域α中所有服务点不平衡度的总和.

$ {W_\alpha }\left( \tau \right) = \sum\limits_{i \in \alpha } {{W_i}\left( \tau \right)} . $

(2)

3) 叶子节点i、j之间的互平衡强度EW_{i, j}(τ).如果服务点i与j间的距离越近、W_i(τ)+W_j(τ)的绝对值越小, 那么这两个节点的互平衡关系越强.D_{i, j}表示服务点i和j之间的距离, 因此EW_{i, j}(τ)的计算可以设计为

$ {\rm{E}}{{\rm{W}}_{i, j}}\left( \tau \right) = \frac{1}{{\left| {{W_i}\left( \tau \right) + {W_j}\left( \tau \right)} \right|\gamma + {D_{i, j}}}}. $

(3)

式中：γ为不平衡度的距离效应转换常数.

4) 枝节点α、β之间的互平衡强度EW_{α, β}(τ).根据分形自相似性可知, 与叶子节点算法相似, 枝节点区域α、β之间的距离越近且不平衡度总和越小, 则说明这两个区域的互平衡关系越强, 计算公式为

$ {\rm{E}}{{\rm{W}}_{\alpha , \beta }}\left( \tau \right) = \frac{1}{{\left| {{W_\alpha }\left( \tau \right) + {W_\beta }\left( \tau \right)} \right|\gamma + {D_{\alpha , \beta }}}}. $

(4)

式中：D_{α, β}为两区域中心点间的距离.区域α的中心点坐标(x_α, y_α)计算公式为

$ {x_\alpha } = \frac{{\sum {_{i \in \alpha }} {x_i}}}{{N\left( \alpha \right)}}, \;\;\;{y_\alpha } = \frac{{\sum {_{i \in \alpha }} {y_i}}}{{N\left( \alpha \right)}}. $

(5)

式中：N(α)为区域α中的服务点个数, (x_i, y_i)为区域α中服务点i的位置坐标.区域β的中心点坐标计算方式相同.

2 BSS调度的分形树自平衡区域划分算法FSPA

构造实现分形树自平衡调度区域划分模型的算法(fractal tree based self-balanced partitioning algorithm, FSPA).FSPA包括调度区域范围设计、基于互平衡强度的动态聚类及引入BSS周转率杠杆的共协矩阵自平衡区域聚类融合.

2.1 考虑快速服务响应的叶子与枝节级调度区域范围设计

调度区域范围设计直接影响BSS的快速响应服务特性, 若叶子级区域范围太大, 则会降低调度响应速度；区域越小, 该区域自平衡往往越困难.因此, 需要确定叶子级调度区域的合理范围.

设定调度车在非服务期间停留在该调度区域的中心位置, 以便及时响应调度需求.设叶子级调度区域覆盖半径为R, 调度车平均行驶速度为v, 调度响应时间为[δ^low, δ^upp], 则在最小响应时间δ^low内调度车能够到达的距离可以作为叶子级调度区域的最小半径R_min, 理想情况下R_min=δ^lowv, 但调度车沿途经过等待调度的服务点时会停留服务.设每个服务点的平均停留服务时间为t, 道路上站点分布密度为ρ, 极限情况下每个沿途服务点都是需调度服务点, 则R_min的计算方程为

$ {R_{\min }} = ({\delta ^{{\rm{low}}}} - {R_{\min }}\rho \bar t)\bar v. $

(6)

求解式(6), 可得

$ {R_{\min }} = \frac{{{\delta ^{{\rm{low}}}}\bar v}}{{1 + \rho \bar t\bar v}}. $

(7)

同理可得R_max, 因此叶子级调度区域面积为[S_min¹, S_max¹], 计算公式为

$ \left. \begin{array}{l} S_{\min }^1 = {\rm{ \mathsf{ π} }}{({R_{\min }})^2} = {\rm{ \mathsf{ π} }}{\left( {\frac{{{\delta ^{{\rm{low}}}}\bar v}}{{1 + \rho \bar t\bar v}}} \right)^2}, {\rm{ }}\\ S_{\max }^1 = {\rm{ \mathsf{ π} }}{({R_{\max }})^2} = {\rm{ \mathsf{ π} }}{\left( {\frac{{{\delta ^{{\rm{upp}}}}\bar v}}{{1 + \rho \bar t\bar v}}} \right)^2}. \end{array} \right\} $

(8)

BSS服务点发出调度需求到获得调度服务的理想时间范围、即[δ^low, δ^upp]与服务级别的对应关系如表 1所示.

对于枝节级调度, 区域的划分需要考虑该层级区域内部的租还需求自平衡和区域面积范围.根据实践经验可知, 上级调度区域负责3~5个下级调度区域, 因此n级调度区域面积范围[S_minⁿ, S_maxⁿ]的计算公式为

$ \left. \begin{array}{l} S_{\min }^n = 3S_{\min }^{n - 1}, \\ S_{\max }^n = 5S_{\max }^{n - 1}. \end{array} \right\} $

(9)

2.2 叶子级自平衡区域的互平衡强度动态聚类

根据某时段τ的各服务点租还数据, 可以计算服务点间的互平衡强度, 开展该时段叶子级调度区域的动态聚类, 算法如下.

1) 将BSS内所有的服务点放入集合C⁰, 服务点数量为N, 这些服务点是分形树的叶子节点.设叶子级调度区域集合为C¹.

2) 计算集合C⁰中所有节点间的互平衡强度, 找出每个节点对应的最大互平衡强度节点组成互补节点对, 如节点i的最大互平衡强度为EW_{i, j}(τ), 则i和j组成互补节点对.

3) 在这些互补节点对中, 去除其中互平衡强度小于平均值的互补节点对；若剩余节点对中有些节点相交, 则去除其中互平衡强度较小的节点对.

4) 将剩余的节点对形成枝节点并计算面积, 枝节点α的面积S_α计算方法为

$ \left. {\begin{array}{*{20}{c}} {{S_\alpha } = [\max ({x_i}) - \min ({x_i})][\max ({y_i}) - \min ({y_i})];}\\ {i \in \alpha .} \end{array}} \right\} $

(10)

若满足S_α>S_min¹, 则将该枝节点放入叶子级调度区域集合C¹；若不满足S_α>S_min¹, 则用α替换集合C⁰中α包含的节点.

5) 重复2)~4), 直到集合C⁰中没有剩余节点.

集合C¹中枝节点是时间段τ的叶子级自平衡区域聚类结果.对于BSS的某一服务周期T, 该周期内可以得到T/τ个聚类结果, 而且不同时段的聚类结果对调度的影响程度不同, 因此需要实现多个聚类结果的融合, 提出导入周转率杠杆的共协矩阵聚类融合改进算法.

2.3 导入周转率杠杆的共协矩阵聚类融合改进算法

将BSS周转率杠杆导入共协矩阵聚类融合算法, 实现多个聚类结果的融合.令m=T/τ为服务周期T内的时段个数, τ_k表示第k个时段(k=1, 2, …, m), N为BSS的服务点数量, 则BSS第k时段的自行车周转率ω_k的计算公式为

$ {\omega _k} = \frac{{\sum {_{i \in {\rm{BSS}}}} Z_i^{{\rm{out}}}({\tau _k}) + \sum {_{i \in {\rm{BSS}}}} Z_i^{{\rm{in}}}({\tau _k})}}{N}. $

(11)

生成N×N维共协矩阵U, 矩阵的元素u_{i, j}表示服务点i和j共同出现在同一个聚类区域中的概率和对应时段周转率的乘积.设置阈值, 若u_{i, j}大于阈值, 则认为服务点i与j属于最终聚类结果中的同一区域.m个不同时段聚类结果形成的序列集合为L={L_k, k=1, 2, …, m}, 其中每一个聚类结果包括多个叶子级调度区域L_k={C_k¹(l), l=1, 2, …, σ_k}, σ_k为第k个时段聚类出的叶子级调度区域个数.U定义为

$ \mathit{\boldsymbol{U}} = {[{u_{i, j}}]^{N \times N}} = \left[ {\begin{array}{*{20}{c}} {{u_{1, 1}}}& \cdots &{{u_{1, N}}}\\ \vdots &{}& \vdots \\ {{u_{N, 1}}}& \cdots &{{u_{N, N}}} \end{array}} \right], $

(12)

$ {u_{i, j}} = \frac{1}{m}\sum\limits_{k = 1}^m {{\omega _k}{\delta _k}\left( {i, j} \right)} . $

(13)

式中：δ_k(i, j)为阶跃函数, 表示服务点i和服务点j在k时段是否属于同一个叶子级调度区域, 判别公式为

$ {\delta _k}\left( {i, j} \right) = \left\{ {\begin{array}{*{20}{c}} {1, }&{i, j \in C_k^1\left( l \right);}\\ {0, }&{i \in C_k^1\left( l \right), j \notin C_k^1\left( l \right).} \end{array}} \right. $

(14)

以周转率均值和平均聚类概率为判断阈值θ, 计算公式为

$ \theta = \frac{{\sum {_{k = 1}^m} {\omega _k}}}{m} \times 50\% . $

(15)

u_{i, j}>θ, 则认为服务点i、j属于同一个叶子级调度区域, 据此可以形成多个新的聚类区域.若有面积小于S_min¹的区域, 则将该区域划分到与其地理位置最近的区域中；最后可以形成服务周期T的叶子级调度区域划分结果：L_new={C_new¹(l), l=1, 2, …, σ_new}.

2.4 枝节级调度区域的构建

在BSS实际运营过程中, 可能会出现某叶子级区域的服务点自行车都被单向大量借出或还入, 早晚高峰时最容易出现, 这需要在不同的叶子级区域之间进行自行车调配, 称为枝节级区域调度, 即叶子级区域内部出现不平衡时需要实施跨区域调度.对于n级自平衡调度区域的构建, 算法如下.

1) 令n=2, 将聚类融合后的叶子级调度区域作为枝节点放入集合C_new^n-1.

2) 令k=1, τ_k表示服务周期T内第k个时段.

3) 计算集合C_new^n-1中枝节点在τ_k时段的互平衡强度EW_{α, β}(τ_k), 找出每个枝节点对应的最大互平衡强度节点组成互补枝节点.

4) 在互补枝节点集合中, 删除其中互平衡强度小于平均值的互补枝节点和相交节点对中互平衡强度较小的互补枝节点.

5) 剩余互补枝节点将成为新的枝节点, 计算区域面积, 计算方法如式(10).若新枝节点面积＞S_minⁿ, 将该枝节点放入n级调度区域集合C_kⁿ；若不满足条件, 则替换集合C_new^n-1中被聚类的节点.

6) 重复3)~5), 直到集合C_new^n-1中没有剩余节点.集合C_kⁿ中的枝节点是τ_k时段的n级自平衡区域.

7) k=k+1, 重复3)~6), 直到周期T内最后一个时段的n级自平衡区域聚类完成.

8) 针对BSS在不同时段聚类产生的n级调度区域C_kⁿ, 采用周转率杠杆共协矩阵聚类融合算法生成最终的n级调度区域聚类结果, 将最终聚类结果放入集合C_newⁿ.

9) n=n+1, 若S_BSS∈[S_minⁿ, S_maxⁿ], 则结束聚类；否则, 重复2)~9).

3 BSS分形树自平衡调度区域划分的实验与分析 3.1 BSS自平衡区域划分实验

以杭州锁桩式BSS下沙地区185个服务点为例, BSS服务级别取二级服务, 对分形树自平衡调度区域的划分方法进行实验验证.

根据杭州BSS历史数据, 调度车每完成一个服务点调度平均行驶0.9 km, 平均调入(或调出)约11量自行车, 因此, 式(3)中γ取0.9/11=0.08.通过式(3)可以计算所有服务点间的互平衡强度, 如表 2所示为某工作日第一个时段(6:00-7:00)的185个服务点与对应的最大互平衡强度服务点信息列表.表中, i为服务点编号, W_i(τ)、W_j(τ)为不平衡度, j为互补服务点, D_{i, j}为服务点间距, EW_{i, j}(τ)为服务点间互平衡强度.

表 2中, 互平衡强度的均值为2.89, 互平衡强度＞2.89的互补叶子节点可以进行聚类区域划分.由历史数据可得, 调度机动车平均行驶速度为20 km/h, 每个服务点的平均服务时间为5.5 min, 线路上站点平均分布密度为2.8个/km, 如BSS服务级别取二级服务, 则根据式(8)可知, 叶子级调度区域面积为[3.71, 8.35].根据该面积范围及叶子级自平衡区域划分的动态聚类算法可知, 6:00-7:00时段的叶子级聚类结果如图 3所示.

图 3中, 圆型图标表示服务点不平衡度W_i(τ) < 0, 需要调出自行车；方型图标表示W_i(τ)>0, 需要调入自行车, 颜色越深则不平衡度绝对值越大.第一轮互补叶子节点聚类结果如图 3的实线矩形框所示.可以看出, 实线框中的聚类节点距离相近, 而且这些节点的自行车调出调入需求相反, 可以互补形成自平衡枝节点.继续计算新枝节点和剩余叶子节点的互平衡强度, 最终聚类形成的叶子级调度区域如图 3的虚线框所示.

如图 4所示为该地区早晚高峰时段互平衡服务点聚类图.可以看出, 在不同时段获取的叶子级自平衡区域聚类结果不同, 因此在获取6:00-22:00之间16个时段的不同聚类结果后, 可以采用引入周转率杠杆的共协矩阵聚类融合算法进行聚类融合.

如表 3所示为根据式(11)计算出的一天的服务时间中各时段的周转率ω_k, 以此作为聚类融合的杠杆系数, 代入各时段聚类结果, 可以建立共协矩阵如下：

$ \begin{array}{l} \mathit{\boldsymbol{U}} = {[{u_{i,{\rm{ }}j}}]^{185 \times 185}} = \\ \left[ {\begin{array}{*{20}{c}} {167.56}&{15.50}&{25.51}& \cdots &0\\ {15.50}&{167.56}&{120.55}& \cdots &0\\ \vdots & \vdots & \vdots &{}& \vdots \\ 0&0&0& \cdots &{167.56} \end{array}} \right]. \end{array} $

根据式(15)计算θ=83.78, 可得融合的叶子级自平衡调度区域C_new¹.如表 4所示为C_new¹各区域服务点及在6:00-7:00的区域不平衡度, 其中共有12个区域, 区域不平衡度为该聚类区域所有服务点的不平衡度总和.可以看出, 有部分区域的不平衡度的总和为0或者较小的数值；还有部分区域的不平衡度较大, 需要跨区域进行调度.

根据每个叶子级区域中心点间的距离和不同时段的不平衡度计算区域互平衡强度, 再采用枝节级区域划分聚类算法, 可得融合的二级自平衡调度区域C_new².以此类推, 最终的三级聚类结果如图 5所示.

图 5中, 下沙地区三级自平衡调度区域划分结果考虑了各服务点及同级调度区域间的需求互补, 将具有自行车调入和调出需求的互补服务点/区域划入同一上级区域, 以达到区域内的调度自平衡.下面对该划分结果进行实验验证和分析.

3.2 BSS自平衡区域划分结果的验证与分析

选择下沙地区BSS调度历史数据进行分析, 如表 5所示为该区域在某工作日早高峰9:00的调度需求量.表中, 正数表示需要调入的自行车数量, 负数表示需要调出的自行车数量, 括号内为需要调度的服务点, 共有62个站点需要调度, 它们分别归属于不同叶子级和枝节级调度区域, 如表 5所示.

从表 5可以看出, BSS两个叶子级调度区域C_new¹(6)、C_new¹(7)中没有需要调度的服务点, 不需要调度；若叶子级区域调度需求量小于区域服务点个数, 则不参与跨区调度, 因此C_new¹(3)、C_new¹(9)和C_new¹(12) 3个叶子级区域只在自己区域内部平衡；其他叶子级区域需要在枝节级区域内部进行跨区自行车调度.由此计算可得枝节级区域C_new²(1)、C_new²(2)与C_new²(3)的调度需求分别为51、13和-45.当BSS的枝节级区域无法自平衡时, 需要跨区调度, 方法与叶子级区域的平衡调度类似；因为C_new²(2)调度需求量小于区域服务点个数, 只在自己区域内部平衡, C_new²(1)与C_new²(3)两个枝节级区域需要跨区调度, 即在根级区域开展自行车平衡调度.

由历史数据的实验分析可知, 采用BSS分形树自平衡调度区域划分方法, 可以使得各级区域内部服务点的租借与归还趋向于平衡.当出现较小不平衡量时, 可以区域内部消化, 例如表 5的3个叶子级区域和1个枝节级区域；当只有出现较大不平衡时, 开始跨区域调度, 这样不仅有助于减少自行车调度的数量, 还可以较大程度地减少调度车的行驶距离.

将研究结果实验应用于杭州下沙BSS的运营管理, 实验效果良好.结果表明, 研究的BSS分形树自平衡区域划分方法切实可行, 既能够较好地提升BSS运行的调度响应速度和服务质量, 又充分利用BSS系统自行车租借归还的空间相关性, 尽量使得各个区域达到自平衡, 从而降低服务成本.

4 结语

科学合理的区域划分和分区调度是大型城市BSS提高车辆调度效果和服务质量的有效方法.提出分形树自平衡调度区域划分模型, 基于服务点租/还需求的互平衡特性进行聚类, 同时采用周转率杠杆共协矩阵聚类融合算法, 分别形成叶子级和枝节级自平衡调度区域.杭州下沙城区共享自行车系统实验结果表明, 自平衡分区算法能够以较低成本满足BSS系统的调度需求.

研究的方法可以推广应用于共享单车, 即无锁桩式BSS.一些城市已经开始试用网格化调度管理共享单车系统, 工作人员各自负责所属网格区域内自行车的有序停放和调度.由于按照行政区划或者经验划分区域, 使得管理效率不高.对于无锁桩式BSS、尤其是电子围栏式BSS的调度分区, 因为这类BSS具有海量出行数据以及与锁桩式BSS相同的调度需求, 采用分形树自平衡调度区域划分模型和方法对管理网格进行划分, 更加有助于提高响应速度和工作效率, 降低调度成本, 从而提升城市管理水平和BSS服务质量.