假设地理空间
为了简化问题,便于验证本代码写作的准确性和有效性,故此做以下假设空间。
论文章节
论文:王露. 城市纯电动汽车快速充电设施的布局选址优化模型研究[D]. 2016.
这篇文章功底很扎实,通篇内容详尽,可见王露同学的深厚理工科背景。
复现代码
下面展示一些 内联代码片
。
import numpy as np
def driving_time_one_moment(W_I,M_J,TRAVELTIME_IJ,Y_IJ):
print(np.sum(np.dot(np.dot(W_I,Y_IJ),TRAVELTIME_IJ.T)))
return np.sum(np.dot(np.dot(W_I,Y_IJ),TRAVELTIME_IJ.T))
def waiting_time_one_moment(W_I,M_J,TRAVELTIME_IJ,Y_IJ):
tc=8 # 服务时间,24小时提供服务;
tf=2 #每辆车充电时长;
# 第一步:计算用户达到率
Tao_J=np.dot(W_I,Y_IJ/tc)
# 第二步:单位时间内充电站平均服务能力
U_J=M_J/tf
# 第三步:充电站排队系统服务强度:,由于ROU_J会存在大于1的情况,从而会促使后面求解
ROU_J=Tao_J/U_J
# 第四步:充电站内充电桩全部空闲概率:
P_J=np.full(shape=(1,M_J.shape[1]),fill_value=0.0)
for j in range(M_J.shape[1]):
temp = 0
for k in range(M_J.shape[1]):
temp+=(np.power(M_J[0,j]*ROU_J[0,j],k))/(np.math.factorial(k))
p_j0=1/(temp+(np.power(M_J[0,j]*ROU_J[0,j],M_J[0,j]))/(np.math.factorial(M_J[0,j])*(1-ROU_J[0,j])))
P_J[0,j]=p_j0
# 第五步:计算排队等候时间期望
W_Jq=np.full(shape=(1,M_J.shape[1]),fill_value=0.0)
for j in range(M_J.shape[1]):
w_jq=(np.power(M_J[0,j],M_J[0,j])*np.power(ROU_J[0,j],M_J[0,j]+1)*P_J[0,j])/(ROU_J[0,j]*np.math.factorial(M_J[0,j])*np.power(1-ROU_J[0,j],2))
W_Jq[0,j]=w_jq
# 第六步:所有用户的总的等待花费时间
T2=0
for j in range(M_J.shape[1]):
T2+=W_Jq[0,j]*Tao_J[0,j]*tc
print(T2)
if __name__ == '__main__':
# 需求量;
W_I = np.array([[10, 20, 30, 50]])
# 充电桩供给量;
M_J = np.array([[5, 10, 15]])
# 行是demand,列是provider
TRAVELTIME_IJ = np.array([[1, 1.5, 2.5],
[1.5, 1, 2],
[1, 1.5, 1],
[2.5, 1.5, 1]])
# 出行时间的阻尼函数,衰减函数
F_DIJ = 1 / TRAVELTIME_IJ
Sum_Dij_I = np.sum(F_DIJ, axis=1)
# 计算选择权重
Y_IJ = np.full(shape=(TRAVELTIME_IJ.shape), fill_value=0.0)
for i in range(W_I.shape[1]):
Y_IJ[i, :] = F_DIJ[i, :] / Sum_Dij_I[i]
# 计算驾车时间
driving_time_one_moment(W_I,M_J,TRAVELTIME_IJ,Y_IJ)
# 计算等候时间
waiting_time_one_moment(W_I,M_J,TRAVELTIME_IJ,Y_IJ)
运行结果
662.9544044665013
-7.133069261540492
问题思考
排队时间会出现负值?这种情况下,应当如何解释?
问题反思与重定位
查到了三个原因;
原因1:以上代码有一点小错误(详细改正见下面的代码);
原因2:因为公式需要出一个u;
原因3:仍旧会存在超出的部分,需要加判断。
原因1和原因3在代码中体现。
查找了排队论PDF之后,发现里面有一句话是这样说的:
所以对求解空闲概率时候的rou做了修改,公式如下:
【该手稿来自Jing Huang同学,以上的改进发现也来自于Jing Huang,非常细致】
在此基础上做测试,发现人会数值不稳定的情况,这是因为还存在被抵消的情况,因此需要对rou_s加上进一步的判断。
def waiting_time_one_moment_2(W_I,M_J,TRAVELTIME_IJ,Y_IJ):
tc=6 # 服务时间,24小时提供服务;
tf=2 #每辆车充电时长;
# 第一步:计算用户达到率
Tao_J=np.dot(W_I,Y_IJ/tc)
# 第二步:单位时间内充电站平均服务能力
U_J=M_J/tf
# 第三步:充电站排队系统服务强度:,由于ROU_J会存在大于1的情况,从而会促使后面求解
ROU_J=Tao_J/U_J
print("ROU_J", ROU_J)
print("ROU_J/M_J",ROU_J/M_J)
# 加一个判断
if np.any(ROU_J/M_J>1):
# 当ROU_J/M_J<1条件下才能得到,即要求顾客的平均到达率小于系统的平均服务率,才能使系统达到统计平衡。
T2=np.inf
print(T2)
return T2
# 第四步:充电站内充电桩全部空闲概率:
P_J=np.full(shape=(1,M_J.shape[1]),fill_value=0.0)
for j in range(M_J.shape[1]):######修改
temp = 0
for k in range(M_J.shape[1]):#修改
temp+=(np.power(ROU_J[0,j],k))/(np.math.factorial(k))
p_j0=1/(temp+(np.power(ROU_J[0,j],M_J[0,j]))/(np.math.factorial(M_J[0,j])*(1-ROU_J[0,j]/M_J[0,j]))) ##修改
P_J[0,j]=p_j0
print("P_J",P_J)
# 第五步:计算排队等候时间期望
W_Jq=np.full(shape=(1,M_J.shape[1]),fill_value=0.0)
for j in range(M_J.shape[1]):
w_jq=(np.power(ROU_J[0,j],M_J[0,j]+1)*P_J[0,j])/(M_J[0,j]*Tao_J[0,j]*np.math.factorial(M_J[0,j])*np.power(1-ROU_J[0,j]/M_J[0,j],2))#修改
W_Jq[0,j]=w_jq
# 第六步:所有用户的总的等待花费时间
T2=0
for j in range(M_J.shape[1]):
T2+=W_Jq[0,j]*Tao_J[0,j]*tc
print(T2)
if T2<0:
T2=np.inf
def calculate_Y_IJ(TRAVELTIME_IJ,W_I):
# 出行时间的阻尼函数,衰减函数
F_DIJ = 1 / TRAVELTIME_IJ
Sum_Dij_I = np.sum(F_DIJ, axis=1)
# 计算选择权重
Y_IJ = np.full(shape=(TRAVELTIME_IJ.shape), fill_value=0.0)
for i in range(W_I.shape[1]):
Y_IJ[i, :] = F_DIJ[i, :] / Sum_Dij_I[i]
return Y_IJ
def single_time():
# 需求量;
W_I = np.array([[10, 20, 30, 50]]) # 4.7
# 充电桩供给量;
M_J = np.array([[1, 20, 3]]) # 会出现负值
M_J = np.array([[1, 1, 1]]) # 是一个BUG,因为除以无效
M_J = np.array([[2, 2, 2]]) # 是一个BUG
M_J = np.array([[4, 4, 10]]) # 是一个边缘值
# 行是demand,列是provider
TRAVELTIME_IJ = np.array([[1, 1.5, 2.5],
[1.5, 1, 2],
[1, 1.5, 1],
[2.5, 1.5, 1]])
# 出行时间的阻尼函数,衰减函数
F_DIJ = 1 / TRAVELTIME_IJ
Sum_Dij_I = np.sum(F_DIJ, axis=1)
# 计算选择权重
Y_IJ = np.full(shape=(TRAVELTIME_IJ.shape), fill_value=0.0)
for i in range(W_I.shape[1]):
Y_IJ[i, :] = F_DIJ[i, :] / Sum_Dij_I[i]
# 计算驾车时间
driving_time_one_moment(W_I, M_J, TRAVELTIME_IJ, Y_IJ)
# 计算等候时间
waiting_time_one_moment_2(W_I, M_J, TRAVELTIME_IJ, Y_IJ)
if __name__ == '__main__':
single_time()