#------------------------------------------------------------------------------------------------------------------------------- # Name: Basic-Data # Purpose: # # Author: Gaoy # # Created: 09/10/2015 # Copyright: (c) Gaoy 2015 # Licence: #------------------------------------------------------------------------------------------------------------------------------- # This file consists of two parts. The first part includes basic parameters, such as station number, services number, # link running times, station dwelling times, extra acc/dec times, headway, passenger arrival rate, OD matrix and so on. # The second part simulates the passenger flow in the disruption period, which generates the numbers of passengers # stranded in stations and in trains in the end of the disruption. # Note that in Python, array index starts from "0"!!! import numpy as np import scipy as sp import xlwt stationNum=13 #Station amount M=5000.0 #A large value M1=5000.0 M2=5000.0 START=0 #the beginning of the recovery is set as the "START" Maxdt=60 Delta1=300 Delta2=45 # slack time in the model W=900 #900 seconds or 15 minutes Gap=0.001 #Gurobi Gap SK1=1 SK2=2 SK3=3 # number of stations skipped by a service xlength=5000 # This is the length of x axis, for plotting H1=0 #disruption starting time: 7:00+H1*15 H2=2 #recovery starting time: 7:00+H2*15 Disruption=(H2-H1)*W #disruption period trainNum=20 # number of trains existingService=20 # there are 10 upstream existing services and 10 downstream existing services serviceNum=35 # we consider (serviceNum-10)=25 services in each direction, that is, 10 existing services and 15 future services sectionlen=[2631,1275,2366,1982,993,1728,1090,1355,2337,2301,2055,1281] # the length of each section accExtra=25 decExtra=25 # extra acceleration/deceleration time headway=90 # minimum headway backturning=120 # minimum back-turning time linkrun10=[190,108,157,135,90,114,103,104,166,150,140,102] # pure running times in upstream direction linkrun20=[100,141,150,162,103,101,111,90,135,157,105,195] # pure running times in downstream direction stationdt=[30,35,30,30,35,30,30,30,30,30,30,35,45] stationdt.extend([40,35,30,30,30,35,30,30,35,30,30,35,30]) # minimum dwelling time in each station N=stationNum # N is short for "stationNum" Capacity=1200 # train capacity #----------------------------------------------Generate passenger arrival rate------------------------------------------- #------------------1----2---3---4---5---6---7---8---9--10--11--12---13 7:00-10:00 PsnVol0=np.array([[315,420,543,630,720,795,618,540,471,420,360,300,300], #1 [495,600,750,732,810,930,939,780,663,540,444,390,402], #2 [486,540,570,630,600,618,570,420,390,354,282,300,288], #3 [1170,1440,1818,2190,2280,2187,1920,1812,1710,1470,1500,1200,972], #4 [756,900,1014,1080,1125,1170,1032,1020,930,780,720,612,510], #5 [315,357,354,390,435,552,495,441,390,381,312,279,228], #6 [273,330,393,480,549,663,510,438,378,339,303,255,216], #7 [222,303,360,420,516,555,450,390,360,339,303,288,243], #8 [957,1020,1170,1344,1416,1518,1440,1341,1170,1059,906,756,642], #9 [2355,2556,2766,2862,2796,2835,2850,2715,2541,2196,1857,1590,1314],#10 [585,672,744,870,846,744,660,591,603,510,426,303,267], #11 [180,222,249,171,147,216,264,159,138,240,168,96,90], #12 [816,897,1077,1170,1134,1209,1245,1119,891,840,702,648,510]]) #13 # number of passengers arriving in 15 minitues PsnVol0=(1./W)*PsnVol0# This is the arrival rate in per second #----------------------Construction of OD matrix------------------------ # 1 2 3 4 5 6 7 8 9 10 11 12 13 OD0=np.array([[0, 5, 5, 5, 5, 5, 5, 5, 5,20,20,20,5], #1 [5, 0, 5, 5, 5, 5, 5, 5, 5,20,20,20,5], #2 [10,5, 0, 5, 5, 5, 5, 5, 5,20,20,20,5], #3 [15,5, 5, 0, 5, 5, 5, 5, 5,20,20,20,5], #4 [15,5, 5, 5, 0, 5, 5, 5, 5,20,20,20,5], #5 [15,5, 5, 5, 5, 0, 5, 5, 5,15,20,20,5], #6 [15,10, 5, 5, 5, 5, 0,5, 5,15,15,20,5], #7 [20,10, 5, 5, 5, 5, 5,0, 5,10,15,20,5], #8 [20,15,5, 5, 5, 5, 5, 5, 0, 5,10,20,5], #9 [20,10,10,5, 5, 5, 5, 5, 5, 0,5,15,10], #10 [20,15, 5,5, 5, 5, 5, 5, 5, 5, 0, 5,5], #11 [20,10, 5,5, 5, 5, 5, 5, 5, 5, 5, 0,5], #12 [20,10,10,5, 5, 5, 5, 5, 5, 5, 5, 5,0]]) #13 #heterogenous OD #OD0=np.ones([N,N]) #for i in range(N): # OD0[i,i]=0 #homogenous OD temp=OD0.sum(1) OD1=np.zeros([N,N]) for i in range(N): OD1[i,]=np.divide(OD0[i,],float(temp[i])) OD2=np.zeros([2*N,2*N]) for i in range(N): for j in range(i,N): OD2[i,j]=OD1[i,j] for i in range(N,2*N): for j in range(i,2*N): OD2[i,j]=OD1[2*N-1-i,2*N-1-j] #-------------------------------------OD2 is the OD matrix we will use------------------------------------ #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- #--------------------------Simulation: generate the value of PsnInTrain and PsnAtStation at the end of disruption----------- #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- #------------------------------------the planned timetable before the disruption------------------------------------ OringinDepInt=210 # departure interval before the disruption PlanedTimetableArr=np.zeros((40,2*stationNum)) PlanedTimetableDep=np.zeros((40,2*stationNum)) StdTimetableArr=np.zeros(2*stationNum) StdTimetableDep=np.zeros(2*stationNum) StdTimetableDep[0]=StdTimetableArr[0]+stationdt[0] for i in range(1,N): StdTimetableArr[i]=StdTimetableDep[i-1]+linkrun10[i-1] StdTimetableDep[i]=StdTimetableArr[i]+stationdt[i] StdTimetableArr[N]=StdTimetableDep[N-1]+backturning StdTimetableDep[N]=StdTimetableArr[N]+stationdt[N] for i in range(N+1,2*N): StdTimetableArr[i]=StdTimetableDep[i-1]+linkrun20[i-1-N] StdTimetableDep[i]=StdTimetableArr[i]+stationdt[i] for k in range(40): for i in range(2*N): PlanedTimetableArr[k,i]=k*210+StdTimetableArr[i] PlanedTimetableDep[k,i]=k*210+StdTimetableDep[i] #----------------------------passenger flow simulation before the disruption------------------------ nWantOn=np.zeros([2*N,2*N])# numbers of passengers arriving at station i with destination j during OringinDepInt(210s) in the normal situation before the disruption nWantOn2=np.zeros([2*N,2*N])#passengers stranded in station i with destination j in the beginning of the recovery period PsnAtStationStart=np.zeros(2*N) # passenger number stranded in stations PsnInTrain=np.zeros([2*N,2*N])# the number of passengers from i to j in a service when the service is departing from a station in the normal situation before the disruption PsnInTrain0=np.zeros([existingService,2*N]) # passengers in existing services in the beginning of the recovery period for i in range(N): for j in range(i+1,N): nWantOn[i,j]=PsnVol0[i,H1]*OringinDepInt*OD2[i,j] for i in range(N,2*N): for j in range(i+1,2*N): nWantOn[i,j]=PsnVol0[2*N-1-i,H1]*OringinDepInt*OD2[i,j] nCanOn=min(Capacity,sum(nWantOn[0]))#passenger number in service when the service departs from station 0 nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[0])# the ratio between "nCanOn" and the sum of "nWantOn" PsnInTrain[0]=[round(nWantOn[0,j]*nCanOn_VS_nWantOn) for j in range(2*N)]# the number of passengers in a service when the service is departing from station 0 for i in range(1,N-1): nOff=PsnInTrain[i-1,i] nCanOn=min(Capacity-sum(PsnInTrain[i-1])+nOff,sum(nWantOn[i])) nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[i]) PsnInTrain[i]=[0*j for j in range(i+1)]+[PsnInTrain[i-1,j]+round(nWantOn[i,j]*nCanOn_VS_nWantOn) for j in range(i+1,2*N)] nCanOn=min(Capacity,sum(nWantOn[N]))#passenger number in train when the train departs from station N nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[N]) PsnInTrain[N]=[round(nWantOn[N,j]*nCanOn_VS_nWantOn) for j in range(2*N)] for i in range(N+1,2*N-1): nOff=PsnInTrain[i-1,i] nCanOn=min(Capacity-sum(PsnInTrain[i-1])+nOff,sum(nWantOn[i])) nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[i]) PsnInTrain[i]=[0*j for j in range(i+1)]+[PsnInTrain[i-1,j]+round(nWantOn[i,j]*nCanOn_VS_nWantOn) for j in range(i+1,2*N)] # In the normal situation, a passenger can always get on the first coming service. #-----------------------------passenger flow simulation during disruption------------------------ startStatus=np.zeros((existingService)) startStatus[:existingService]=[ 25, 24, 23, 22, 21, 19, 18, 16, 15, 13, 12, 11, 10, 9, 8, 6, 5, 3, 2, 1] # 'startStatus' is the locations (station or next station) of the existing services when the disruption is over timeInstance=np.ones(existingService)*(-15) # timeInstrance[k]<0, means that the train has arrived at corresponding station abs(timeInstrance[k]) ago. # In the experiments, it is assumed that all the existing services are dwelling at stations during the disruption period. # So, timeInstrance[k] can be any negative value. #---------------------------downstream existing services--------------------------- temp=2*N-2 # station for k in range(10): if startStatus[int(k)]==25: continue; for i in range(temp,int(startStatus[k]),-1): for j in range(i+2,2*N):#pay attention to i+2 nWantOn[i,j]=PsnVol0[2*N-1-i,H1]*Disruption*OD2[i,j] nWantOn2[i,j]=nWantOn[i,j] PsnAtStationStart[i]=sum(nWantOn2[i]) temp=int(startStatus[k])-1 i=int(startStatus[k]) for j in range(i+2,2*N): nWantOn[i,j]=PsnVol0[2*N-1-i,H1]*Disruption*OD2[i,j] if i<2*N-1: nOff=PsnInTrain[i-1,i]+PsnInTrain[i-1,i+1] else: nOff=PsnInTrain[i-1,i] #number of passengers getting off nCanOn=min(Capacity-sum(PsnInTrain[i-1])+nOff,sum(nWantOn[i])) nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[i]) PsnInTrain0[k]=[0*j for j in range(i+2)]+[PsnInTrain[i-1,j]+round(nWantOn[i,j]*nCanOn_VS_nWantOn) for j in range(i+2,2*N)] nWantOn2[i]=[nWantOn[i,j]*(1-nCanOn_VS_nWantOn) for j in range(2*N)] PsnAtStationStart[i]=sum(nWantOn2[i]) if temp>=13: for i in range(temp,12,-1): for j in range(i+2,2*N): nWantOn[i,j]=PsnVol0[2*N-1-i,H1]*Disruption*OD2[i,j] nWantOn2[i,j]=nWantOn[i,j] PsnAtStationStart[i]=sum(nWantOn2[i]) #--------------------------------upstream existing services----------------------- temp=N-2 for k in range(10,existingService): if startStatus[k]==13: continue; if startStatus[k]==12: continue; for i in range(temp,int(startStatus[k]),-1): for j in range(i+2,2*N): nWantOn[i,j]=PsnVol0[i,H1]*Disruption*OD2[i,j] nWantOn2[i,j]=nWantOn[i,j] PsnAtStationStart[i]=sum(nWantOn[i]) temp=int(startStatus[k])-1 i=int(startStatus[k]) for j in range(i+2,N): nWantOn[i,j]=PsnVol0[i,H1]*Disruption*OD2[i,j] if i=0: for i in range(temp,-1,-1): for j in range(i+2,N): nWantOn[i,j]=PsnVol0[i,H1]*Disruption*OD2[i,j] nWantOn2[i,j]=nWantOn[i,j] PsnAtStationStart[i]=sum(nWantOn2[i])