求解驾驶时间和等待时间，复现【城市纯电动汽车快速充电设施的布局选址优化模型研究】论文中的排队模型

假设地理空间

为了简化问题，便于验证本代码写作的准确性和有效性，故此做以下假设空间。

论文章节

论文：王露. 城市纯电动汽车快速充电设施的布局选址优化模型研究[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()