公共自行车站点管理

公共自行车服务点管理

摘要

本题主要是针对公共自行车服务点出现空架情况而利用数学原理建立可靠模型。通过对服务点在借、还车高峰期的观测以及工作人员安排调度等各方面因素的综合考虑，建立了基于中心极限定理的服务点车量数的指标等级预测模型。在调度问题上，要使得公共自行车服务点尽可能保持理想状态，建立了最优回路的模型。

第一问：根据题目要求采用了数据拟合法，对选定的具有代表性的服务点的一天的借、还数量进行统计，得出的数据进行拟合。根据最终的拟合曲线可得出，工作人员应当在在8点30至9点时间段对车辆进行上架操作，17点30分至18点进行下架操作。接着我们用SAS求出样本均值和标准差。再利用中心极限定理将指标等级划分为：低、较低、理想、较高、高五个等级，最后利用MATLAB编程对上架下架车数量进行处理得到最佳上架数量为4，下架数量为3。

第二问：根据杭州市公共自行车借还的实际情况，将需要调整的站点根据分布密集度分为三个区域，每个区域分配两名人手；然后将公共自行车服务点间自行车的调度问题视为一个特殊的旅行商问题模型，再将各站点位置按照已划分好的区域位置在地图上进行标注。再利用MATLAB编程最终确定出最优回路分别为：T=5158-508 8-5065-5371-5211-5198-5287-5131-5169-5147-5365-5307 1

-5287-5073-5158

T=5332-5354-5298-5346-5327-5339-5332

2

T=5239-5251-5100-5180-5121-5107-5239

3

最终建立优化模型求解出每个站点上架下架车的数量。

关键词：数据拟合法、最优回路、TSP、区域划分、旅行商问题

一、问题重述

2008年5月，杭州在全国首推公共自行车服务系统，以缓解行路停车难、解决“最后一公里”交通问题。截至2012年底，杭州共有公共自行车 69750 辆，公共自行车服务点 2962 个，遍布杭城大街小巷，日均租车量达到 25.75 万人次。

每个服务点完成公共自行车的借车及还车环节。每个服务点设有若干车架，每一车架只能锁住一辆自行车。当服务点所有车架满架时，只能借车而不能还车；当服务点所有车架空架时，不能借车而只能还车。当车架既不满架也不空架时即为理想状态。当车架不为理想状态时，需要人工对车辆进行调配，以满足出行人的需要。因此需要为工作人员配置一台调度车，通过上、下架使得服务点恢复理想状态。

公共自行车是一项公益性事业，在向市民提供优质、低价服务的同时必须考虑运营成本。对于无专人值守的公共自行车服务系统，可以通过工作人员驾驶调度车巡逻的方式对服务点进行管理，使车架在尽可能多的时段内处于理想状态，提升市民的满意度。试建立数学模型回答以下问题：

1．假设某服务点车架数已知，通过实地调查统计掌握了一段时期内从早上6 时至晚上 9 时间隔 10 分钟的借车、还车需求数据。试给出巡逻人员到该服务点的时刻及上架或下架自行车的数量，使得一天内需对该点进行管理的次数最少。若规定一天内到该服务点的次数不超过二次，如何才能使服务点处于理想状态的时段尽可能长。

2．以杭州西湖区为例，请根据杭州公交出行实时信息服务系统，确定巡逻管理该区域的工作人员数量及调度方案，使该区域内的公共自行车服务点尽可能保持理想状态，且人力成本较小。

二、问题分析

问题一：

根据题目要求，已知服务点为选定的一点，结合杭州市租车点较多、日均租车量较高等因素进行综合考虑，由此我们决定以住宅区附近的服务点为例，进行自行车流动量统计。由此根据杭州公交出行实时信息服务系统选定了编号为1015服务点为调查研究对象，对该服务点从早上6点至晚上9点进行借车、还车数量的需求量统计。统计结果如下：

表2.1 服务点当天的借车、还车需求数据

第一问同时要求同时结合上述数据和人员调配，给出一个成本低廉、服务点处于理想状态最大得出在最佳时间点的上、下架车辆数，使得该服务点处于理想

化的流动方案。根据对上表中数据的处理及分析后，可利用SAS求解，同时通过马尔科夫链数学原理进行等级预测将指标等级划分，再结合MATLAB进行计算，状态的时间尽可能长，且成本低廉。

问题二：

题目中规定以杭州西湖区为例，根据杭州公交出行实时信息服务系统查询到该片区的信息，确定巡逻管理该区域的工作人员数量及调度方案，使该区域内的

公共自行车服务点尽可能长时间保持理想状态，且人力成本较小。

根据查询的资料可知，杭州西湖区共有501处设有公共自行车站点，我们选

取了某时刻杭州西湖区各个站点的实时可借、还数据，统计了该时刻不处于理想

状态的站点，并按照一定的比例，对真实环境距离进行按照一定比例还原，将公共自行车服务点间自行车的调度视为一个特殊的旅行商问题的模型即：在一个固定区域选择一个运输车，从一个固定点出发，经过每个站点收集多余的自行车或投放补给一定数量的自行车，最终回到原点的一条长度最短的回路T，然后建立优化模型求出上架下架数量，使西湖区各服务点尽可能长时间的在保持理想状态。

三、模型假设

1、假设自行车的流通量不因天气变化而变化；

2、公共自行车无损坏、被偷等情况出现；

3、每位市民是否租借自行车相互独立；

四、定义与符号说明

x样本均值，

S 标准差

b>表示服务点i上架i b量自行车

i

b<表示服务点i下架i b量自行车

i

ω表示上架下架车的总次数

五、模型的建立与求解

(一)问题一

模型准备

根据题目要求，对某一定点的一天的租车流动量进行统计，得到的数据经过分析后，利用MATLAB进行分析处理，得到该服务点一天内租车流动量拟合曲线，根据图像我们可以看出借车、还车环节在一天内各自的曲线走势。由此，我们可以通过计算，大致得出最佳上、下架时间点。

同样对于服务点上、下架的车辆数的求解，可利用SAS进行计算。利用马尔科夫链的指标等级预测将指标等级划分，再利用MATLAB进行进一步计算，得出在最佳时间点的上、下架车俩数。由此使得服务点处于理想状态的时间尽可能长。模型建立与求解

1、工作人员巡逻时间

利用MATLAB通过对调查数据分析后得到的数据进行拟合，得到如下所示的拟合曲线:

图1一天时间内不同时刻可租车数量

图2一天时间内不同时刻可还车数量

图3可借车可还车交叉图

经过对图1、图2、图3的分析，我们可以看出该服务点在一天中的借车、还车环节都各自存在高峰期。分析图1、图2 ，我们可以将一天借、还车环节的曲线走势大致区分为：早高峰为借车环节，晚高峰为还车环节。根据杭州当地普通居民的上班、下班时间综合考虑，不难判断出，在早晨7点30分到9点这个时间段内，借车数量急剧上增，即可借车数量急剧下降，并在9点后大致趋于稳定；在晚上17点30分到18点30分这个时间段内，还车数量急剧上增，即可还车数量急剧下降，并在19点左右趋势逐渐趋于平缓。根据这一现象，我们可以对工作人员对该点的巡逻时间进行合理安排，即在8点30至9点时间段对车辆进行上架操作，17点30分至18点进行下架操作。

2、上、下架车辆数量

通过SAS 计算服务点的样本均值x 、标准差 S ,再利用中心极限定理划分指标值序列，将其化为低、较低、理想、较高、高五个等级，确定理想状态的范围，再通过MATLAB 计算出上下架自行车的数量。

用SAS 求得1290x x x L 指标序列的相关数据如下

样本均值9.17x = ，标准差S=6.70 划分为五组：

11222211(,),(,),(,),(,),(,)x s x s x s x s x s x s x s x s αααααααα?∞????+++++∞ 1α可在[]1.0,1.5中取值，2α可在[]0.3,0.6中取值。在此处我们取1α=1.0 2α=0.3

则指标分级序列为：

表5.1指标分级序列

x 为服务点的车辆数 指标分级

x<3.10 低 3.10

高

由表可知：最理想状态为7.16

利用MATLAB 编程得到上架下架变更后的图像，在最佳时间点进行上架下架

操作，更改上架下架数量，观察数量不同下的图像变化趋势，能够使服务点的车辆数尽可能长时间的维持在理想状态。

我们得到在8点对车辆进行上架操作，上架数量为4，在17点30分对车辆进行下架操作，下架数量为3。得到如下图像：

图4变更后的曲线图

从图4我们可以看出，在8点对车辆进行上架操作上架数量为4，在17点30分对车辆进行下架操作下架数量为3后我们可以看出整个时间段绝大部分的车辆数都在3.10-15.87之间得到我们想要的理想区域。

(二)问题二

通过查阅杭州公交出行实时信息服务系统，我们统计了某时刻西湖区公共自行车服务点的位置、总车架数等信息以及有车车架数和无车车架数等实时数据（见附录一）。从中我们找出了未处于理想状态的站点，如下表所示：

表5.2 为处于理想状态站点

根据上表我们确定了需要调度的站点，根据实际观测各个站点的具体位置，

我们大致分为三个区域，这样我们给每个区域匹配两个人手，再根据杭州公交出行实时信息服务系统按照一定的比例画出这几个站点间的比例图形。

第一个区域是由5065,5073,5088,5131,5147,5158,5169,5198,5211,5225,5287，5307,5365,5371构成。比例图如下：

图5 区域一比例图

第二个区域是5298,5327,5332,5339,5346,5354构成。比例图如下：

图6 区域二比例图

第三个区域是5100,5107,5121,5180,5239,5251构成。比例图如下：

图7 区域三比例图

以上就是三个区域的结构图，根据各个站点以及他们之间的距离比例关系，我们利用图论中的最短路原理用MATLAB 编程，我们确定一条长度最短的回路T,即从某一点出发，遍及所有服务点当且仅当一次的最短距离为T。

第一个区域的回路1T 为：

5158-5088-5065-5371-5211-5198-5287-5131-5169-5147-5365-5307-5287-5073-5158

第一个区域的可借可还信息如下：

表5.3第一个区域的可借可还信息

第二个区域的回路2T 为：

5332-5354-5298-5346-5327-5339-5332

第二个区域的可借可还信息如下：

表5.4第二个区域的可借可还信息

第三个区域的回路3T 为：

5239-5251-5100-5180-5121-5107-5239

第三个区域的可借可还信息如下：

表5.5第三个区域的可借可还信息

对于每个区域回路各个站点的上架下架车数量我们建立优化模型。对于各服务点的集合为{}12,,......n v v v v =回路的顺序为(1,2,...,1)i n =，用0i b >表示服务

点i 上架i b 量自行车；0i b <表示服务点i 下架i b 量自行车。

则对每个区域而言，我们为了使人力资源最少，我们对上架下架自行车的数量进行控制，对此建立目标函数：

1

min n

i

b ω=∑

约束条件：

设每个站点的自行车的初值为i x ，为了使每个站点尽可能处于理想状态，根据问题一中的指标分级序列，我们对进行上架下架车后每个站点的自行车的总和进行限制，即316i i x b ≤+≤

运输车上的自行车量数应该大于等于下一站的补给需求，则1

1i i i b b ?≥∑

最终回到始发点的时候运输车上的自行车数为零，则所有服务点i b 满足

1

0n

i i b ==∑ 在三个区域分别针对此模型进行求解,我们得到每个区域的每个站点的上架下架车的数量如下表所示:

表5.6第一个区域各个站点的上架下架车数量

表5.7

第二个区域各个站点的上架下架车数量

表5.8第三个区域各个站点的上架下架车数量

六、模型评价

在第一问中我们运用了基于中心极限定理的服务点车量数的指标等级预测模型，该模型能够在已知数据较少的情况下较为准确的表达出该点的不同时刻的需求量和理想状态下的可租车两辆数，切准确度相对较高，从而判断出上下架车辆的标准。但是由于数据过少，且不能对实际的服务点分布进行更加详细的统计，因此不能更宏观的表现出该站点的理性状态的具体数值。

但该模型运用范围较广，延伸性较大，可以运用到最佳灾情巡视路线，售货员送货的最佳路线等。

七、参考文献

[1] 柳祖鹏，丁卫东，程逸旻．公共自行车系统站间调度优化研究［Ｊ］．城市

公共交通，2011（1）：39-42；．

[2] 姚 遥，周扬军．杭州市公共自行车系统规划［Ｊ］．城市交通．

[3] 郭敏辉，钟 明．上海市公共自行车系统规划与实践［Ｊ］．城市交通.

[4] 李黎辉，陈 华，孙小丽．武汉市公共自行车租赁点布局规划［Ｊ］．城市

交通，2009，7（4）：39-44，38．

[5] 田贵超,黎明，韦雪洁．旅行商问题（TSP）的几种求解方法［Ｊ］． 计算

机仿真，2006（8）：153－157．

[6] 周 品，赵新芬.MATLAB数学建模与仿真，北京市海淀区，国防工业出版社，

2009年

八、附录

附录一

附录二

第一题时间数量关系：

clc;

clear;

for i=1:90

t(i)=i;

end

s=[19,18,17,15,15,16,21,21,21,18,18,10,7,6,4,1,2,2,2,2,2,2,2,2,3,4,2,3,4,4,4,3,3,3,3,3, 1,1,1,1,1,2,2,5,7,7,6,6,7,7,7,7,7,7,7,3,8,6,6,8,8,6,6,6,6,6,5,8,10,10,13,15,17,18,19,21, 18,17,16,18,19,16,17,18,19,17,18,18,19,19];

ti=1:0.01:90;

% plot(t,s,'.');

% grid;

% si=interp1(t,s,ti,'spline');

% hold on

% plot(ti,si);

p=polyfit(t,s,8);

ti1=1:0.01:13;

si1=polyval(p,ti1);

ti2=13:0.01:70;

si2=polyval(p,ti2)+5;

ti3=70:0.01:90;

si3=polyval(p,ti3)-5;

hold on

plot(ti1,si1,ti2,si2,ti3,si3);

TSP问题遗传算法通用Matlab程序

附录三

程序一：主程序

%TSP问题（又名：旅行商问题，货郎担问题）遗传算法通用matlab程序

%D是距离矩阵，n为种群个数

%参数a是中国31个城市的坐标

%C为停止代数，遗传到第C代时程序停止,C的具体取值视问题的规模和耗费的时间而定

%m为适应值归一化淘汰加速指数,最好取为1,2,3,4,不宜太大

%alpha为淘汰保护指数,可取为0~1之间任意小数,取1时关闭保护功能,建议取0.8~1.0之间的值

%R为最短路径,Rlength为路径长度

function [R,Rlength]=geneticTSP(D,a,n,C,m,alpha)

[N,NN]=size(D);

farm=zeros(n,N);%用于存储种群

for i=1:n

farm(i,:)=randperm(N);%随机生成初始种群

end

R=farm(1,:);

subplot(1,3,1)

scatter(a(:,1),a(:,2),'x')

pause(1)

subplot(1,3,2)

plotaiwa(a,R)

pause(1)

farm(1,:)=R;

len=zeros(n,1);%存储路径长度

fitness=zeros(n,1);%存储归一化适应值

counter=0;

while counter for i=1:n

len(i,1)=myLength(D,farm(i,:));%计算路径长度end

maxlen=max(len);

minlen=min(len);

fitness=fit(len,m,maxlen,minlen);%计算归一化适应值rr=find(len==minlen);

R=farm(rr(1,1),:);%更新最短路径

FARM=farm;%优胜劣汰，nn记录了复制的个数

nn=0;

for i=1:n

if fitness(i,1)>=alpha*rand

nn=nn+1;

FARM(nn,:)=farm(i,:);

end

end

FARM=FARM(1:nn,:);

[aa,bb]=size(FARM);%交叉和变异

while aa if nn<=2

nnper=randperm(2);

else

nnper=randperm(nn);

end

A=FARM(nnper(1),:);

B=FARM(nnper(2),:);

[A,B]=intercross(A,B);

FARM=[FARM;A;B];

[aa,bb]=size(FARM);

end

if aa>n

FARM=FARM(1:n,:);%保持种群规模为n

end

farm=FARM;

clear FARM

counter=counter+1

end

Rlength=myLength(D,R);

subplot(1,3,3)

plotaiwa(a,R)

附录四

程序二：计算邻接矩阵

%输入参数a是中国31个城市的坐标

%输出参数D是无向图的赋权邻接矩阵

function D=ff01(a)

[c,d]=size(a);

D=zeros(c,c);

for i=1:c

for j=i:c

bb=(a(i,1)-a(j,1)).^2+(a(i,2)-a(j,2)).^2;

D(i,j)=bb^(0.5);

D(j,i)=D(i,j);

end

end

附录五

程序三：计算归一化适应值

%计算归一化适应值的子程序

function fitness=fit(len,m,maxlen,minlen)

fitness=len;

for i=1:length(len)

fitness(i,1)=(1-((len(i,1)-minlen)/(maxlen-minlen+0.0001))).^m;

end

附录六

程序四：交叉和变异的子程序

%交叉算法采用的是由Goldberg和Lingle于1985年提出的PMX(部分匹配交叉)

function [a,b]=intercross(a,b)

L=length(a);

if L<=10%确定交叉宽度

W=9;

elseif ((L/10)-floor(L/10))>=rand&&L>10

W=ceil(L/10)+8;

else

W=floor(L/10)+8;

end

p=unidrnd(L-W+1);%随机选择交叉范围，从p到p+W for i=1:W%交叉

x=find(a==b(1,p+i-1));

y=find(b==a(1,p+i-1));

[a(1,p+i-1),b(1,p+i-1)]=exchange(a(1,p+i-1),b(1,p+i-1));

[a(1,x),b(1,y)]=exchange(a(1,x),b(1,y));

end

function [x,y]=exchange(x,y)

temp=x;

x=y;

y=temp;

附录七

程序五: 计算路径的子程序

%该路径长度是一个闭合的路径的长度

function len=myLength(D,p)

[N,NN]=size(D);

len=D(p(1,N),p(1,1));

for i=1:(N-1)

len=len+D(p(1,i),p(1,i+1));

end

附录八

程序六：用于绘制路径示意图的程序

function plotaiwa(a,R)

scatter(a(:,1),a(:,2),'x')

hold on

plot([a(R(1),1),a(R(31),1)],[a(R(1),2),a(R(31),2)])

hold on

for i=2:length(R)

x0=a(R(i-1),1);

y0=a(R(i-1),2);

x1=a(R(i),1);

y1=a(R(i),2);

xx=[x0,x1];

yy=[y0,y1];

plot(xx,yy)

hold on

end

运算结果：

clear;

clc;

% a=[0,0;1.3,0;1.3,2.5;1.5,3.2;1.6,2.9;1.9,4]; %区域一

% a=[0,0;1.5,-1;1.5,1.2;2.1,0.8;3,-0.7;5,3]; %区域二

% a=[3,0;0.5,3;4,1;2,6;3,2;3,6.5;2,1;3.5,8;1.5,0;2,9;4,10;2.5,11;2,12;1,10]; %区域三D=ff01(a);

n=200;

C=100000;

m=2;

alpha=0.8;

[R,Rlength]=geneticTSP(D,a,n,C,m,alpha);