CN115936240B - 共享单车需求预测与投放调度方法 - Google Patents
共享单车需求预测与投放调度方法 Download PDFInfo
- Publication number
- CN115936240B CN115936240B CN202211665231.5A CN202211665231A CN115936240B CN 115936240 B CN115936240 B CN 115936240B CN 202211665231 A CN202211665231 A CN 202211665231A CN 115936240 B CN115936240 B CN 115936240B
- Authority
- CN
- China
- Prior art keywords
- station
- site
- demand
- stations
- particles
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 60
- 238000012384 transportation and delivery Methods 0.000 title claims abstract description 24
- 238000003066 decision tree Methods 0.000 claims abstract description 38
- 238000005192 partition Methods 0.000 claims abstract description 18
- 238000012549 training Methods 0.000 claims abstract description 13
- 238000005457 optimization Methods 0.000 claims abstract description 7
- 239000002245 particle Substances 0.000 claims description 108
- 230000006870 function Effects 0.000 claims description 25
- 230000008569 process Effects 0.000 claims description 22
- 238000012937 correction Methods 0.000 claims description 15
- 238000004422 calculation algorithm Methods 0.000 claims description 10
- 230000035772 mutation Effects 0.000 claims description 9
- 238000012546 transfer Methods 0.000 claims description 9
- 229910052799 carbon Inorganic materials 0.000 claims description 7
- 238000001556 precipitation Methods 0.000 claims description 7
- 230000002093 peripheral effect Effects 0.000 claims description 6
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 230000003044 adaptive effect Effects 0.000 claims description 3
- 230000002068 genetic effect Effects 0.000 claims description 3
- 238000003780 insertion Methods 0.000 claims description 3
- 230000037431 insertion Effects 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 230000004931 aggregating effect Effects 0.000 claims description 2
- 239000000446 fuel Substances 0.000 claims description 2
- 238000013461 design Methods 0.000 abstract description 3
- 230000009286 beneficial effect Effects 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 4
- 238000007726 management method Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 241000512668 Eunectes Species 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 210000000349 chromosome Anatomy 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 239000012634 fragment Substances 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- KRTSDMXIXPKRQR-AATRIKPKSA-N monocrotophos Chemical compound CNC(=O)\C=C(/C)OP(=O)(OC)OC KRTSDMXIXPKRQR-AATRIKPKSA-N 0.000 description 1
- 238000013439 planning Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000013316 zoning Methods 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及共享单车需求预测与投放调度方法，所述方法包括以下步骤：S1、建立XGBoost决策树，将相似和相邻的站点聚集成集群；S2、对单个站点的真实需求进行纠偏，并带入到XGBoost决策树中经训练后得到优化的站点聚类结果；S3、根据训练好的XGBoost决策树预测每个站点的借/还车需求；S4、考虑每个站点容量限制和需求到达分布，计算单车的初始投放量；S5、划分城市的调度分区；S6、建立区域内单车调度模型，获得调度路径的最优选择。本发明从共享单车的需求预测和调度优化着手，考虑多种现实情况和现实问题，设计模型和求解方法，为城市共享单车***的运营决策提供依据，有利于快速实现城市居民的借/还车需求，方便生活，健康出行。
Description
技术领域
本发明涉及共享单车运营的技术领域，具体而言，涉及共享单车需求预测与投放调度方法。
背景技术
目前，共享单车已经融入日常生活，需要配合城市治理，成为公共交通服务的一部分。部分地区政府发布的相关管理措施主要作用表现在：(1)各地政府限制运营商投放共享单车，为了保持用户的满意度，企业只能想办法增加车辆的利用率；(2)政策规定共享单车平台不得收取用户押金，共享单车的金融产品特征被削弱，市场竞争的焦点转移到了精细化管理、智能化出行体验以及打造优质服务上；(3)政府将押扣违规停放的共享单车，为了尽量减少赎金损失，共享单车平台标注可停放区域和禁停区，并且需要及时将违规停放车辆调转至可停区。
因此，现存单车企业的管理模式已从原来通过规模扩张模式转变为了精细化管理。如今的单车市场，只有提升运营效率，提高用户满意度，控制运维成本才能使得企业继续得到发展。做到这些关键是要解决站点单车资源分布不均衡问题和违规停放问题，这需要对共享单车进行科学合理的再分配，通过制定合理的方案，把单车从供应点运到需求点，从禁停点运到可停点，匹配用户需求。事实上，共享单车***中单车利用率不高的问题普遍存在，企业调度均衡仍有较大空间。
发明内容
本发明的主要目的在于提供共享单车需求预测与投放调度方法，可以预测不同时间的城市内用户借/还车需求分布情况，并设计方案决定调运单车的地点、时间、数量和路线，保证单车的供需匹配和***的运营效率，从而解决上述背景技术中提出的技术问题。
为达到上述目的，本发明提供了共享单车需求预测与投放调度方法，包括以下步骤：
S1、建立XGBoost决策树，将相似和相邻的站点聚集成集群；
S2、对单个站点的真实需求进行纠偏，并带入到XGBoost决策树中经训练后得到优化的站点聚类结果；
S3、根据训练好的XGBoost决策树预测每个站点的借/还车需求；
S4、考虑每个站点容量限制和需求到达分布，计算单车的初始投放量；
S5、划分城市的调度分区；
S6、建立区域内单车调度模型，获得调度路径的最优选择。
优选地，步骤S1中所述建立XGBoost决策树，将相似和相邻的站点聚集成集群包括：
分别统计过去d个工作日内每个站点h时段内净车流数y和特征变量loni、lati、temi、raini、speedi，共收集样本N个，该XGBoost决策树模型描述为:
式(1)中，子树的数量用k表示，表示站点i的客户借/还车需求预测值，loni表示站点i的经度，lati表示站点i的纬度，temi表示站点i的温度，raini表示站点i的降水量，speedi表示站点i的风速，该XGBoost决策树模型的目标函数为：
根据泰勒展开理论，式(2)可以转化为：
式(3)中，T为叶子节点的总数，Glf为lf叶子节点下样本的之和，Hlf为lf叶子节点下样本的/>之和；
未知fk(x)＝ωq(x)中，ω表示子树k的叶子节点值，q(x)表示子树k的结构；将式(3)看作是关于变量ωlf.的二次函数，令一阶导数为零，得到叶节点lf的最优解ωlf *和最优目标值Objlf *如下：
确定q(x)即树的结构：从根节点开始，每次判断是否分割节点以及使用哪个特征变量作为分割条件，都取决于增益Gain，Gain的计算如下：
遍历计算各种***方式的增益，选择Gain最大的方式将一个叶节点向下分为左右两个节点，当所有Gain≤0则不再***。
优选地，将待预测日期的站点的经度、纬度、站点地区h时段的温度、降水量、风速输入到XGBoost决策树模型训练所得的K颗子树中，统计相邻的站点j和站点j’的相似度指标其中/>为布尔值，j站点和j’站点在子树k中落入同一个叶子节点则/>否则/>ωk为子树k的权重，0≤Mjj’≤1；
用1-Mjj’表示两个站点之间的距离，进行距离参数为Ma的层级聚类得到站点聚类结果。
优选地，步骤S2中所述对单个站点的真实需求进行纠偏，并带入到XGBoost决策树中经训练后得到优化的站点聚类结果包括：
对于给定时段h和***中的n个站点，设站点i的真实需求为xi,i＝1,2,…,n，***的记录数为fi,i＝1,2,…,n；定义α来表示用户遇到失衡站点时的转移意愿，wij＝α/Dij表示第i站与第j站之间的转移率，其中Dij表示站点i与站点j之间的距离；
S1表示时段h内发生过失衡的站点集，S2表示时段h内从未发生失衡的站点集，|S1|＝n1，|S2|＝n2，n1+n1＝n，S1∪S2＝S；
属于S2的站点，xi与fi的关系如下：
属于S1的站点，xi与fi的关系如下：
在式(6)和(7)中，trans(i,j)表示从站点i转移到站点j的借/还车需求，bool(i,j)表示若站点j是失衡站点i在ti内距离最近的正常工作的站点，bool(i,j)＝1，否则bool(i,j)＝0,站点i∈S1；ti表示站点i处于失衡状态时区；
对于每个属于S1的站点有：
约束条件如下：
fj-xj≥0 j∈S2 (10)
xi≥0, α≥0 i∈S (12)
目标函数如下：
在式(13)中，Sc表示聚类结果中集群C所包含的站点集,|Sc|＝nc。
优选地，步骤S3中根据训练好的XGBoost决策树预测每个站点的借/还车需求包括：
将各站点的记录数量fi替换为xi，得到纠偏后的需求数据，使用纠偏后的需求数据作为标签，对XGBoost决策树模型进行训练，在训练后的模型中通过输入特征变量预测站点i在h时段的真实需求如下：
优选地，步骤S4中所述考虑每个站点容量限制和需求到达分布，计算单车的初始投放量包括：
某站点单个时段内被拒需求数量的估计：
给定站点j和时段h＝[t0 h,ts h]，按照泊松分布分别由h时段的借车需求量概率分布和还车需求量概率分布产生随机值和/>则h内单车净变化量/>设站点容量为C，在初始时刻t0 h站点单车数量为q的条件下，h时段内产生的被拒借车需求数量的期望值/>表示为：
在式(15)中，u1表示被拒借车需求数量；
h时段内产生的被拒还车需求数量的期望值表示为：
在式(16)中，u2表示被拒还车需求数量；
在初始单车投放量为q的条件下，h时段内被拒需求数量的期望值为：
分别计算q从0到C的取/>最小值对应的q为h时段站点j的最优初始投放量。
优选地，某站点连续时段内被拒需求数量的估计：
若测算连续n个时段h1,..hm,..hn内站点j产生被拒需求数量的期望值，hm的初始库存是hm-1的末点库存，末点库存有三种情况：
h时段站点的单车库存量ph由q变为0的概率：
h时段站点的单车库存量ph由q变为C的概率：
h时段站点的单车库存量ph由q变为q’的概率：
初始单车投放量为q的条件下，站点在连续n个时段内产生的被拒需求数量的期望值为E(δf)：
取E(δf)最小值对应的q为该站点连续n个时段的最优初始投放量。
优选地，步骤S5中所述划分城市的调度分区包括：
S51、首先筛选出净借需求量绝对值过大的站点，即极端需求站点，依次以这些站点为中心，搜索其附近1公里以内的全部站点，分别计算这些站点与中心站点间的距离，从最近的站点开始收纳入“零需”站点集，当站点集的净借需求量总和已经为零或者略大于零，则终止搜寻，为下一个极端需求站点构建“零需”站点集；
S52、将周围站点全部吸纳仍未达到终止条件，则以当前站点集中心位置为虚拟站点，以集内站点的净借需求量总和为虚拟站点的净借需求量，重复上述过程，吸纳新的站点加入；
S53、若站点集容量达到阈值仍然未达成终止条件，则释放这个站点集，直接为下一个极端需求站点构建“零需”站点集；若某站点附近1公里以内并没有还未加入站点集的站点，则释放这个站点集，将这个站点加入“落单”站点集，直接进行下一个“零需”站点集的构建；
S54、完成极端需求站点的站点集之后，在剩下的站点中找出地理位置上的***站点，依次为其构建“零需”站点集，循环这一过程，直到未加入站点集的站点数量与落单站点集数量相同；
S55、将落单站点视为一个“零需”站点集，计算所有“零需”集站点集的地理中心位置，用一个虚拟站点代替一个站点集，以经度最高的虚拟站点为始，与周围的虚拟站点合并，当实际站点数超过40且区域内所有站点的单车库存之和处于该区域最优单车投放量的上界和下界之间时，停止合并，作为一个调度分区；
S56、重复S51-S54，完成城市内所有站点的调度分区划分。
优选地，步骤S6中所述建立区域内单车调度模型，获得调度路径的最优选择包括：
规定区域内有若干共享单车站点，用集合S1＝{1，2，……，N}表示，由标准运输货车进行调度作业，运输货车的出发点用S2＝{0}表示。每辆货车的载量上限为Q。在整个调度作业过程中，每辆作业货车的实时载量e必须满足0≤e≤Q；
目标函数1：调运成本和碳排放的最小化：
目标函数2：用户不满惩罚最小化：
约束条件：
Lij＝Lki+qi,k∈{0,…,N-1},Rki＝1,i∈S1,j∈{2,…,N}∪{0}(27)
Lij≤M·Rij(28)
L0j ＝Lj0 ＝ 0 (30)
qi为整数(36)
S表示调度货车的调度路径，S＝S1∪S2＝{0,1,2,…,N}，其中0表示调度货车的出发点，其他表示各站点；qi表示在站点投放的单车的数量，qi<0时表示在站点i取走-qi辆单车；Ri,j为0-1变量，如果从站点i到站点集j在货车的调度路线中，Ri,j＝1，否则Ri,j＝0；gi表示站点i在调度工作进行前的单车数量；fi表示站点i在调度工作完成后的单车数量；F(x,i)表示表示站点i的单车数量为x时，预期会产生的损失需求的数量；Lij表示调度货车从站点i行驶到站点j过程中所载单车数量；Q表示调度货车的的容量；Dij表示站点i到站点j之间的运输距离；Ce表示调度货车单位公里单位千克耗用的燃油量；Cd表示每产生一个被拒绝的需求的不满惩罚；Ci表示站点i的容量；M表示一个无限大的正值；B表示消除子回路的辅助变量。
优选地，通过混合多目标粒子群优化算法获得调度路径的最优选择，包括：
S61、生成初始解集：计算区域内站点间的距离矩阵，随机选择一个站点作为路线首个站点，然后在距离该站点最近的4个站点中以等概率随机选取任意一个作为下一个目的地，然后重复这一过程，避免选择已被选择的站点，直到将所有站点加入路线；将生成的路线作为解Rij输入式(22)-(35)，调用Python的scipy.optimize求解得到各个站点处的最优装/卸单车数量qi，与Rij构成一个完整的初始解；以相同的方式形成其他初始解，构建初始解集；
S62、适应度函数设计及个体最优选择：通过比较粒子的适应度值来获得粒子的优势关系，适应度函数来自目标函数1和目标函数2，指定如下：
S63、构建非劣解集：定义帕累托支配关系:对于粒子z1和z2，如f1(z1)≤f1(z2)且f2(z1)<f2(z2)，或者f1(z1)<f1(z2)且f2(z1)≤f2(z2)，则粒子z1支配粒子z2；同样，如果f1(z2)≤f1(z1)且f2(z2)<f2(z1)，或者f1(z2)<f1(z1)且f2(z2)≤f2(z1)，则粒子z2支配粒子z1；否则，粒子z1和z2没有支配关系；
对于任意粒子z*，当且仅当f1(z*)≤f1(z)且f2(z*)<f2(z)或者f1(z*)<f1(z)且f2(z*)≤f2(z)时，粒子z*被称为非支配粒子或帕累托粒子，所有这些非支配粒子都记录在非劣解集REP中；
S64、选择全局最优选择：全局最佳粒子gbest是从REP中选择的，REP保存搜索过程中找到的所有非支配解；在搜索开始时，将在初始化阶段构造的非支配粒子全部添加到REP；在搜索过程中，将每次迭代中发现的当前非支配粒子与REP中的所有解一一比较；如果一个当前粒子被REP中的粒子支配，则将其丢弃；否则，可以将此类粒子附加到REP；若REP中有粒子被新成员支配，则将这些粒子从REP中删除；当REP超过其最大容量时，用自适应网格方法删除多余的解，保持解决方案在REP中分布良好；对于gbest的选择，将解决方案空间划分为多个相等的网格，其中具有较少粒子的网格具有较高的被选择机会；设调整参数为β，网格数为n，网格j包含的粒子数量由nj表示，则第i个网格中的解被选择为gbest的概率为：
S65、交叉算子设计：使用遗传算法的交叉算子来更新算法中的粒子：首先由pbest更新粒子，然后由gbest更新粒子，粒子交叉算子其中Xi(t+1)和Xi(t)分别表示第i个粒子在当前迭代(t)和下一次迭代(t+1)中的位置；pbest(t)是第i个粒子的历史最优位置，gbest(t)是总体的全局最优位置；/>是交叉运算符；
S66、变异算子设计：使用结合VNS的变异算子对粒子进行更新，采用三种邻域结构：
(1)节点***：随机选择一个站点并将其从其原始位置移除，然后将其***另一个随机位置；
(2)节点交换：随机选择两个站点，并交换其位置；
(3)随机选择两个站点，将它们之间的顺序以相反的顺序排列；
以上三个邻域结构随机执行，进行变异操作后，如果新粒子支配了原始粒子，将替换原始粒子；如果原始粒子占主导地位，则丢弃新粒子；如果这两个粒子之间没有支配关系，则随机选择其中之一。
综上，本发明提供的共享单车需求预测与投放调度方法从共享单车的需求预测和调度优化着手，考虑多种现实情况和现实问题，设计模型和求解方法，为城市共享单车***的运营决策提供依据，有利于快速实现城市居民的借/还车需求，方便生活，健康出行；提高了共享单车平台的服务水平，提升了企业的品牌口碑和客户粘性，并且降低了控制成本，提高了企业竞争力；同时有利于减少城市碳排放，保护环境，建设智慧城市。
附图说明
图1为本发明的实施方式提供的XGBoost决策树对站点聚类的示意图；
图2是本发明的实施方式提供的共享单车需求预测与投放调度方法的预测模型的全流程图；
图3为预测模型中不同参数max_dept设置下的拟合度变化情况；
图4为预测模型中不同参数n_estimator设置下的拟合误差变化情况；
图5为以cover为指标的各个特征在预测模型中的重要程度；
图6为2022年1月24日8:00-9:00共享单车***借车需求的站点聚类结果；
图7为第一个集群站点的地理图像；
图8为2022年1月25日的h时段内净需求量预测结果分布；
图9为估算连续n个时段被拒需求数量期望值的流程图；
图10为调度分区示意图；
图11为交叉算子示例；
图12为三种邻域结构示意图；
图13为根据调度模型在最佳效益下获得的调度路径。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的实施方式提供了一种共享单车需求预测与投放调度方法，包括以下步骤：
S1、建立XGBoost决策树，将相似和相邻的站点聚集成集群；
S2、对单个站点的真实需求进行纠偏，并带入到XGBoost决策树中经训练后得到优化的站点聚类结果；
S3、根据训练好的XGBoost决策树预测每个站点的借/还车需求；
S4、考虑每个站点容量限制和需求到达分布，计算单车的初始投放量；
S5、划分城市的调度分区；
S6、建立区域内单车调度模型，获得调度路径的最优选择。
本发明首先通过XGBoost将相似和相邻的站点聚集成集群。在真实需求的估计中，站点与站点之间的相互作用是不可忽视的。当一个车站没有可用的车桩/单车，而顾客在附近的站点借了一辆，这个满的/空的车站会产生需求外溢到邻近的车站。因此，在不考虑溢出效应的情况下，单独处理***中的站点，会导致对需求的低估和高估。这意味着我们甚至不能将始终运行正常(没有发生失衡)的站点的借/还车记录视为真正的需求，因为它们的部分记录可能来自相邻的满/空站。此外，对车站的需求不仅受到自然因素(如天气条件、人口、土地利用等)的影响，还受到周围车站分布的影响。例如，对一个车站的需求将会由于附近新车站的建立而减少。因此，关注一个区域的需求量而不是单个车站有助于避免站点间转移和交互的影响，更精确地捕捉当地的出行模式，这就需要站点聚类。
若固定一天内的某个时段h(一个小时)，站点间需求的差异可以由地理位置和天气条件来解释。图1为本发明的实施方式提供的XGBoost决策树对站点聚类的示意图，如图1所示，本发明步骤S1中所述建立XGBoost决策树，将相似和相邻的站点聚集成集群包括：
分别统计过去d个工作日内每个站点h时段内净车流数y和特征变量loni、lati、temi、raini、speedi，共收集样本N个，该XGBoost决策树模型描述为:
式(1)中，子树的数量用k表示，表示站点i的客户借/还车需求预测值，loni表示站点i的经度，lati表示站点i的纬度，temi表示站点i的温度，raini表示站点i的降水量，speedi表示站点i的风速，该XGBoost决策树模型的目标函数为：
根据泰勒展开理论，式(2)可以转化为：
式(3)中，T为叶子节点的总数，Glf为lf叶子节点下样本的之和，Hlf为lf叶子节点下样本的/>之和；
未知fk(x)＝ωq(x)中，ω表示子树k的叶子节点值，q(x)表示子树k的结构；将式(3)看作是关于变量ωlf.的二次函数，令一阶导数为零，得到叶节点lf的最优解ωlf *和最优目标值Objlf *如下：
确定q(x)即树的结构：从根节点开始，每次判断是否分割节点以及使用哪个特征变量作为分割条件，都取决于增益Gain，Gain的计算如下：
遍历计算各种***方式的增益，选择Gain最大的方式将一个叶节点向下分为左右两个节点，当所有Gain≤0则不再***。
进一步地，将待预测日期的站点的经度、纬度、站点地区h时段的温度、降水量、风速输入到XGBoost决策树模型训练所得的K颗子树中，统计相邻的站点j和站点j’的相似度指标其中/>为布尔值，j站点和j’站点在子树k中落入同一个叶子节点则/>否则/>ωk为子树k的权重，0≤Mjj’≤1；
用1-Mjj’表示两个站点之间的距离，进行距离参数为Ma的层级聚类得到站点聚类结果。
需要注意的是训练子树的最大深度max_dept和子树的个数n_estimator是需要谨慎确定的参数，过小影响模型精度，过大容易发生过拟合。
上述XGBoost决策树模型训练以位置和天气信息为输入，借/还车记录数量为目标，训练XGBoost模型。上述过程是本发明提供的预测模型的第一阶段，即提出模型结构，根据训练好的模型结构对站点进行聚类，然而这些作为标签的值并不是真正的需求。在站点失衡(空站或者满站)的情况下，用户无法借车或者还车，此时***的交易记录与用户初始的借/还车需求不符，即一部分借/还车需求损失了。因此，需要估计每个站点用户的原始需求，对***记录数据进行纠偏，要使结果接近真正的用户需求。
纠偏模型中设置的变量设置如下:
在本发明的一种实施方式中，步骤S2中所述对单个站点的真实需求进行纠偏，并带入到XGBoost决策树中经训练后得到优化的站点聚类结果包括：
对于给定时段h和***中的n个站点，设站点i的真实需求为xi,i＝1,2,…,n，***的记录数为fi,i＝1,2,…,n；定义α来表示用户遇到失衡站点时的转移意愿，wij＝α/Dij表示第i站与第j站之间的转移率，其中Dij表示站点i与站点j之间的距离；
属于S2的站点，xi与fi的关系如下：
属于S1的站点，xi与fi的关系如下：
对于每个属于S1的站点有：
约束条件如下：
fj-xj≥0 j∈S2 (10)
xi≥0, α≥0 i∈S (12)
目标函数如下：
上述过程是本发明提供的预测模型的第二阶段，即对于真实借车需求的估计，由式(13)、(6)，(7)和(9)-(12)建立非线性规划模型，该模型基于一个假设：当遭遇站点失衡时，用户将转移到在ti期间始终保持运行的车站。虽然可能会有一个更近的车站，在ti的部分时间内还没有失衡，但它不会被选择，因为它得库存已经处于危险水平，用户担心在步行期间它变成空站。
对于还车需求的估计，与借车需求的模型略有不同。当用户遇到一个满站时，他们必须归还自行车以避免更多的付款，这意味着w＝1。因此，将借车需求纠偏的模型删去式(9)，并且忽略α/Dij，就可以用于还车需求量的纠偏工作。
图2是本发明的实施方式提供的共享单车需求预测与投放调度方法的预测模型的全流程图。如图3所示，本发明步骤S3中根据训练好的XGBoost决策树预测每个站点的借/还车需求包括：
将各站点的记录数量fi替换为xi，得到纠偏后的需求数据，使用纠偏后的需求数据作为标签，对XGBoost决策树模型进行训练，在训练后的模型中通过输入特征变量预测站点i在h时段的真实需求如下：
上述过程是本发明提供的预测模型的第三阶段，即通过纠偏后的需求数据对XGBoost决策树模型进行训练，然后利用训练后的模型预测单个站点的真实需求。
以下为XGBoost决策树模型参数调试和聚类结果示例：
本示例中数据处理与分析工作使用Python(Anaconda 3.7)编码和运行，在IntelCore i7-7500 [email protected]上计算实验，其中Xgboost模型由TensorFlow CPU(1.14.0)训练，设置模型的学习率eta＝0.2，复杂度惩罚(γ＝0.5，λ＝1)，尝试不同的max_dept和n_estimator调整参数。本文研究样本选自2021年1月1日到2022年1月25日的所有工作日，以20个工作日为窗口，采取滚动预测的方式。以预测2022年1月25日的借/还车需求为例，模型的训练数据应取自2021年12月17日到2022年1月24日。
图3为预测模型中不同参数max_dept设置下的拟合度变化情况。树的最大深度max_dept越大，模型的拟合度越高，同时也越容易发生过拟合。以样本集训练xgboost模型时，分别设置不同max_dept，其对应的拟合度如图3所示，拟合度在深度6之后增长缓慢，因此设置参数max_dept＝6。
图4为预测模型中不同参数n_estimator设置下的拟合误差变化情况。子树的个数n_estimator越大，模型的拟合误差越小，但也同时越容易发生过拟合。以样本集训练xgboost模型时，分别设置不同的学习器个数，其对应的均方根误差(rmse)如图4所示，rmse在子树个数为8之后增长缓慢，因此设置参数n_estimator＝8。
前8棵子树对应的rmse分别为5.9573，5.6791，5.4581，5.3188，5.1856，5.1139，5.0279和4.9403。按四分位法统计标签数据的Q1和Q3，得四分位距IQR＝Q3-Q1＝9。这里认为在没有任何预测基础上随意取值，和目标值的平均差距水平为8。那么第1个子树的贡献为8-5.9573＝3.0427，其他的贡献依次为0.2782，0.2210，0.1393，0.1332，0.0717，0.086，0.0876。进而计算前八棵子树的权重分别为74.95％，6.85％，5.44％，3.43％，3.28％，1.77％，2.12％和2.16％，用于层级聚类。
训练好的XGboost模型可以由下表表示，f0～f4分别是纬度、经度、温度、风速、降水5个特征。表中信息可知，所有8个子树的总节点数量为949，其中第八棵子树有120个树节点(根节点、***节点和叶子节点的总和)。第一行的含义为：第一个树的根节点，其***条件为“纬度<40.767029”,若满足条件，进入第一棵子树的第二个节点(ID为0-1的节点)，否则进入第一棵子树的第三个节点(ID为0-2的节点)，在这个节点以这一条件进行样本划分的增益Gain为5425.8652(Gain的计算见3.1.2节式子3-5)。由于叶子节点已经不再***，因此***条件部分都是空值NaN。
表1训练好的xgboost模型的部分展示
| Tree | Node | ID | Feature | Split | Yes | No | Gain | |
| 0 | 0 | 0 | 0-0 | f0 | 40.767029 | 0-1 | 0-2 | 5425.8652 |
| 1 | 0 | 1 | 0-1 | f0 | 40.736374 | 0-3 | 0-4 | 14389.0361 |
| 2 | 0 | 2 | 0-2 | f1 | -73.96813 | 0-5 | 0-6 | 10376.0088 |
| 3 | 0 | 3 | 0-3 | f0 | 40.731632 | 0-7 | 0-8 | 3230.9951 |
| 4 | 0 | 4 | 0-4 | f1 | -73.98779 | 0-9 | 0-10 | 11749.0039 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 945 | 7 | 116 | 7-116 | Leaf | NaN | NaN | NaN | 0.3539 |
| 946 | 7 | 117 | 7-117 | Leaf | NaN | NaN | NaN | -1.4300 |
| 947 | 7 | 118 | 7-118 | Leaf | NaN | NaN | NaN | -0.3808 |
| 948 | 7 | 119 | 7-119 | Leaf | NaN | NaN | NaN | 0.1083 |
| 949 | 7 | 120 | 7-120 | Leaf | NaN | NaN | NaN | -0.0112 |
在xgboost模块中，衡量特征的重要性程度有3个指标。Weight表示在所有树中，某特征被用来***节点的次数。cover表示在所有树中，某特征在***节点时处理(覆盖)的样本的平均数量。Gain表示某特征在每次***节点时带来的平均增益。图5为以cover为指标的各个特征在预测模型中的重要程度，如图5所示，可见位置信息是显著重要的，因为它包含了土地利用、人口密度、建筑设施等多个因素，第二重要的是降雨因素，然后是温度，最后是风速。
图6为2022年1月24日8:00-9:00共享单车***借车需求的站点聚类结果。由训练好的xgboost模型对站点聚类，以2022年1月24日上午8:00-9:00为例，当天的活跃站点由997个，在DataSet中取出当天的样本构成测试集输入模型，得到相似度矩阵，以距离阈值参数Ma＝0.05进行层级聚类，获得在2022-01-24时段h8的天气条件下的站点聚类结果，如图6所示。
以下为需求量的纠偏和预测示例：
由于需求量的纠偏以整个***为一个规划问题，以第一个集群站点为例，展示纠偏工作得过程分析。图7为第一个集群站点的地理图像，如图7所示，第一个集群站点包括ID为7713.11,7727.08,7692.11,7713.01,7680.03,7650.05,7646.04,和7634.01的8个站点。在本节编码为站点1～8。
在2022-01-24的时段[8:00,9:00]内，各个站点的单车库存变化情况展示于表2。
表2第一集群中个站点在h＝[8:00,9:00]内的库存跟踪
由表2可知，站点1在8:00-8:53之间处于失衡状态(空站)停止服务，故t1＝[8:00,8:53]，同理，t2＝[8:00,9:00]，t5＝[8:37,8:48]，t6＝[8:52,9:00]，t7＝[8:48,9:00]。站点3、4、8在h内始终保持可用。然会根据3.2.2节，描述八个站点xi和fi之间的关系如下：
f3-x3＝trans(2,3)
f4-x4＝trans(1,4)
f8-x8＝trans(6,8)
/>
根据***历史记录，给定(f1＝0,f2＝0,f3＝2,f4＝3,f5＝3,f6＝4,f7＝4,f8＝1)和(D1,4＝241.02,D2,3＝225.22,D5,6＝247.89,D6,8＝231.39,D7,5＝474.86)。使用SciPy.Optimize模块得到解(α＝93.4857,x1＝3,x2＝2,x3＝1,x4＝2,x5＝3,x6＝4,x7＝4x8＝1)。这表明在t1内，1名用户从站点1转移到站点4，2名用户在站点1放弃了租借共享单车的想法；在t1内，1名用户从站点2转移到3站，1名顾客在站点2放弃租车；在t7内，没有顾客到达站点7；在t5内，没有顾客到达站点5；在t6内，没有顾客到达站点6。
由于本例中没有站点发生全满的失衡状态，因此可直接将***记录中的还车记录数量视为各站点的真实还车需求。如果有满站的情况发生，则用同样的方法(忽略α/Dij)对各站点的还车记录进行校正。在整个***中，2022年1月24日当天有41个站点发生满站失衡，67个站点发生空站失衡，两种站点在失衡期间涉及到的溢出站点有86个。受到溢出效应影响的站点数小于发生失衡的站点数，是因为可能出现两个失衡站点对同一个站点产生溢出效应，另外空的失衡站点和满的失衡站点之间可能出现需求的相互溢出效应。建立借车需求的纠偏优化模型时，待求量有68个(67个站点的真实需求xi和转移系数α)，模型求解用时16秒；建立还车需求的纠偏优化模型时，待求量有41个(41个站点的真实需求xi)，模型求解用时7秒。
图8为2022年1月25日的h时段内净需求量预测结果分布。分别完成20各工作日在h＝[8:00,9:00]内的数据矫正，用新数据xi更新原数据集DataSet的标签，再次训练XGBoost模型，得到2022年1月25日08:00-09:00之间的客户需求预测结果。由各站点的借车需求预测值减去还车需求预测值，得到净借车需求量，如图8所示。
由图8可以发现呈现净借需求量分布明显左偏。因为在工作日的上午8:00～9:00，有个别地方如写字楼附近会有集中的、大量的还车需求，而借车需求可能分散在各个居民区，没有非常明显的借车需求。
用平均绝对百分比误差devi衡量***中所有站点的预测值和记录值的平均偏差水平，计算得2.8836，计算公式如下：
在本发明的一种实施方式中，进一步考虑每个站点容量限制和需求到达分布，讨论单车的初始投放量方案。假设单位时间内需求的到达服从泊松分布，使用前文所得站点各时段的借/还车需求量预测值，可以得到单位时间内到达的借/还车需求量的一组随机值。多次模拟借/还车需求到达站点的情况，跟踪站点的单车和车桩数量，然后统计固定时段内发生的被拒借车需求数量δw和被拒还车需求数量δr，由二者多次模拟结果均值的合表示被拒需求数量的期望δf＝δw+δr，度量站点产生的用户不满意。
具体地，步骤S4中所述考虑每个站点容量限制和需求到达分布，计算单车的初始投放量包括：
一方面，某站点单个时段内被拒需求数量的估计：
给定站点j和时段按照泊松分布分别由h时段的借车需求量概率分布和还车需求量概率分布产生随机值/>和/>则h内单车净变化量/>设站点容量为C，在初始时刻t0 h站点单车数量为q的条件下，h时段内产生的被拒借车需求数量的期望值/>表示为：
在式(15)中，u1表示被拒借车需求数量；
h时段内产生的被拒还车需求数量的期望值表示为：
在式(16)中，u2表示被拒还车需求数量；
在初始单车投放量为q的条件下，h时段内被拒需求数量的期望值为：
分别计算q从0到C的取/>最小值对应的q为h时段站点j的最优初始投放量。
另一方面，某站点连续时段内被拒需求数量的估计：
若测算连续n个时段h1,..hm,..hn内站点j产生被拒需求数量的期望值，hm的初始库存是hm-1的末点库存，末点库存有三种情况：
h时段站点的单车库存量ph由q变为0的概率：
h时段站点的单车库存量ph由q变为C的概率：
h时段站点的单车库存量ph由q变为q’的概率：
图9为估算连续n个时段被拒需求数量期望值的流程图。hm的初始单车量可能是0～C的任意值，并且有对应概率，对于每个初始单车量的情况，hm时段末点也都对应有0～C库存量的离散概率分布(4.1.1节的过程)。hm时段末点库存为各个值的可能性由全概率公式计算，得到的概率分布将作为hm+1时段的初始条件。迭代过程如图9所示。
初始单车投放量为q的条件下，站点在连续n个时段内产生的被拒需求数量的期望值为E(δf)：
取E(δf)最小值对应的q为该站点连续n个时段的最优初始投放量。
共享单车主要解决“最后”一公里问题，大部分人的骑行距离在500-1000m之间，90％的用户在两公里之内骑行。因此单车在城市间的流动和大概率在小范围内进行，一个区域中的单车数量在较长一段时间内可以保持粗略平衡。基于这一点，本发明提出城市调度分区的概念。理想情况下，一个共享单车***的总的供需关系是平衡的，单车“有借必有还”，***中有单车紧缺站点也必定有单车过多站点。调度分区的思路是在获得各个站点的净借需求量预测值后，将城市在区域上分为多块，尽量保证每个区块内站点的净借需求量之和为零，即每个区块包括的站点集内部可以自给自足。理论上，为每个人调度分区配备1辆调度货车，该车在作业初始时无需携带单车，通过完成分区内站点间的相互调度实现再平衡，作业完成后也不会携带多余单车。
具体地，在本发明的一种实施方式中，步骤S5中所述划分城市的调度分区包括：
S51、首先筛选出净借需求量绝对值过大的站点，即极端需求站点，依次以这些站点为中心，搜索其附近1公里以内的全部站点，分别计算这些站点与中心站点间的距离，从最近的站点开始收纳入“零需”站点集，当站点集的净借需求量总和已经为零或者略大于零，则终止搜寻，为下一个极端需求站点构建“零需”站点集；
S52、将周围站点全部吸纳仍未达到终止条件，则以当前站点集中心位置为虚拟站点，以集内站点的净借需求量总和为虚拟站点的净借需求量，重复上述过程，吸纳新的站点加入；
S53、若站点集容量达到阈值仍然未达成终止条件，则释放这个站点集，直接为下一个极端需求站点构建“零需”站点集；若某站点附近1公里以内并没有还未加入站点集的站点，则释放这个站点集，将这个站点加入“落单”站点集，直接进行下一个“零需”站点集的构建；
S54、完成极端需求站点的站点集之后，在剩下的站点中找出地理位置上的***站点，依次为其构建“零需”站点集，循环这一过程，直到未加入站点集的站点数量与落单站点集数量相同；
S55、将落单站点视为一个“零需”站点集，计算所有“零需”集站点集的地理中心位置，用一个虚拟站点代替一个站点集，以经度最高的虚拟站点为始，与周围的虚拟站点合并，当实际站点数超过40且区域内所有站点的单车库存之和处于该区域最优单车投放量的上界和下界之间时，停止合并，作为一个调度分区；
S56、重复S51-S54，完成城市内所有站点的调度分区划分。
图10为根据上述方法实施的调度分区示意图。
为每个分区分别配备一辆调度货车，则各个分区的调度工作都是个独立的问题，在本发明的实施方式中，重点描述一个分区的调度问题。
调度模型中设置的变量设置如下:
步骤S6中所述建立区域内单车调度模型，获得调度路径的最优选择包括：
规定区域内有若干共享单车站点，用集合S1＝{1，2，……，N}表示，由标准运输货车进行调度作业，运输货车的出发点用S2＝{0}表示。每辆货车的载量上限为Q。在整个调度作业过程中，每辆作业货车的实时载量e必须满足0≤e≤Q；
目标函数1：调运成本和碳排放的最小化：
目标函数2：用户不满惩罚最小化：
约束条件：
/>
Lij＝Lki+qi,k∈{0,…,N-1},Rki＝1,i∈S1,j∈{2,…,N}∪{0}(27)
Lij≤M·Rij(28)
L0j ＝Lj0 ＝ 0 (30)
约束条件(24)到(26)限制了货车在站点的取车量上限为站点的初始库存，放车量的上限不能超过该站点剩余容量。约束条件(27)到(29)保证运输车辆所载单车数量为非负且不超过货车最大容量限制。约束条件(30)约束运输车辆空载出发并且空载返回，保证该区域内的单车总数量载调度前后一样。约束条件(31)确保表示所有的站点最多被访问一次，以满足需求。约束条件(32)确保路线在站点间前后连接形成。约束条件(33)确保重平衡车辆从某地出发并且返回该地。约束条件(34)是消除子回路约束。约束条件(35)和(36)定义了变量的定义域。
在本发明的一种实施方式中，通过混合多目标粒子群优化算法获得调度路径的最优选择，包括：
S61、生成初始解集：计算区域内站点间的距离矩阵，随机选择一个站点作为路线首个站点，然后在距离该站点最近的4个站点中以等概率随机选取任意一个作为下一个目的地，然后重复这一过程，避免选择已被选择的站点，直到将所有站点加入路线；将生成的路线作为解Rij输入式(22)-(35)，调用Python的scipy.optimize求解得到各个站点处的最优装/卸单车数量qi，与Rij构成一个完整的初始解；以相同的方式形成其他初始解，构建初始解集；
S62、适应度函数设计及个体最优选择：通过比较粒子的适应度值来获得粒子的优势关系，适应度函数来自目标函数1和目标函数2，指定如下：
S63、构建非劣解集：定义帕累托支配关系:对于粒子z1和z2，如f1(z1)≤f1(z2)且f2(z1)<f2(z2)，或者f1(z1)<f1(z2)且f2(z1)≤f2(z2)，则粒子z1支配粒子z2；同样，如果f1(z2)≤f1(z1)且f2(z2)<f2(z1)，或者f1(z2)<f1(z1)且f2(z2)≤f2(z1)，则粒子z2支配粒子z1；否则，粒子z1和z2没有支配关系；
对于任意粒子z*，当且仅当f1(z*)≤f1(z)且f2(z*)<f2(z)或者f1(z*)<f1(z)且f2(z*)≤f2(z)时，粒子z*被称为非支配粒子或帕累托粒子，所有这些非支配粒子都记录在非劣解集REP中；
S64、选择全局最优选择：全局最佳粒子gbest是从REP中选择的，REP保存搜索过程中找到的所有非支配解；在搜索开始时，将在初始化阶段构造的非支配粒子全部添加到REP；在搜索过程中，将每次迭代中发现的当前非支配粒子与REP中的所有解一一比较；如果一个当前粒子被REP中的粒子支配，则将其丢弃；否则，可以将此类粒子附加到REP；若REP中有粒子被新成员支配，则将这些粒子从REP中删除；当REP超过其最大容量时，用自适应网格方法删除多余的解，保持解决方案在REP中分布良好；对于gbest的选择，将解决方案空间划分为多个相等的网格，其中具有较少粒子的网格具有较高的被选择机会；设调整参数为β，网格数为n，网格j包含的粒子数量由nj表示，则第i个网格中的解被选择为gbest的概率为：
S65、交叉算子设计：使用遗传算法的交叉算子来更新算法中的粒子：首先由pbest更新粒子，然后由gbest更新粒子，粒子交叉算子其中Xi(t+1)和Xi(t)分别表示第i个粒子在当前迭代(t)和下一次迭代(t+1)中的位置；pbest(t)是第i个粒子的历史最优位置，gbest(t)是总体的全局最优位置；/>是交叉运算符；
图11为交叉算子示例。如图11中(a)所示，首先选择Xi(t)及其pbesti作为父代P1和P2，然后通过交换它们随机选择的染色体片段来更新车辆路线，将路线中重复的站点用就近原则替换为遗失的站点，然后调整各站点装/卸单车数量；在交叉操作之后，生成两个后代粒子O1和O2，如图11中(b)所示；将它们通过由式(37)计算的适应度值进行比较，选择不被支配的一个以执行后续操作；如果它们没有支配关系，则随机选择其中之一；图11中，后代O2被选作新的父代P1'，gbest是另一个父代P2'；以相同的方式处理P1'和P2'产生后代粒子O1'和O2'，如图11(c)所示；同样，将非支配的个体作为更新的粒子Xi(t+1)。
S66、变异算子设计：使用结合VNS的变异算子对粒子进行更新，采用三种邻域结构，图12为三种邻域结构示意图，如图12所示：
(1)节点***：随机选择一个站点并将其从其原始位置移除，然后将其***另一个随机位置；
(2)节点交换：随机选择两个站点，并交换其位置；
(3)随机选择两个站点，将它们之间的顺序以相反的顺序排列；
以上三个邻域结构随机执行，进行变异操作后，如果新粒子支配了原始粒子，将替换原始粒子；如果原始粒子占主导地位，则丢弃新粒子；如果这两个粒子之间没有支配关系，则随机选择其中之一。
图13为根据调度模型在最佳效益下获得的调度路径。由于分区后每个区块的站点数量仍然较多，为了方便展示，本发明取了四个“零需”站点集为例，分别进行调度求解。取算法运行10次后的最优结果作为每个算例的求解结果，各个算例的最佳效益路径如图13所示。
四个算例的站点数量分别15，20，25和30，对应的求解运行时长分别为18.37秒，27.7秒，42.38秒和59.38秒。
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。
Claims (5)
1.一种共享单车需求预测与投放调度方法，包括以下步骤：
S1、建立XGBoost决策树，将相似和相邻的站点聚集成集群；
S2、对单个站点的真实需求进行纠偏，并带入到XGBoost决策树中经训练后得到优化的站点聚类结果；
S3、根据训练好的XGBoost决策树预测每个站点的借/还车需求；
S4、考虑每个站点容量限制和需求到达分布，计算单车的初始投放量；
S5、划分城市的调度分区；
S6、建立区域内单车调度模型，获得调度路径的最优选择；
步骤S1中所述建立XGBoost决策树，将相似和相邻的站点聚集成集群包括：
分别统计过去d个工作日内每个站点h时段内净车流数y和特征变量loni、lati、temi、raini、speedi，共收集样本N个，该XGBoost决策树模型描述为:
式(1)中，子树的数量用k表示，表示站点i的客户借/还车需求预测值，loni表示站点i的经度，lati表示站点i的纬度，temi表示站点i的温度，raini表示站点i的降水量，speedi表示站点i的风速，该XGBoost决策树模型的目标函数为：
根据泰勒展开理论，式(2)可以转化为：
式(3)中，T为叶子节点的总数，Glf为lf叶子节点下样本的'之和，Hlf为lf叶子节点下样本的/>”之和；
未知fk(x)＝ωq(x)中，ω表示子树k的叶子节点值，q(x)表示子树k的结构；将式(3)看作是关于变量ωlf.的二次函数，令一阶导数为零，得到叶节点lf的最优解ωlf *和最优目标值Objlf *如下：
确定q(x)即树的结构：从根节点开始，每次判断是否分割节点以及使用哪个特征变量作为分割条件，都取决于增益Gain，Gain的计算如下：
遍历计算各种***方式的增益，选择Gain最大的方式将一个叶节点向下分为左右两个节点，当所有Gain≤0则不再***；
将待预测日期的站点的经度、纬度、站点地区h时段的温度、降水量、风速输入到XGBoost决策树模型训练所得的K颗子树中，统计相邻的站点j和站点j’的相似度指标其中/>为布尔值，j站点和j’站点在子树k中落入同一个叶子节点则/>否则/>ωk为子树k的权重，0≤Mjj’≤1；
用1-Mjj’表示两个站点之间的距离，进行距离参数为Ma的层级聚类得到站点聚类结果；
步骤S2中所述对单个站点的真实需求进行纠偏，并带入到XGBoost决策树中经训练后得到优化的站点聚类结果包括：
对于给定时段h和***中的n个站点，设站点i的真实需求为xi,i＝1,2,…,n，***的记录数为fi,i＝1,2,…,n；定义α来表示用户遇到失衡站点时的转移意愿，wij＝α/Dij表示第i站与第j站之间的转移率，其中Dij表示站点i与站点j之间的距离；
S1表示时段h内发生过失衡的站点集，S2表示时段h内从未发生失衡的站点集，|S1|＝n1，|S2|＝n2，n1+n1＝n，S1∪S2＝S；
属于S2的站点，xi与fi的关系如下：
属于S1的站点，xi与fi的关系如下：
在式(6)和(7)中，trans(i,j)表示从站点i转移到站点j的借/还车需求，bool(i,j)表示若站点j是失衡站点i在ti内距离最近的正常工作的站点，bool(i,j)＝1，否则bool(i,j)＝0,站点i∈S1；ti表示站点i处于失衡状态时区；
对于每个属于S1的站点有：
约束条件如下：
fj-xj≥0 j∈S2 (10)
xi≥0,α≥0 i∈S (12)
目标函数如下：
在式(13)中，Sc表示聚类结果中集群C所包含的站点集，|Sc|＝nc；
步骤S3中根据训练好的XGBoost决策树预测每个站点的借/还车需求包括：
将各站点的记录数量fi替换为xi，得到纠偏后的需求数据，使用纠偏后的需求数据作为标签，对XGBoost决策树模型进行训练，在训练后的模型中通过输入特征变量预测站点i在h时段的真实需求如下：
步骤S4中所述考虑每个站点容量限制和需求到达分布，计算单车的初始投放量包括：
某站点单个时段内被拒需求数量的估计：
给定站点j和时段h＝[t0 h,ts h]，按照泊松分布分别由h时段的借车需求量概率分布和还车需求量概率分布产生随机值和/>则h内单车净变化量/>设站点容量为C，在初始时刻t0 h站点单车数量为q的条件下，h时段内产生的被拒借车需求数量的期望值表示为：
在式(15)中，u1表示被拒借车需求数量；
h时段内产生的被拒还车需求数量的期望值表示为：
在式(16)中，u2表示被拒还车需求数量；
在初始单车投放量为q的条件下，h时段内被拒需求数量的期望值为：
分别计算q从0到C的取/>最小值对应的q为h时段站点j的最优初始投放量。
2.如权利要求1所述的共享单车需求预测与投放调度方法，其特征在于：
某站点连续时段内被拒需求数量的估计：
若测算连续n个时段h1,..hm,..hn内站点j产生被拒需求数量的期望值，hm的初始库存是hm-1的末点库存，末点库存有三种情况：
h时段站点的单车库存量ph由q变为0的概率：
h时段站点的单车库存量ph由q变为C的概率：
h时段站点的单车库存量ph由q变为q’的概率：
初始单车投放量为q的条件下，站点在连续n个时段内产生的被拒需求数量的期望值为E(δf)：
取E(δf)最小值对应的q为该站点连续n个时段的最优初始投放量。
3.如权利要求1所述的共享单车需求预测与投放调度方法，其特征在于：步骤S5中所述划分城市的调度分区包括：
S51、首先筛选出净借需求量绝对值过大的站点，即极端需求站点，依次以这些站点为中心，搜索其附近1公里以内的全部站点，分别计算这些站点与中心站点间的距离，从最近的站点开始收纳入“零需”站点集，当站点集的净借需求量总和已经为零或者略大于零，则终止搜寻，为下一个极端需求站点构建“零需”站点集；
S52、将周围站点全部吸纳仍未达到终止条件，则以当前站点集中心位置为虚拟站点，以集内站点的净借需求量总和为虚拟站点的净借需求量，重复上述过程，吸纳新的站点加入；
S53、若站点集容量达到阈值仍然未达成终止条件，则释放这个站点集，直接为下一个极端需求站点构建“零需”站点集；若某站点附近1公里以内并没有还未加入站点集的站点，则释放这个站点集，将这个站点加入“落单”站点集，直接进行下一个“零需”站点集的构建；
S54、完成极端需求站点的站点集之后，在剩下的站点中找出地理位置上的***站点，依次为其构建“零需”站点集，循环这一过程，直到未加入站点集的站点数量与落单站点集数量相同；
S55、将落单站点视为一个“零需”站点集，计算所有“零需”集站点集的地理中心位置，用一个虚拟站点代替一个站点集，以经度最高的虚拟站点为始，与周围的虚拟站点合并，当实际站点数超过40且区域内所有站点的单车库存之和处于该区域最优单车投放量的上界和下界之间时，停止合并，作为一个调度分区；
S56、重复S51-S54，完成城市内所有站点的调度分区划分。
4.如权利要求1所述的共享单车需求预测与投放调度方法，其特征在于：步骤S6中所述建立区域内单车调度模型，获得调度路径的最优选择包括：
规定区域内有若干共享单车站点，用集合S1＝{1，2，……，N}表示，由标准运输货车进行调度作业，运输货车的出发点用S2＝{0}表示，每辆货车的载量上限为Q，在整个调度作业过程中，每辆作业货车的实时载量e必须满足0≤e≤Q；
目标函数1：调运成本和碳排放的最小化：
目标函数2：用户不满惩罚最小化：
约束条件：
Lij＝Lki+qi,k∈{0,…,N-1},Rki＝1,i∈S1,j∈{2,…,N}∪{0}(27)
Lij≤M·Rij(28)
L0j＝Lj0＝0(30)
S表示调度货车的调度路径，S＝S1∪S2＝{0,1,2,…,N}，其中0表示调度货车的出发点，其他表示各站点；qi表示在站点投放的单车的数量，qi<0时表示在站点i取走-qi辆单车；Ri,j为0-1变量，如果从站点i到站点集j在货车的调度路线中，Ri,j＝1，否则Ri,j＝0；gi表示站点i在调度工作进行前的单车数量；fi表示站点i在调度工作完成后的单车数量；F(x,i)表示表示站点i的单车数量为x时，预期会产生的损失需求的数量；Lij表示调度货车从站点i行驶到站点j过程中所载单车数量；Q表示调度货车的的容量；Dij表示站点i到站点j之间的运输距离；Ce表示调度货车单位公里单位千克耗用的燃油量；Cd表示每产生一个被拒绝的需求的不满惩罚；Ci表示站点i的容量；M表示一个无限大的正值；B表示消除子回路的辅助变量。
5.如权利要求4所述的共享单车需求预测与投放调度方法，其特征在于：通过混合多目标粒子群优化算法获得调度路径的最优选择，包括：
S61、生成初始解集：计算区域内站点间的距离矩阵，随机选择一个站点作为路线首个站点，然后在距离该站点最近的4个站点中以等概率随机选取任意一个作为下一个目的地，然后重复这一过程，避免选择已被选择的站点，直到将所有站点加入路线；将生成的路线作为解Rij输入式(22)-(35)，调用Python的scipy.optimize求解得到各个站点处的最优装/卸单车数量qi，与Rij构成一个完整的初始解；以相同的方式形成其他初始解，构建初始解集；
S62、适应度函数设计及个体最优选择：通过比较粒子的适应度值来获得粒子的优势关系，适应度函数来自目标函数1和目标函数2，指定如下：
S63、构建非劣解集：定义帕累托支配关系:对于粒子z1和z2，如f1(z1)≤f1(z2)且f2(z1)<f2(z2)，或者f1(z1)<f1(z2)且f2(z1)≤f2(z2)，则粒子z1支配粒子z2；同样，如果f1(z2)≤f1(z1)且f2(z2)<f2(z1)，或者f1(z2)<f1(z1)且f2(z2)≤f2(z1)，则粒子z2支配粒子z1；否则，粒子z1和z2没有支配关系；
对于任意粒子z*，当且仅当f1(z*)≤f1(z)且f2(z*)<f2(z)或者f1(z*)<f1(z)且f2(z*)≤f2(z)时，粒子z*被称为非支配粒子或帕累托粒子，所有这些非支配粒子都记录在非劣解集REP中；
S64、选择全局最优选择：全局最佳粒子gbest是从REP中选择的，REP保存搜索过程中找到的所有非支配解；在搜索开始时，将在初始化阶段构造的非支配粒子全部添加到REP；在搜索过程中，将每次迭代中发现的当前非支配粒子与REP中的所有解一一比较；如果一个当前粒子被REP中的粒子支配，则将其丢弃；否则，可以将此类粒子附加到REP；若REP中有粒子被新成员支配，则将这些粒子从REP中删除；当REP超过其最大容量时，用自适应网格方法删除多余的解，保持解决方案在REP中分布良好；对于gbest的选择，将解决方案空间划分为多个相等的网格，其中具有较少粒子的网格具有较高的被选择机会；设调整参数为β，网格数为n，网格j包含的粒子数量由nj表示，则第i个网格中的解被选择为gbest的概率为：
S65、交叉算子设计：使用遗传算法的交叉算子来更新算法中的粒子：首先由gbest更新粒子，然后由gbest更新粒子，粒子交叉算子其中Xi(t+1)和Xi(t)分别表示第i个粒子在当前迭代(t)和下一次迭代(t+1)中的位置；pbest(t)是第i个粒子的历史最优位置，gbest(t)是总体的全局最优位置；/>是交叉运算符；
S66、变异算子设计：使用结合VNS的变异算子对粒子进行更新，采用三种邻域结构：
(1)节点***：随机选择一个站点并将其从其原始位置移除，然后将其***另一个随机位置；
(2)节点交换：随机选择两个站点，并交换其位置；
(3)随机选择两个站点，将它们之间的顺序以相反的顺序排列；
以上三个邻域结构随机执行，进行变异操作后，如果新粒子支配了原始粒子，将替换原始粒子；如果原始粒子占主导地位，则丢弃新粒子；如果这两个粒子之间没有支配关系，则随机选择其中之一。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211665231.5A CN115936240B (zh) | 2022-12-23 | 2022-12-23 | 共享单车需求预测与投放调度方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211665231.5A CN115936240B (zh) | 2022-12-23 | 2022-12-23 | 共享单车需求预测与投放调度方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115936240A CN115936240A (zh) | 2023-04-07 |
| CN115936240B true CN115936240B (zh) | 2024-03-15 |
Family
ID=86555414
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211665231.5A Active CN115936240B (zh) | 2022-12-23 | 2022-12-23 | 共享单车需求预测与投放调度方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115936240B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116186571B (zh) * | 2023-04-25 | 2023-07-07 | 北京阿帕科蓝科技有限公司 | 车辆聚类方法、装置、计算机设备及存储介质 |
| CN116824861B (zh) * | 2023-08-24 | 2023-12-05 | 北京亦庄智能城市研究院集团有限公司 | 基于城市大脑平台多维数据的共享单车调度方法和*** |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107146334A (zh) * | 2017-05-09 | 2017-09-08 | 重庆壹元电科技有限公司 | 共享移动电源的供需平衡调度方法及*** |
| CN110119838A (zh) * | 2019-04-17 | 2019-08-13 | 成都信息工程大学 | 一种共享单车需求预测***、方法及装置 |
| CN111598481A (zh) * | 2020-05-26 | 2020-08-28 | 山东交通学院 | 共享单车流动***、基于子区划分的自动调度***及方法 |
| CN113642758A (zh) * | 2021-06-08 | 2021-11-12 | 上海汽车集团股份有限公司 | 一种共享电动汽车的供需预测及调度方法和装置 |
| CN113869715A (zh) * | 2021-09-26 | 2021-12-31 | 上海钧正网络科技有限公司 | 站点分组方法、归属站点确定方法、装置和计算机设备 |
| CN114240180A (zh) * | 2021-12-20 | 2022-03-25 | 北京工业大学 | 一种基于混合蚁群算法的多企业共享单车联合调度方法 |
| CN114266168A (zh) * | 2021-12-28 | 2022-04-01 | 武汉杰创达科技有限公司 | 一种共享节点负荷调度参数优选方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112907188B (zh) * | 2021-03-12 | 2022-05-24 | 北京化工大学 | 一种基于自适应邻域搜索算法的共享单车搬运优化方法 |
-
2022
- 2022-12-23 CN CN202211665231.5A patent/CN115936240B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107146334A (zh) * | 2017-05-09 | 2017-09-08 | 重庆壹元电科技有限公司 | 共享移动电源的供需平衡调度方法及*** |
| CN110119838A (zh) * | 2019-04-17 | 2019-08-13 | 成都信息工程大学 | 一种共享单车需求预测***、方法及装置 |
| CN111598481A (zh) * | 2020-05-26 | 2020-08-28 | 山东交通学院 | 共享单车流动***、基于子区划分的自动调度***及方法 |
| CN113642758A (zh) * | 2021-06-08 | 2021-11-12 | 上海汽车集团股份有限公司 | 一种共享电动汽车的供需预测及调度方法和装置 |
| CN113869715A (zh) * | 2021-09-26 | 2021-12-31 | 上海钧正网络科技有限公司 | 站点分组方法、归属站点确定方法、装置和计算机设备 |
| CN114240180A (zh) * | 2021-12-20 | 2022-03-25 | 北京工业大学 | 一种基于混合蚁群算法的多企业共享单车联合调度方法 |
| CN114266168A (zh) * | 2021-12-28 | 2022-04-01 | 武汉杰创达科技有限公司 | 一种共享节点负荷调度参数优选方法 |
Non-Patent Citations (2)
| Title |
|---|
| 共享单车需求预测及调度方法研究;史越;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;第13页第1段至第83页最后1段 * |
| 基于用地分布特征的共享汽车需求估计模型;孙立山;王顺超;罗薇;姚丽亚;王琨;;重庆交通大学学报(自然科学版);20200515(第05期);第2页第1栏第3段至第6页第1栏第2段 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115936240A (zh) | 2023-04-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115936240B (zh) | 共享单车需求预测与投放调度方法 | |
| CN111862579B (zh) | 一种基于深度强化学习的出租车调度方法及*** | |
| Simic et al. | An integrated CRITIC and MABAC based type-2 neutrosophic model for public transportation pricing system selection | |
| CN103632212B (zh) | 一种时变用户均衡动态网络演化客流预测***和方法 | |
| Pourghader Chobar et al. | A novel multi-objective model for hub location problem considering dynamic demand and environmental issues | |
| Caggiani et al. | A modular soft computing based method for vehicles repositioning in bike-sharing systems | |
| Caggiani et al. | A dynamic simulation based model for optimal fleet repositioning in bike-sharing systems | |
| CN110288212B (zh) | 基于改进的mopso的电动出租车新建充电站选址方法 | |
| CN112434849B (zh) | 基于改进多目标算法的危险品运输路径动态规划方法 | |
| CN108256969B (zh) | 一种公共自行车租赁点调度区域划分方法 | |
| CN108038578B (zh) | 基于需求预测和中心辐射网络的公共自行车静态调度方法 | |
| Guo et al. | Industrial information integration method to vehicle routing optimization using grey target decision | |
| CN104933474B (zh) | 危险化学品运输的模糊双层优化方法 | |
| Xu et al. | Dynamic bicycle scheduling problem based on short-term demand prediction | |
| CN111882909A (zh) | 基于双层规划的应急救援调度与动态路径集成方法 | |
| CN114139794A (zh) | 基于停车大数据的区域公共停车场充电需求预测方法 | |
| Liu et al. | Exploring the impact of spatiotemporal granularity on the demand prediction of dynamic ride-hailing | |
| Musolino et al. | A modelling framework to simulate paths and routes choices of freight vehicles in sub-urban areas | |
| Ensafian et al. | Cost-optimal deployment of autonomous mobile lockers co-operating with couriers for simultaneous pickup and delivery operations | |
| Jun et al. | Evolutionary neural network for learning of scalable heuristics for pickup and delivery problems with time windows | |
| CN117408598A (zh) | 基于动态路网的高速公路大件运输路径快速决策方法 | |
| Li et al. | A Simulation-optimization framework of vehicle relocation for one-way electric carsharing systems | |
| Chen et al. | Determining the optimal layout design for public bicycle system within the attractive scope of a metro station | |
| CN111709663A (zh) | 一种基于大数据的电动汽车充电场站选址方法 | |
| CN111091286A (zh) | 一种公共自行车调度模型与求解方法 |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |