# In order to implement this programming, Python 2.6 should be installed with package numpy, scipy and matplotlib.
# Besides, Gurobi should be installed and configured.

# All the parameters and data are stored or generated in file Basicdata.py.
# Running this file, the rescheduled timetable will be obtained.
# Besides, the passenger volumn evolution in the upstream direction will be plotted.

# the downstream services and upstream services are numbered in different ways.
# downstream services are numbered as 0,1,2,...,24
# upstream services are numbered as 10, 11, 12,...,34

from gurobipy import *
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import math

from Basicdata import *# Basicdata.py must be put in the same folder.
from time import clock

# Creat optimization model
m = Model('timetable')


upOpt=-1
downOpt=0
# service in each direction

downOpt0=0
#the first service  in downstream direction

for k in range(existingService-1):
    if (startStatus[k]>stationNum or (startStatus[k]==stationNum and timeInstance[k]<0)) and (startStatus[k+1]<stationNum
    or (startStatus[k+1]==stationNum and timeInstance[k+1]>0)):
        upOpt=k+1
        break
upOpt0=upOpt
#determine the first service in upstream direction

upxk0=np.ones(stationNum)
upxk1=np.ones(stationNum)

downxk0=np.ones(stationNum)
downxk1=np.ones(stationNum)
#correspond to constraints x^{k}_{i}+x^{k+1}_{i}+x^{k+2}_{i}>=2

PreX=np.ones((serviceNum,2*stationNum))
PretimetableArr=np.zeros((serviceNum,2*stationNum))
PretimetableDep=np.zeros((serviceNum,2*stationNum))
#timetable with the shortest travel time

#-----------------------------------------------------------------------------------------------------------------------
#------------------------------------results obtained by the optimization solver----------------------------------------
FinalX=np.ones((serviceNum,2*stationNum))
FinaltimetableArr=np.zeros((serviceNum,2*stationNum))
FinaltimetableDep=np.zeros((serviceNum,2*stationNum))# timetable result

Finalw2=np.zeros([serviceNum,2*stationNum])
Finalw3=np.zeros([serviceNum,2*stationNum,2*stationNum])
Finalwr2=np.zeros([serviceNum,2*stationNum])
Finalwr3=np.zeros([serviceNum,2*stationNum,2*stationNum])
Finaln2=np.zeros([serviceNum,2*stationNum])
Finalnb3=np.zeros([serviceNum,2*stationNum,2*stationNum])
Finalnb2=np.zeros([serviceNum,2*stationNum])
Finalwb3=np.zeros([serviceNum,2*stationNum,2*stationNum])
Finalwb2=np.zeros([serviceNum,2*stationNum])# passenger flow obtained by Gurobi
#-----------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------------------------------------
#-----------------------------------results obtained by the simulation--------------------------------------------------
StandardtimetableDep=np.zeros((serviceNum,2*stationNum))
Standardw3=np.zeros((serviceNum,2*stationNum,2*stationNum))
StandardTotal=np.zeros((2,serviceNum))

SFinalw2=np.zeros([serviceNum,2*stationNum])
SFinalw3=np.zeros([serviceNum,2*stationNum,2*stationNum])
SFinalwr2=np.zeros([serviceNum,2*stationNum])
SFinalwr3=np.zeros([serviceNum,2*stationNum,2*stationNum])
SFinaln2=np.zeros([serviceNum,2*stationNum])
SFinalna2=np.zeros([serviceNum,2*stationNum])
SFinalnb3=np.zeros([serviceNum,2*stationNum,2*stationNum])
SFinalnb2=np.zeros([serviceNum,2*stationNum])
SFinalwb3=np.zeros([serviceNum,2*stationNum,2*stationNum])
SFinalwb2=np.zeros([serviceNum,2*stationNum])
# passenger flow obtained by the simulation
#-----------------------------------------------------------------------------------------------------------------------



while True:

#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------upstream----------------------------------------------------
#-------------------------------------------------------------------------------------------------------------

    if upOpt<serviceNum:#upstream

        k=upOpt

#----------------------------------------upstream existing service--------------------------------------------
#-------------------------------------------------------------------------------------------------------------

        if k<existingService and stationNum>int(startStatus[k]):

#---------------------------------calculate the shortest travel time for existing service k-------------------
            m.reset()
            x = {}
            for i in range(2*stationNum):
                x[i] = m.addVar(vtype=GRB.BINARY)

            timetableArr = {}
            for i in range(2*stationNum):
                timetableArr[i] = m.addVar(lb=-40)

            timetableDep = {}
            for i in range(2*stationNum):
                timetableDep[i] = m.addVar(lb=-40)

            m.update()

            m.addConstr(timetableArr[int(startStatus[k])]==START+timeInstance[k])#starting time constraints

            m.addConstr(x[int(startStatus[k])]==1)#the first station can not be skipped
            if k>upOpt0:
                for i in range(int(startStatus[k-1]),stationNum):
                    m.addConstr(timetableArr[i]-FinaltimetableDep[k-1,i]>=headway, name='headway_%s_%s'% (k,i))
                    #headway constraints
            if stationNum>int(startStatus[k])>=0.5*stationNum:
                m.addConstr(quicksum(x[i] for i in range(int(startStatus[k]),stationNum))>=stationNum-
                int(startStatus[k])-SK1)#skipping constraints
            elif 0.5*stationNum>int(startStatus[k]):
                m.addConstr(quicksum(x[i] for i in range(int(startStatus[k]),stationNum))>=stationNum-
                int(startStatus[k])-SK2)#skipping constraints

            for i in range(int(startStatus[k]),stationNum):
                m.addConstr(Maxdt>=timetableDep[i]-timetableArr[i])
                m.addConstr(Maxdt*x[i]>=timetableDep[i]-timetableArr[i])
                m.addConstr(timetableDep[i]-timetableArr[i]>=stationdt[i]*x[i])#dwelling time constraints

            for i in range(int(startStatus[k])+1,stationNum):
                m.addConstr(timetableArr[i]==timetableDep[i-1]+linkrun10[i-1]-(1-x[i-1])*accExtra-(1-x[i])*decExtra)#running time constraints

            for i in range(int(startStatus[k]),stationNum):
                if upxk0[i]==0 or upxk1[i]==0:
                    m.addConstr(x[i]==1)
            m.update()


            obj = LinExpr()


            obj.add(timetableDep[stationNum-1]-timetableArr[int(startStatus[k])])
            m.setObjective(obj, GRB.MINIMIZE)
            m.params.MIPGap=0.001
            m.optimize()

            T=obj.getValue()

            for i in range(int(startStatus[k]),stationNum):
                PreX[k,i]=x[i].X
                PretimetableArr[k,i]=timetableArr[i].X
                PretimetableDep[k,i]=timetableDep[i].X

#--------------------------------------------------------------------------------------------------------------------------------------------------------
#------------------------------------------the main programming for upstream existing services-----------------------------------------------------------
            m.reset()

            x = {}
            for i in range(2*stationNum):
                x[i] = m.addVar(vtype=GRB.BINARY)


            timetableArr = {}
            for i in range(2*stationNum):
                timetableArr[i] = m.addVar(lb=-40)

            timetableDep = {}
            for i in range(2*stationNum):
                timetableDep[i] = m.addVar(lb=-40)

            z={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    z[i,j] = m.addVar(vtype=GRB.BINARY, name='z_%s_%s_%s' % (k, i,j))

            w2 = {}
            for i in range(2*stationNum):
                w2[i] = m.addVar(lb=0)

            w3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    w3[i,j] = m.addVar(lb=0)


            wb2 = {}
            for i in range(2*stationNum):
                wb2[i] = m.addVar(lb=0)

            wb3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    wb3[i,j] = m.addVar(lb=0)

            wr3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    wr3[i,j] = m.addVar(lb=0)

            n2 = {}
            for i in range(2*stationNum):
                n2[i] = m.addVar(lb=0,ub=Capacity)

            na2 = {}
            for i in range(2*stationNum):
                na2[i] = m.addVar(lb=0,ub=Capacity)

            nb2 = {}
            for i in range(2*stationNum):
                nb2[i] = m.addVar(lb=0,ub=Capacity)

            nbx2={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    nbx2[i,j] = m.addVar(lb=0,ub=Capacity)

            nb3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    nb3[i,j] = m.addVar(lb=0,ub=Capacity)

            indicator1={}
            for i in range(2*stationNum):
                indicator1[i] = m.addVar(vtype=GRB.BINARY)

            indicator2={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    indicator2[i,j] = m.addVar(vtype=GRB.BINARY)

            m.update()

#----------------------------------------------timetable part----------------------------------
            m.addConstr(timetableArr[int(startStatus[k])]==START+timeInstance[k])#starting time constraints

            m.addConstr(x[int(startStatus[k])]==1)#the first station can not be skipped
            if k>upOpt0:
                for i in range(int(startStatus[k-1]),stationNum):
                    m.addConstr(timetableArr[i]-FinaltimetableDep[k-1,i]>=headway, name='headway_%s_%s'% (k,i))
                    #headway constraints
            if stationNum>int(startStatus[k])>=0.5*stationNum:
                m.addConstr(quicksum(x[i] for i in range(int(startStatus[k]),stationNum))>=stationNum-
                int(startStatus[k])-SK1)#skipping constraints
                PsnInTraink=np.zeros(2*N)
                for i in range(2*N):
                    PsnInTraink[i]=PsnInTrain0[k,i]#initialize passenger number in existing service
            elif 0.5*stationNum>int(startStatus[k]):
                m.addConstr(quicksum(x[i] for i in range(int(startStatus[k]),stationNum))>=stationNum-
                int(startStatus[k])-SK2)#skipping constraints
                PsnInTraink=np.zeros(2*N)
                if int(startStatus[k])>0:
                    for i in range(2*N):
                        PsnInTraink[i]=PsnInTrain0[k,i]#initialize passenger number in existing service

            for i in range(int(startStatus[k]),stationNum):
                m.addConstr(Maxdt>=timetableDep[i]-timetableArr[i])
                m.addConstr(Maxdt*x[i]>=timetableDep[i]-timetableArr[i])
                m.addConstr(timetableDep[i]-timetableArr[i]>=stationdt[i]*x[i])#dwelling time constraints

            for i in range(int(startStatus[k])+1,stationNum):
                m.addConstr(timetableArr[i]==timetableDep[i-1]+linkrun10[i-1]-(1-x[i-1])*accExtra-(1-x[i])*decExtra)#running time constraints

            for i in range(int(startStatus[k]),stationNum):
                if upxk0[i]==0 or upxk1[i]==0:
                    m.addConstr(x[i]==1)
                if i==N/2:
                    m.addConstr(timetableDep[i]-PretimetableDep[k,i]<=Delta1)
                if i==N-1:
                    m.addConstr(timetableDep[i]-PretimetableDep[k,i]<=Delta2)


            StandardtimetableDep[k,int(startStatus[k])]=START+timeInstance[k]+stationdt[int(startStatus[k])]
            for i in range(int(startStatus[k])+1,stationNum):
                StandardtimetableDep[k,i]+=StandardtimetableDep[k,i-1]+linkrun10[i-1]+stationdt[i]#standard departure time

#----------------------------------------------Passenger flow part----------------------------------

            for i in range(int(startStatus[k]),stationNum-1):
                for j in range(i+1,stationNum):
                    m.addConstr(wb3[i,j]<=w3[i,j])#no. of psgers from i to j
                    m.addConstr(wb3[i,j]<=M*x[i])
                    m.addConstr(wb3[i,j]<=M*x[j])
                    m.addConstr(wb3[i,j]>=w3[i,j]-M*(2-x[i]-x[j]))


            for i in range(int(startStatus[k]),stationNum-1):
                m.addConstr(quicksum(w3[i,j] for j in range(i+1,stationNum))==w2[i])#number of passengers on platform
                m.addConstr(quicksum(wb3[i,j] for j in range(i+1,stationNum))==wb2[i])#number of psgers wanting to get on
                m.addConstr(quicksum(nb3[i,j] for j in range(i+1,stationNum))==nb2[i])#no. of psgers can get on

            remainPsn=0
            temp=int(START)/W+H2

            if stationNum-1>int(startStatus[k])>=0.5*stationNum:

                i=int(startStatus[k])
                m.addConstr(na2[i]==PsnInTraink[i]+(1-x[i+1])*PsnInTraink[i+1])
                for j in range(i+1,stationNum):
                    Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]
                    m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[i,temp]*(timetableDep[i]-START)*OD2[i,j])
                    if j==i+1:
                        m.addConstr(wr3[i,j]==w3[i,j]+(1-x[j])*PsnInTraink[j]-nb3[i,j])#some passenger get off at int(startStatus[k])
                    else:
                        m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])#some passenger get off at int(startStatus[k])

                for i in range(int(startStatus[k])+1,stationNum-1):
                    m.addConstr(na2[i]==x[i]*PsnInTraink[i]+z[i,i+1]*PsnInTraink[i+1]+quicksum(nb3[j,i] for j in range(int(startStatus[k]),i)))
                    m.addConstr(2*z[i,i+1]>=x[i]-x[i+1]-0.5)
                    m.addConstr(2*z[i,i+1]<=x[i]-x[i+1]+1.5)
                    if k==upOpt0 or int(startStatus[k-1])>i:
                        temp=START/W+H2
                        for j in range(i+1,stationNum):
                            Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]
                            m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[i,temp]*(timetableDep[i]-START)*OD2[i,j])
                            if j==i+1:
                                m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#number of passengers on platform
                            else:
                                m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])#number of passengers on platform
                    else:
                        temp=FinaltimetableDep[k-1,i]/W+H2
                        temp=int(temp)
                        if temp>12:
                            temp=12
                        for j in range(i+1,stationNum):
                            Standardw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[i,temp]*(StandardtimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]
                            m.addConstr(w3[i,j]==SFinalwr3[k-1,i,j]+PsnVol0[i,temp]*(timetableDep[i]-FinaltimetableDep[k-1,i])*OD2[i,j])
                            if j==i+1:
                                m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#number of passengers on platform
                            else:
                                m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])
            elif 0.5*stationNum>int(startStatus[k]):
                i=int(startStatus[k])
                m.addConstr(na2[i]==PsnInTraink[i]+(z[i,i+1])*PsnInTraink[i+1]+(z[i,i+2])*PsnInTraink[i+2])
                m.addConstr(2*z[i,i+1]>=x[i]-x[i+1]-0.5)
                m.addConstr(2*z[i,i+1]<=x[i]-x[i+1]+1.5)
                m.addConstr(3*z[i,i+2]>=x[i]-x[i+1]-x[i+2]-0.5)
                m.addConstr(3*z[i,i+2]<=x[i]-x[i+1]-x[i+2]+2.5)
                for j in range(i+1,stationNum):
                    Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]
                    m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[i,temp]*(timetableDep[i]-START)*OD2[i,j])

                    if j<=i+2:
                        m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#some passenger get off at int(startStatus[k])
                    else:
                        m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])
                for i in range(int(startStatus[k])+1,stationNum-1):
                    if i<stationNum-2:
                        m.addConstr(na2[i]==x[i]*PsnInTraink[i]+(z[i,i+1])*PsnInTraink[i+1]+(z[i,i+2])*PsnInTraink[i+2]
                                                +quicksum(nb3[j,i] for j in range(int(startStatus[k]),i)))
                        m.addConstr(2*z[i,i+1]>=x[i]-x[i+1]-0.5)
                        m.addConstr(2*z[i,i+1]<=x[i]-x[i+1]+1.5)
                        m.addConstr(3*z[i,i+2]>=x[i]-x[i+1]-x[i+2]-0.5)
                        m.addConstr(3*z[i,i+2]<=x[i]-x[i+1]-x[i+2]+2.5)
                    elif i==stationNum-2:
                        m.addConstr(na2[i]==x[i]*PsnInTraink[i]+(z[i,i+1])*PsnInTraink[i+1]
                                                +quicksum(nb3[j,i] for j in range(int(startStatus[k]),i)))
                        m.addConstr(2*z[i,i+1]>=x[i]-x[i+1]-0.5)
                        m.addConstr(2*z[i,i+1]<=x[i]-x[i+1]+1.5)
                    if k==upOpt0 or int(startStatus[k-1])>i:
                        temp=START/W+H2
                        for j in range(i+1,stationNum):
                            Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]
                            m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[i,temp]*(timetableDep[i]-START)*OD2[i,j])
                            if j<=i+2:
                                m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#number of passengers on platform
                            else:
                                m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])#number of passengers on platform
                    else:
                        temp=FinaltimetableDep[k-1,i]/W+H2
                        temp=int(temp)
                        if temp>12:
                            temp=12
                        for j in range(i+1,stationNum):
                            Standardw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[i,temp]*(StandardtimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]
                            m.addConstr(w3[i,j]==SFinalwr3[k-1,i,j]+PsnVol0[i,temp]*(timetableDep[i]-FinaltimetableDep[k-1,i])*OD2[i,j])
                            if j<=i+2:
                                m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#number of passengers on platform
                            else:
                                m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])

            i=int(startStatus[k])
            m.addConstr(n2[i]==sum(PsnInTraink)-na2[i]+nb2[i])
            m.addConstr(nb2[i]<=wb2[i])
            m.addConstr(nb2[i]<=Capacity-sum(PsnInTraink)+na2[i])
            m.addConstr(nb2[i]>=wb2[i]-M2*(1-indicator1[i]))
            m.addConstr(nb2[i]>=Capacity-sum(PsnInTraink)+na2[i]-M2*(indicator1[i]))
            m.addConstr(indicator1[i]<=1+(-wb2[i]+Capacity-sum(PsnInTraink)+na2[i])/M1)
            m.addConstr(indicator1[i]>=(-wb2[i]+Capacity-sum(PsnInTraink)+na2[i])/M1)

            for i in range(int(startStatus[k])+1,stationNum-1):
                m.addConstr(n2[i]==n2[i-1]-na2[i]+nb2[i])
                m.addConstr(nb2[i]<=wb2[i])
                m.addConstr(nb2[i]<=Capacity-n2[i-1]+na2[i])
                m.addConstr(nb2[i]>=wb2[i]-M2*(1-indicator1[i]))
                m.addConstr(nb2[i]>=Capacity-n2[i-1]+na2[i]-M2*(indicator1[i]))
                m.addConstr(indicator1[i]<=1+(-wb2[i]+Capacity-n2[i-1]+na2[i])/M1)
                m.addConstr(indicator1[i]>=(-wb2[i]+Capacity-n2[i-1]+na2[i])/M1)

                # no. of psnger in train when departing from station i

            for i in range(int(startStatus[k]),stationNum-1):
                temp=sum(Standardw3[k,i])
                for j in range(i+1,stationNum):
                    m.addConstr(nb3[i,j]<=wb3[i,j])
                    m.addConstr(nb3[i,j]>=wb3[i,j]-M*indicator2[i,j])
                    m.addConstr(nb3[i,j]>=nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01))-M*(1-indicator2[i,j]))
                    m.addConstr(indicator2[i,j]>=(wb3[i,j]-nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01)))/M)
                    m.addConstr(indicator2[i,j]<=1+(wb3[i,j]-nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01)))/M)
                    m.addConstr(nb3[i,j]<=nbx2[i,j])
                    m.addConstr(nbx2[i,j]<=nb2[i])
                    m.addConstr(nbx2[i,j]<=M*x[j])
                    m.addConstr(nbx2[i,j]>=nb2[i]-M*(1-x[j]))

            m.update()

#------------------------------------------optimization part---------------------------------------------

            for i in range(int(startStatus[k]),stationNum-1):
                for j in range(i+1,stationNum):
                    StandardTotal[0,k]+=Standardw3[k,i,j]

            obj = LinExpr()

            for i in range(int(startStatus[k]),stationNum-1):
                for j in range(i+1,stationNum):
                    obj+=wr3[i,j]

            m.setObjective(obj, GRB.MINIMIZE)
            m.params.MIPGap=Gap

            startuptime=clock()
            m.optimize()
            enduptime=clock()

            for i in range(N):
                upxk0[i]=upxk1[i]
                upxk1[i]=1
            for i in range(int(startStatus[upOpt]),stationNum):
                FinaltimetableArr[upOpt,i]=int(timetableArr[i].X)
                FinaltimetableDep[upOpt,i]=int(timetableDep[i].X)
                FinalX[upOpt,i]=int(x[i].X)
                upxk1[i]=int(x[i].X)
                Finalw2[upOpt,i]=int(w2[i].X)
                Finaln2[upOpt,i]=int(n2[i].X)
                Finalnb2[upOpt,i]=int(nb2[i].X)
                Finalwb2[upOpt,i]=int(wb2[i].X)
                if i<stationNum-1:
                    for j in range(i+1,stationNum):
                        Finalw3[upOpt,i,j]=round(w3[i,j].X)
                        Finalwr3[upOpt,i,j]=round(wr3[i,j].X)
                        Finalwb3[upOpt,i,j]=round(wb3[i,j].X)
                        Finalnb3[upOpt,i,j]=round(nb3[i,j].X)

#-------------------------------------------------------------------------------------------------------------


        elif k<existingService and stationNum==int(startStatus[k]):
            FinaltimetableDep[k,stationNum-1]=START-timeInstance[k]
#-------------------------------------------------------------------------------------------------------------


#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------upstream future services------------------------------------------------

        elif k>=existingService:

#-----------------------------------calculate the shortest travel time for future service k-------------------
            m.reset()
            x = {}
            for i in range(2*stationNum):
                x[i] = m.addVar(vtype=GRB.BINARY)


            timetableArr = {}
            for i in range(2*stationNum):
                timetableArr[i] = m.addVar(lb=-40)

            timetableDep = {}
            for i in range(2*stationNum):
                timetableDep[i] = m.addVar(lb=-40)

            m.update()

            m.addConstr(quicksum(x[i] for i in range(stationNum))>=stationNum-SK3)#skipping constraints

            m.addConstr(timetableArr[0]-FinaltimetableDep[k-existingService,2*stationNum-1]>=backturning)#backturning constraints

            for i in range(stationNum):
                m.addConstr(timetableArr[i]-FinaltimetableDep[k-1,i]>=headway, name='headway_%s_%s'% (k,i))#headway

            for i in range(stationNum):
                m.addConstr(Maxdt*x[i]>=timetableDep[i]-timetableArr[i])
                m.addConstr(timetableDep[i]-timetableArr[i]>=stationdt[i]*x[i])#dwelling time constraints

            for i in range(1,stationNum):
                m.addConstr(timetableArr[i]==timetableDep[i-1]+linkrun10[i-1]-(1-x[i-1])*accExtra-(1-x[i])*decExtra)#running time constraints

            for i in range(stationNum):
                if upxk0[i]==0 or upxk1[i]==0:
                    m.addConstr(x[i]==1)

            obj = LinExpr()

            obj.add(timetableDep[stationNum-1]-FinaltimetableDep[k-serviceNum,2*stationNum-1])
            m.update()
            m.setObjective(obj, GRB.MINIMIZE)
            m.params.MIPGap=0.001
            startuptime=clock()
            m.optimize()
            enduptime=clock()
            print 'Pre timetable Computation Time: ',enduptime-startuptime

            T=obj.getValue()

            for i in range(stationNum):
                PreX[k,i]=x[i].X
                PretimetableArr[k,i]=timetableArr[i].X
                PretimetableDep[k,i]=timetableDep[i].X

#--------------------------------------------------------------------------------------------------------------------------------------------------------
#------------------------------------------the main programming for upstream future services-------------------------------------------------------------
            m.reset()
            x = {}
            for i in range(2*stationNum):
                x[i] = m.addVar(vtype=GRB.BINARY)


            timetableArr = {}
            for i in range(2*stationNum):
                timetableArr[i] = m.addVar(lb=-40)

            timetableDep = {}
            for i in range(2*stationNum):
                timetableDep[i] = m.addVar(lb=-40)

            z={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    z[i,j] = m.addVar(vtype=GRB.BINARY, name='z_%s_%s_%s' % (k, i,j))

            w2 = {}
            for i in range(2*stationNum):
                w2[i] = m.addVar(lb=0)

            w3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    w3[i,j] = m.addVar(lb=0)


            wb2 = {}
            for i in range(2*stationNum):
                wb2[i] = m.addVar(lb=0)

            wb3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    wb3[i,j] = m.addVar(lb=0)

            wr3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    wr3[i,j] = m.addVar(lb=0)

            n2 = {}
            for i in range(2*stationNum):
                n2[i] = m.addVar(lb=0,ub=Capacity)

            na2 = {}
            for i in range(2*stationNum):
                na2[i] = m.addVar(lb=0,ub=Capacity)

            nb2 = {}
            for i in range(2*stationNum):
                nb2[i] = m.addVar(lb=0,ub=Capacity)

            nbx2={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    nbx2[i,j] = m.addVar(lb=0,ub=Capacity)

            nb3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    nb3[i,j] = m.addVar(lb=0,ub=Capacity)

            indicator1={}
            for i in range(2*stationNum):
                indicator1[i] = m.addVar(vtype=GRB.BINARY)

            indicator2={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    indicator2[i,j] = m.addVar(vtype=GRB.BINARY)

            m.update()

#--------------------------------timetable part--------------------------------------
            m.addConstr(quicksum(x[i] for i in range(stationNum))>=stationNum-SK3)#skipping constraints

            m.addConstr(timetableArr[0]-FinaltimetableDep[k-existingService,2*stationNum-1]>=backturning)#backturning constraints


            for i in range(stationNum):
                m.addConstr(timetableArr[i]-FinaltimetableDep[k-1,i]>=headway, name='headway_%s_%s'% (k,i))#headway

            for i in range(stationNum):
                m.addConstr(Maxdt>=timetableDep[i]-timetableArr[i])
                m.addConstr(Maxdt*x[i]>=timetableDep[i]-timetableArr[i])
                m.addConstr(timetableDep[i]-timetableArr[i]>=stationdt[i]*x[i])#dwelling time constraints

            for i in range(1,stationNum):
                m.addConstr(timetableArr[i]==timetableDep[i-1]+linkrun10[i-1]-(1-x[i-1])*accExtra-(1-x[i])*decExtra)#running time constraints

            for i in range(stationNum):
                if upxk0[i]==0 or upxk1[i]==0:
                    m.addConstr(x[i]==1)
                if i==N/2:
                    m.addConstr(timetableDep[i]-PretimetableDep[k,i]<=Delta1)
                if i==N-1:
                    m.addConstr(timetableDep[i]-PretimetableDep[k,i]<=Delta2)

            StandardtimetableDep[k,0]=FinaltimetableDep[k-existingService,2*stationNum-1]+backturning
            for i in range(1,stationNum):
                StandardtimetableDep[k,i]+=StandardtimetableDep[k,i-1]+linkrun10[i-1]+stationdt[i]#standard departure time

#------------------------------------ passenger flow part--------------------------------
            for i in range(stationNum-1):
                for j in range(i+1,stationNum):
                    m.addConstr(wb3[i,j]<=w3[i,j])#no. of psgers from i to j
                    m.addConstr(wb3[i,j]<=M*x[i])
                    m.addConstr(wb3[i,j]<=M*x[j])
                    m.addConstr(wb3[i,j]>=w3[i,j]-M*(2-x[i]-x[j]))

            for i in range(stationNum-1):
                m.addConstr(quicksum(w3[i,j] for j in range(i+1,stationNum))==w2[i])#number of passengers on platform
                m.addConstr(quicksum(wb3[i,j] for j in range(i+1,stationNum))==wb2[i])#number of psgers wanting to get on
                m.addConstr(quicksum(nb3[i,j] for j in range(i+1,stationNum))==nb2[i])#no. of psgers can get on

            remainPsn=0
            for i in range(stationNum-1):
                if k==existingService and int(startStatus[k-1])>i:
                    temp=START/W+H2
                    for j in range(i+1,stationNum):
                        Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]

                        m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[i,temp]*(timetableDep[i]-START)*OD2[i,j])#number of passengers on platform
                else:
                    temp=FinaltimetableDep[k-1,i]/W+H2
                    temp=int(temp)
                    if temp>12:
                        temp=12
                    for j in range(i+1,stationNum):
                        Standardw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[i,temp]*(StandardtimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]

                        m.addConstr(w3[i,j]==SFinalwr3[k-1,i,j]+
                                                PsnVol0[i,temp]*(timetableDep[i]-FinaltimetableDep[k-1,i])*OD2[i,j])#number of passengers on platform

            for i in range(1,stationNum):
                m.addConstr(na2[i]==quicksum(nb3[j,i] for j in range(i)))
                #number of passenger getting off

            i=0
            m.addConstr(n2[i]==nb2[i])
            m.addConstr(nb2[i]<=wb2[i])
            m.addConstr(nb2[i]<=Capacity)
            m.addConstr(nb2[i]>=wb2[i]-M2*(1-indicator1[i]))
            m.addConstr(nb2[i]>=Capacity-M2*(indicator1[i]))
            m.addConstr(indicator1[i]<=1+(-wb2[i]+Capacity)/M1)
            m.addConstr(indicator1[i]>=(-wb2[i]+Capacity)/M1)

            for i in range(1,stationNum-1):
                m.addConstr(n2[i]==n2[i-1]-na2[i]+nb2[i])
                m.addConstr(nb2[i]<=wb2[i])
                m.addConstr(nb2[i]<=Capacity-n2[i-1]+na2[i])
                m.addConstr(nb2[i]>=wb2[i]-M2*(1-indicator1[i]))
                m.addConstr(nb2[i]>=Capacity-n2[i-1]+na2[i]-M2*(indicator1[i]))
                m.addConstr(indicator1[i]<=1+(-wb2[i]+Capacity-n2[i-1]+na2[i])/M1)
                m.addConstr(indicator1[i]>=(-wb2[i]+Capacity-n2[i-1]+na2[i])/M1)
                # no. of psnger in train when departing from station i

            for i in range(stationNum-1):
                temp=sum(Standardw3[k,i])
                for j in range(i+1,stationNum):
                    m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])# remaining passengers in platform
                    m.addConstr(nb3[i,j]<=wb3[i,j])
                    m.addConstr(nb3[i,j]>=wb3[i,j]-M*indicator2[i,j])
                    m.addConstr(nb3[i,j]>=nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01))-M*(1-indicator2[i,j]))
                    m.addConstr(indicator2[i,j]>=(wb3[i,j]-nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01)))/M)
                    m.addConstr(indicator2[i,j]<=1+(wb3[i,j]-nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01)))/M)
                    m.addConstr(nb3[i,j]<=nbx2[i,j])
                    m.addConstr(nbx2[i,j]<=nb2[i])
                    m.addConstr(nbx2[i,j]<=M*x[j])
                    m.addConstr(nbx2[i,j]>=nb2[i]-M*(1-x[j]))

#-----------------------------------------------optimization------------------------------------
            for i in range(stationNum-1):
                for j in range(i+1,stationNum):
                    StandardTotal[0,k]+=Standardw3[k,i,j]

            obj=LinExpr();
            for i in range(stationNum-1):
                for j in range(i+1,stationNum):
                    obj+=wr3[i,j]


            m.update()
            m.setObjective(obj, GRB.MINIMIZE)
            m.params.MIPGap=Gap

            startuptime=clock()
            m.optimize()
            enduptime=clock()

            for i in range(N):
                upxk0[i]=upxk1[i]
                upxk1[i]=1
            for i in range(stationNum):
                FinaltimetableArr[upOpt,i]=int(timetableArr[i].X)
                FinaltimetableDep[upOpt,i]=int(timetableDep[i].X)
                FinalX[upOpt,i]=int(x[i].X)
                upxk1[i]=int(x[i].X)
                Finalw2[upOpt,i]=int(w2[i].X)
                Finaln2[upOpt,i]=int(n2[i].X)
                Finalnb2[upOpt,i]=int(nb2[i].X)
                Finalwb2[upOpt,i]=int(wb2[i].X)
                if i <stationNum-1:
                    for j in range(i+1,stationNum):
                        Finalw3[upOpt,i,j]=round(w3[i,j].X)
                        Finalwr3[upOpt,i,j]=round(wr3[i,j].X)
                        Finalwb3[upOpt,i,j]=round(wb3[i,j].X)
                        Finalnb3[upOpt,i,j]=round(nb3[i,j].X)

#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#----------------------------Passenger flow simulation in the upstream direction------------------------------

    if upOpt<serviceNum:#upstream
        k=upOpt

#--------------------------------------------upstream existing services-------------------------------
        if k<existingService and stationNum>int(startStatus[k]):

            temp=int(START)/W+H2
            if stationNum>int(startStatus[k])>=0.5*stationNum:
                PsnInTraink=np.zeros(2*N)
                for i in range(2*N):
                    PsnInTraink[i]=PsnInTrain0[k,i]#initialize passenger number in existing service
            elif 0.5*stationNum>int(startStatus[k]):
                PsnInTraink=np.zeros(2*N)
                if int(startStatus[k])>0:
                    for i in range(2*N):
                        PsnInTraink[i]=PsnInTrain0[k,i]#initialize passenger number in existing service


            if stationNum-1>int(startStatus[k])>=0.5*stationNum:
                i=int(startStatus[k])
                SFinalna2[k,i]=PsnInTraink[i]+(1-FinalX[k,i+1])*PsnInTraink[i+1]#some passenger get off at int(startStatus[k])
                for j in range(i+1,stationNum):
                    SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]
                    SFinalw2[k,i]+=SFinalw3[k,i,j]
                    SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                    SFinalwb2[k,i]+=SFinalwb3[k,i,j]

                SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-sum(PsnInTraink)+SFinalna2[k,i])
                SFinaln2[k,i]=sum(PsnInTraink)-SFinalna2[k,i]+SFinalnb2[k,i]

                for j in range(i+1,stationNum):
                    SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                    if j==i+1:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+(1-FinalX[k,j])*PsnInTraink[j]
                    else:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]
                    SFinalwr2[k,i]+=SFinalwr3[k,i,j]

                for i in range(int(startStatus[k])+1,stationNum-1):
                    if k==upOpt0 or int(startStatus[k-1])>i:
                        temp=START/W+H2
                        for j in range(i+1,stationNum):
                            SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]#number of passengers on platform
                            SFinalw2[k,i]+=SFinalw3[k,i,j]
                            SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                            SFinalwb2[k,i]+=SFinalwb3[k,i,j]
                    else:
                        temp=FinaltimetableDep[k-1,i]/W+H2
                        temp=int(temp)
                        if temp>12:
                            temp=12
                        for j in range(i+1,stationNum):
                            SFinalw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[i,temp]*(FinaltimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]#number of passengers on platform
                            SFinalw2[k,i]+=SFinalw3[k,i,j]
                            SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                            SFinalwb2[k,i]+=SFinalwb3[k,i,j]

                    SFinalna2[k,i]=FinalX[k,i]*PsnInTraink[i]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[i+1]
                    for j in range(int(startStatus[k]),i):
                        SFinalna2[k,i]+=SFinalnb3[k,j,i]

                    SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-SFinaln2[k,i-1]+SFinalna2[k,i])
                    SFinaln2[k,i]=SFinaln2[k,i-1]-SFinalna2[k,i]+SFinalnb2[k,i]

                    for j in range(i+1,stationNum):
                        SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                        if j==i+1:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,j])*PsnInTraink[j]
                        else:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]
                        SFinalwr2[k,i]+=SFinalwr3[k,i,j]

            elif 0.5*stationNum>int(startStatus[k]):
                i=int(startStatus[k])
                SFinalna2[k,i]=PsnInTraink[i]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[i+1]+FinalX[k,i]*(1-FinalX[k,i+1])*(1-FinalX[k,i+2])*PsnInTraink[i+2]
                for j in range(i+1,stationNum):
                    SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]#some passenger get off at int(startStatus[k])
                    SFinalw2[k,i]+=SFinalw3[k,i,j]
                    SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                    SFinalwb2[k,i]+=SFinalwb3[k,i,j]

                SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-sum(PsnInTraink)+SFinalna2[k,i])
                SFinaln2[k,i]=sum(PsnInTraink)-SFinalna2[k,i]+SFinalnb2[k,i]
                for j in range(i+1,stationNum):
                    SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                    if j==i+1:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[j]
                    elif j==i+2:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,i+1])*(1-FinalX[k,i+2])*PsnInTraink[i+2]
                    else:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]
                    SFinalwr2[k,i]+=SFinalwr3[k,i,j]

                for i in range(int(startStatus[k])+1,stationNum-1):
                    if k==upOpt0 or int(startStatus[k-1])>i:
                        temp=START/W+H2
                        for j in range(i+1,stationNum):
                            SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]#some passenger get off at int(startStatus[k])
                            SFinalw2[k,i]+=SFinalw3[k,i,j]
                            SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                            SFinalwb2[k,i]+=SFinalwb3[k,i,j]
                    else:
                        temp=FinaltimetableDep[k-1,i]/W+H2
                        temp=int(temp)
                        if temp>12:
                            temp=12
                        for j in range(i+1,stationNum):
                            SFinalw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[i,temp]*(FinaltimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]#number of passengers on platform
                            SFinalw2[k,i]+=SFinalw3[k,i,j]
                            SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                            SFinalwb2[k,i]+=SFinalwb3[k,i,j]

                    if i<stationNum-2:
                        SFinalna2[k,i]=FinalX[k,i]*PsnInTraink[i]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[i+1]+FinalX[k,i]*(1-FinalX[k,i+1])*(1-FinalX[k,i+2])*PsnInTraink[i+2]
                    elif i==stationNum-2:
                        SFinalna2[k,i]=FinalX[k,i]*PsnInTraink[i]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[i+1]
                    for j in range(int(startStatus[k]),i):
                        SFinalna2[k,i]+=SFinalnb3[k,j,i]

                    SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-SFinaln2[k,i-1]+SFinalna2[k,i])
                    SFinaln2[k,i]=SFinaln2[k,i-1]-SFinalna2[k,i]+SFinalnb2[k,i]
                    for j in range(i+1,stationNum):
                        SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                        if j==i+1:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[j]
                        elif j==i+2:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,i+1])*(1-FinalX[k,i+2])*PsnInTraink[i+2]
                        else:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]
                        SFinalwr2[k,i]+=SFinalwr3[k,i,j]

#--------------------------------------------upstream future services-------------------------------
        elif k>=existingService:
            for i in range(stationNum-1):
                if k==existingService and int(startStatus[k-1])>i:
                    temp=START/W+H2
                    for j in range(i+1,stationNum):
                        SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]#number of passengers on platform
                        SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                        SFinalw2[k,i]+=SFinalw3[k,i,j]
                        SFinalwb2[k,i]+=SFinalwb3[k,i,j]
                else:
                    temp=FinaltimetableDep[k-1,i]/W+H2
                    temp=int(temp)
                    if temp>12:
                        temp=12
                    for j in range(i+1,stationNum):
                        SFinalw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[i,temp]*(FinaltimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]
                        SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                        SFinalw2[k,i]+=SFinalw3[k,i,j]
                        SFinalwb2[k,i]+=SFinalwb3[k,i,j]

            i=0
            SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity)
            SFinaln2[k,i]=SFinalnb2[k,i]
            for j in range(i+1,stationNum):
                SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]# remaining passengers in platform
                SFinalwr2[k,i]+=SFinalwr3[k,i,j]


            for i in range(1,stationNum-1):
                for j in range(i):
                    SFinalna2[k,i]+=SFinalnb3[k,j,i]
                SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-SFinaln2[k,i-1]+SFinalna2[k,i])
                SFinaln2[k,i]=SFinaln2[k,i-1]-SFinalna2[k,i]+SFinalnb2[k,i]
                for j in range(i+1,stationNum):
                    SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                    SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]# remaining passengers in platform
                    SFinalwr2[k,i]+=SFinalwr3[k,i,j]

#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
        print 'upOpt: ',upOpt
        upOpt=upOpt+1


#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#---------------------------------------#downstream-----------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------

    if downOpt<serviceNum:
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------

        k=downOpt

#----------------------------------------downstream existing service------------------------------------------
#-------------------------------------------------------------------------------------------------------------

        if k<upOpt0 and 0<int(startStatus[k]):#downstream service, before entering backturning section

#---------------------------------calculate the shortest travel time for existing service k-------------------
            m.reset()
            x = {}
            for i in range(2*stationNum):
                x[i] = m.addVar(vtype=GRB.BINARY)


            timetableArr = {}
            for i in range(2*stationNum):
                timetableArr[i] = m.addVar(lb=-40)

            timetableDep = {}
            for i in range(2*stationNum):
                timetableDep[i] = m.addVar(lb=-40)

            m.update()

            m.addConstr(timetableArr[int(startStatus[k])]==START+timeInstance[k])#starting time constraints


            m.addConstr(x[int(startStatus[k])]==1)#the first station can not be skipped
            if k>downOpt0:
                for i in range(int(startStatus[k-1]),2*stationNum):
                    m.addConstr(timetableArr[i]-FinaltimetableDep[k-1,i]>=headway, name='headway_%s_%s'% (k,i))
                    #headway constraints
            if 2*stationNum>int(startStatus[k])>=1.5*stationNum:
                m.addConstr(quicksum(x[i] for i in range(int(startStatus[k]),2*stationNum))>=2*stationNum-
                int(startStatus[k])-SK1)#skipping constraints
                PsnInTraink=np.zeros(2*N)
                for i in range(2*N):
                    PsnInTraink[i]=PsnInTrain[int(startStatus[k])-1,i]#initialize passenger number in existing service
            elif 1.5*stationNum>int(startStatus[k]):
                m.addConstr(quicksum(x[i] for i in range(int(startStatus[k]),2*stationNum))>=2*stationNum-
                int(startStatus[k])-SK2)#skipping constraints
                PsnInTraink=np.zeros(2*N)
                if int(startStatus[k])>stationNum:
                    for i in range(2*N):
                        PsnInTraink[i]=PsnInTrain[int(startStatus[k])-1,i]#initialize passenger number in existing service

            for i in range(int(startStatus[k]),2*stationNum):
                m.addConstr(Maxdt*x[i]>=timetableDep[i]-timetableArr[i])
                m.addConstr(timetableDep[i]-timetableArr[i]>=stationdt[i]*x[i])#dwelling time constraints

            for i in range(int(startStatus[k])+1,2*stationNum):
                m.addConstr(timetableArr[i]==timetableDep[i-1]+linkrun20[i-N-1]-(1-x[i-1])*accExtra-(1-x[i])*decExtra)#running time constraints

            for i in range(int(startStatus[k]),2*stationNum):
                if downxk0[i-N]==0 or downxk1[i-N]==0:
                    m.addConstr(x[i]==1)

            obj = LinExpr()

            obj.add(timetableDep[2*stationNum-1]-timetableArr[int(startStatus[k])])
            m.update()
            m.setObjective(obj, GRB.MINIMIZE)
            m.params.MIPGap=0.001
            m.optimize()

            T=obj.getValue()

            for i in range(int(startStatus[k]),2*stationNum):
                PreX[k,i]=x[i].X
                PretimetableArr[k,i]=timetableArr[i].X
                PretimetableDep[k,i]=timetableDep[i].X


#--------------------------------------------------------------------------------------------------------------------------------------------------------
#------------------------------------------the main programming for downstream existing services---------------------------------------------------------

            m.reset()
            x = {}
            for i in range(2*stationNum):
                x[i] = m.addVar(vtype=GRB.BINARY)


            timetableArr = {}
            for i in range(2*stationNum):
                timetableArr[i] = m.addVar(lb=-40)

            timetableDep = {}
            for i in range(2*stationNum):
                timetableDep[i] = m.addVar(lb=-40)

            z={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    z[i,j] = m.addVar(vtype=GRB.BINARY, name='z_%s_%s_%s' % (k, i,j))

            w2 = {}
            for i in range(2*stationNum):
                w2[i] = m.addVar(lb=0)

            w3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    w3[i,j] = m.addVar(lb=0)


            wb2 = {}
            for i in range(2*stationNum):
                wb2[i] = m.addVar(lb=0)

            wb3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    wb3[i,j] = m.addVar(lb=0)

            wr3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    wr3[i,j] = m.addVar(lb=0)

            n2 = {}
            for i in range(2*stationNum):
                n2[i] = m.addVar(lb=0,ub=Capacity)

            na2 = {}
            for i in range(2*stationNum):
                na2[i] = m.addVar(lb=0,ub=Capacity)

            nb2 = {}
            for i in range(2*stationNum):
                nb2[i] = m.addVar(lb=0,ub=Capacity)

            nbx2 = {}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    nbx2[i,j] = m.addVar(lb=0,ub=Capacity)

            nb3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    nb3[i,j] = m.addVar(lb=0,ub=Capacity)


            indicator1 = {}
            for i in range(2*stationNum):
                indicator1[i] = m.addVar(vtype=GRB.BINARY)

            indicator2={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    indicator2[i,j] = m.addVar(vtype=GRB.BINARY)


            m.update()

#----------------------------------------timetable part----------------------------------

            m.addConstr(timetableArr[int(startStatus[k])]==START+timeInstance[k])#starting time constraints

            StandardtimetableDep[k,int(startStatus[k])]=START+timeInstance[k]+stationdt[int(startStatus[k])]
            for i in range(int(startStatus[k])+1,2*stationNum):
                StandardtimetableDep[k,i]+=StandardtimetableDep[k,i-1]+linkrun20[i-N-1]+stationdt[i]#standard departure time

            m.addConstr(x[int(startStatus[k])]==1)#the first station can not be skipped
            if k>downOpt0:
                for i in range(int(startStatus[k-1]),2*stationNum):
                    m.addConstr(timetableArr[i]-FinaltimetableDep[k-1,i]>=headway, name='headway_%s_%s'% (k,i))
                    #headway constraints
            if 2*stationNum>int(startStatus[k])>=1.5*stationNum:
                m.addConstr(quicksum(x[i] for i in range(int(startStatus[k]),2*stationNum))>=2*stationNum-
                int(startStatus[k])-SK1)#skipping constraints
                PsnInTraink=np.zeros(2*N)
                for i in range(2*N):
                    PsnInTraink[i]=PsnInTrain0[k,i]#initialize passenger number in existing service
            elif 1.5*stationNum>int(startStatus[k]):
                m.addConstr(quicksum(x[i] for i in range(int(startStatus[k]),2*stationNum))>=2*stationNum-
                int(startStatus[k])-SK2)#skipping constraints
                PsnInTraink=np.zeros(2*N)
                if int(startStatus[k])>stationNum:
                    for i in range(2*N):
                        PsnInTraink[i]=PsnInTrain0[k,i]#initialize passenger number in existing service

            for i in range(int(startStatus[k]),2*stationNum):
                m.addConstr(Maxdt>=timetableDep[i]-timetableArr[i])
                m.addConstr(Maxdt*x[i]>=timetableDep[i]-timetableArr[i])
                m.addConstr(timetableDep[i]-timetableArr[i]>=stationdt[i]*x[i])#dwelling time constraints

            for i in range(int(startStatus[k])+1,2*stationNum):
                m.addConstr(timetableArr[i]==timetableDep[i-1]+linkrun20[i-N-1]-(1-x[i-1])*accExtra-(1-x[i])*decExtra)#running time constraints

            for i in range(int(startStatus[k]),2*stationNum):
                if downxk0[i-N]==0 or downxk1[i-N]==0:
                    m.addConstr(x[i]==1)
                if i==N+N/2:
                    m.addConstr(timetableDep[i]-PretimetableDep[k,i]<=Delta1)
                if i==2*N-1:
                    m.addConstr(timetableDep[i]-PretimetableDep[k,i]<=Delta2)

#----------------------------- passenger flow part-------------------------------------------

            for i in range(int(startStatus[k]),2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    m.addConstr(wb3[i,j]<=w3[i,j])#no. of psgers from i to j
                    m.addConstr(wb3[i,j]<=M*x[i])
                    m.addConstr(wb3[i,j]<=M*x[j])
                    m.addConstr(wb3[i,j]>=w3[i,j]-M*(2-x[i]-x[j]))


            for i in range(int(startStatus[k]),2*stationNum-1):
                m.addConstr(quicksum(w3[i,j] for j in range(i+1,2*stationNum))==w2[i])#number of passengers on platform
                m.addConstr(quicksum(wb3[i,j] for j in range(i+1,2*stationNum))==wb2[i])#number of psgers wanting to get on
                m.addConstr(quicksum(nb3[i,j] for j in range(i+1,2*stationNum))==nb2[i])#no. of psgers can get on

            temp=START/W+H2

            if 2*stationNum-1>int(startStatus[k])>=1.5*stationNum:
                i=int(startStatus[k])
                m.addConstr(na2[i]==PsnInTraink[i]+(1-x[i+1])*PsnInTraink[i+1])
                for j in range(i+1,2*stationNum):
                    Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]
                    m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(timetableDep[i]-START)*OD2[i,j])#some passenger get off at int(startStatus[k])

                    if j==i+1:
                        m.addConstr(wr3[i,j]==w3[i,j]+(1-x[j])*PsnInTraink[j]-nb3[i,j])#some passenger get off at int(startStatus[k])
                    else:
                        m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])#some passenger get off at int(startStatus[k])

                for i in range(int(startStatus[k])+1,2*stationNum-1):
                    m.addConstr(na2[i]==x[i]*PsnInTraink[i]+(z[i,i+1])*PsnInTraink[i+1]+quicksum(nb3[j,i] for j in range(int(startStatus[k]),i)))
                    m.addConstr(2*z[i,i+1]>=x[i]-x[i+1]-0.5)
                    m.addConstr(2*z[i,i+1]<=x[i]-x[i+1]+1.5)
                    if k==downOpt0 or int(startStatus[k-1])>i:
                        temp=START/W+H2
                        for j in range(i+1,2*stationNum):
                            Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]
                            m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(timetableDep[i]-START)*OD2[i,j])

                            if j==i+1:
                                m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#number of passengers on platform
                            else:
                                m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])#number of passengers on platform
                    else:
                        temp=FinaltimetableDep[k-1,i]/W+H2
                        temp=int(temp)
                        if temp>12:
                            temp=12
                        for j in range(i+1,2*stationNum):
                            Standardw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[2*N-1-i,temp]*(StandardtimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]
                            m.addConstr(w3[i,j]==SFinalwr3[k-1,i,j]+PsnVol0[2*N-1-i,temp]*(timetableDep[i]-FinaltimetableDep[k-1,i])*OD2[i,j])
                            if j==i+1:
                                m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#number of passengers on platform
                            else:
                                m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])
            elif 1.5*stationNum>int(startStatus[k]):
                i=int(startStatus[k])
                m.addConstr(na2[i]==PsnInTraink[i]+(z[i,i+1])*PsnInTraink[i+1]+(z[i,i+2])*PsnInTraink[i+2])
                m.addConstr(2*z[i,i+1]>=x[i]-x[i+1]-0.5)
                m.addConstr(2*z[i,i+1]<=x[i]-x[i+1]+1.5)
                m.addConstr(3*z[i,i+2]>=x[i]-x[i+1]-x[i+2]-0.5)
                m.addConstr(3*z[i,i+2]<=x[i]-x[i+1]-x[i+2]+2.5)
                for j in range(i+1,2*stationNum):
                    Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]
                    m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(timetableDep[i]-START)*OD2[i,j])#some passenger get off at int(startStatus[k])

                    if j<=i+2:
                        m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#some passenger get off at int(startStatus[k])
                    else:
                        m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])#some passenger get off at int(startStatus[k])
                for i in range(int(startStatus[k])+1,2*stationNum-1):
                    if i<2*stationNum-2:
                        m.addConstr(na2[i]==x[i]*PsnInTraink[i]+(z[i,i+1])*PsnInTraink[i+1]+(z[i,i+2])*PsnInTraink[i+2]
                                                +quicksum(nb3[j,i] for j in range(int(startStatus[k]),i)))
                        m.addConstr(2*z[i,i+1]>=x[i]-x[i+1]-0.5)
                        m.addConstr(2*z[i,i+1]<=x[i]-x[i+1]+1.5)
                        m.addConstr(3*z[i,i+2]>=x[i]-x[i+1]-x[i+2]-0.5)
                        m.addConstr(3*z[i,i+2]<=x[i]-x[i+1]-x[i+2]+2.5)
                    else:
                        m.addConstr(na2[i]==x[i]*PsnInTraink[i]+(z[i,i+1])*PsnInTraink[i+1]
                                                +quicksum(nb3[j,i] for j in range(int(startStatus[k]),i)))
                        m.addConstr(2*z[i,i+1]>=x[i]-x[i+1]-0.5)
                        m.addConstr(2*z[i,i+1]<=x[i]-x[i+1]+1.5)

                    if k==upOpt0 or int(startStatus[k-1])>i:
                        temp=START/W+H2
                        for j in range(i+1,2*stationNum):
                            Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]
                            m.addConstr(w3[i,j]==nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(timetableDep[i]-START)*OD2[i,j])#number of passengers on platform

                            if j<=i+2:
                                m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#number of passengers on platform
                            else:
                                m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])#number of passengers on platform
                    else:
                        temp=FinaltimetableDep[k-1,i]/W+H2
                        temp=int(temp)
                        if temp>12:
                            temp=12
                        for j in range(i+1,2*stationNum):
                            Standardw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[2*N-1-i,temp]*(StandardtimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]
                            m.addConstr(w3[i,j]==SFinalwr3[k-1,i,j]+PsnVol0[2*N-1-i,temp]*(timetableDep[i]-FinaltimetableDep[k-1,i])*OD2[i,j])#number of passengers on platform

                            if j<=i+2:
                                m.addConstr(wr3[i,j]==w3[i,j]+(z[i,j])*PsnInTraink[j]-nb3[i,j])#number of passengers on platform
                            else:
                                m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])

            i=int(startStatus[k])
            m.addConstr(n2[i]==sum(PsnInTraink)-na2[i]+nb2[i])
            m.addConstr(nb2[i]<=wb2[i])
            m.addConstr(nb2[i]<=Capacity-sum(PsnInTraink)+na2[i])
            m.addConstr(nb2[i]>=wb2[i]-M2*(1-indicator1[i]))
            m.addConstr(nb2[i]>=Capacity-sum(PsnInTraink)+na2[i]-M2*(indicator1[i]))
            m.addConstr(indicator1[i]<=1+(-wb2[i]+Capacity-sum(PsnInTraink)+na2[i])/M1)
            m.addConstr(indicator1[i]>=(-wb2[i]+Capacity-sum(PsnInTraink)+na2[i])/M1)

            for i in range(int(startStatus[k])+1,2*stationNum-1):
                m.addConstr(n2[i]==n2[i-1]-na2[i]+nb2[i])
                m.addConstr(nb2[i]<=wb2[i])
                m.addConstr(nb2[i]<=Capacity-n2[i-1]+na2[i])
                m.addConstr(nb2[i]>=wb2[i]-M2*(1-indicator1[i]))
                m.addConstr(nb2[i]>=Capacity-n2[i-1]+na2[i]-M2*(indicator1[i]))
                m.addConstr(indicator1[i]<=1+(-wb2[i]+Capacity-n2[i-1]+na2[i])/M1)
                m.addConstr(indicator1[i]>=(-wb2[i]+Capacity-n2[i-1]+na2[i])/M1)
                # no. of psnger in train when departing from station i

            for i in range(int(startStatus[k]),2*stationNum-1):
                temp=sum(Standardw3[k,i])
                for j in range(i+1,2*stationNum):
                    m.addConstr(nb3[i,j]<=wb3[i,j])
                    m.addConstr(nb3[i,j]>=wb3[i,j]-M*indicator2[i,j])
                    m.addConstr(nb3[i,j]>=nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01))-M*(1-indicator2[i,j]))
                    m.addConstr(indicator2[i,j]>=(wb3[i,j]-nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01)))/M)
                    m.addConstr(indicator2[i,j]<=1+(wb3[i,j]-nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01)))/M)
                    m.addConstr(nb3[i,j]<=nbx2[i,j])
                    m.addConstr(nbx2[i,j]<=nb2[i])
                    m.addConstr(nbx2[i,j]<=M*x[j])
                    m.addConstr(nbx2[i,j]>=nb2[i]-M*(1-x[j]))

#------------------------------------------------optimization part-------------------------------------

            for i in range(int(startStatus[k]),2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    StandardTotal[1,k]+=Standardw3[k,i,j]

            obj = LinExpr()

            for i in range(int(startStatus[k]),2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    obj+=wr3[i,j]

            m.update()
            m.setObjective(obj, GRB.MINIMIZE)
            m.params.MIPGap=Gap

            startuptime=clock()
            m.optimize()
            enduptime=clock()


            for i in range(N):
                downxk0[i]=downxk1[i]
                downxk1[i]=1

            for i in range(int(startStatus[downOpt]),2*stationNum):
                FinaltimetableArr[downOpt,i]=int(timetableArr[i].X)
                FinaltimetableDep[downOpt,i]=int(timetableDep[i].X)
                FinalX[downOpt,i]=int(x[i].X)
                downxk1[i-N]=int(x[i].X)
                Finalw2[downOpt,i]=int(w2[i].X)
                Finaln2[downOpt,i]=int(n2[i].X)
                Finalnb2[downOpt,i]=int(nb2[i].X)
                Finalwb2[downOpt,i]=int(wb2[i].X)
                if i<2*stationNum-1:
                    for j in range(i+1,2*stationNum):
                        Finalw3[downOpt,i,j]=round(w3[i,j].X)
                        Finalwr3[downOpt,i,j]=round(wr3[i,j].X)
                        Finalwb3[downOpt,i,j]=round(wb3[i,j].X)
                        Finalnb3[downOpt,i,j]=round(nb3[i,j].X)

#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------


        elif k<upOpt0 and 0==int(startStatus[k]):
            FinaltimetableDep[k,2*stationNum-1]=START-timeInstance[k]

#-------------------------------------------------------------------------------------------------------------

#-------------------------------------------------------------------------------------------------------------
#----------------------------------------downstream future services-------------------------------------------

        elif k>=upOpt0:
#------------------------------------calculate the shortest travel time of service k-------------------
            m.reset()
            x = {}
            for i in range(2*stationNum):
                x[i] = m.addVar(vtype=GRB.BINARY)


            timetableArr = {}
            for i in range(2*stationNum):
                timetableArr[i] = m.addVar(lb=-40)

            timetableDep = {}
            for i in range(2*stationNum):
                timetableDep[i] = m.addVar(lb=-40)

            m.update()


            m.addConstr(quicksum(x[i] for i in range(stationNum,2*stationNum))>=stationNum-SK3)#skipping constraints
            m.addConstr(timetableArr[stationNum]-FinaltimetableDep[k,stationNum-1]>=backturning)#backturning constraints


            for i in range(stationNum,2*stationNum):
                m.addConstr(timetableArr[i]-FinaltimetableDep[k-1,i]>=headway, name='headway_%s_%s'% (k,i))

            for i in range(stationNum,2*stationNum):
                m.addConstr(Maxdt*x[i]>=timetableDep[i]-timetableArr[i])
                m.addConstr(timetableDep[i]-timetableArr[i]>=stationdt[i]*x[i])#dwelling time constraints

            for i in range(stationNum+1,2*stationNum):
                m.addConstr(timetableArr[i]==timetableDep[i-1]+linkrun20[i-N-1]-(1-x[i-1])*accExtra-(1-x[i])*decExtra)#running time constraints

            for i in range(stationNum,2*stationNum):
                if downxk0[i-N]==0 or downxk1[i-N]==0:
                    m.addConstr(x[i]==1)

            obj = LinExpr()

            obj.add(timetableDep[2*stationNum-1]-FinaltimetableDep[k,stationNum-1])
            m.update()
            m.setObjective(obj, GRB.MINIMIZE)
            m.params.MIPGap=0.001
            m.optimize()

            T=obj.getValue()

            for i in range(stationNum,2*stationNum):
                PreX[k,i]=x[i].X
                PretimetableArr[k,i]=timetableArr[i].X
                PretimetableDep[k,i]=timetableDep[i].X


#--------------------------------------------------------------------------------------------------------------------------------------------------------
#--------------------------------------------the main programming for downstream future services---------------------------------------------------------


            m.reset()
            x = {}
            for i in range(2*stationNum):
                x[i] = m.addVar(vtype=GRB.BINARY)


            timetableArr = {}
            for i in range(2*stationNum):
                timetableArr[i] = m.addVar(lb=-40)

            timetableDep = {}
            for i in range(2*stationNum):
                timetableDep[i] = m.addVar(lb=-40)

            z={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    z[i,j] = m.addVar(vtype=GRB.BINARY, name='z_%s_%s_%s' % (k, i,j))

            w2 = {}
            for i in range(2*stationNum):
                w2[i] = m.addVar(lb=0)

            w3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    w3[i,j] = m.addVar(lb=0)


            wb2 = {}
            for i in range(2*stationNum):
                wb2[i] = m.addVar(lb=0)

            wb3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    wb3[i,j] = m.addVar(lb=0)

            wr3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    wr3[i,j] = m.addVar(lb=0)

            n2 = {}
            for i in range(2*stationNum):
                n2[i] = m.addVar(lb=0,ub=Capacity)

            na2 = {}
            for i in range(2*stationNum):
                na2[i] = m.addVar(lb=0,ub=Capacity)

            nb2 = {}
            for i in range(2*stationNum):
                nb2[i] = m.addVar(lb=0,ub=Capacity)


            nbx2 = {}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    nbx2[i,j] = m.addVar(lb=0,ub=Capacity)

            nb3={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    nb3[i,j] = m.addVar(lb=0,ub=Capacity)


            indicator1 = {}
            for i in range(2*stationNum):
                indicator1[i] = m.addVar(vtype=GRB.BINARY)

            indicator2={}
            for i in range(2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    indicator2[i,j] = m.addVar(vtype=GRB.BINARY)


            m.update()


#------------------------------------------ timetable part---------------------------
            m.addConstr(quicksum(x[i] for i in range(stationNum,2*stationNum))>=stationNum-SK3)#skipping constraints
            m.addConstr(timetableArr[stationNum]-FinaltimetableDep[k,stationNum-1]>=backturning)#backturning constraints


            StandardtimetableDep[k,stationNum]=FinaltimetableDep[k,stationNum-1]+backturning+stationdt[stationNum]
            for i in range(stationNum+1,2*stationNum):
                StandardtimetableDep[k,i]=StandardtimetableDep[k,i-1]+linkrun20[i-N-1]+stationdt[i]

            for i in range(stationNum,2*stationNum):
                m.addConstr(timetableArr[i]-FinaltimetableDep[k-1,i]>=headway, name='headway_%s_%s'% (k,i))

            for i in range(stationNum,2*stationNum):
                m.addConstr(Maxdt>=timetableDep[i]-timetableArr[i])
                m.addConstr(Maxdt*x[i]>=timetableDep[i]-timetableArr[i])
                m.addConstr(timetableDep[i]-timetableArr[i]>=stationdt[i]*x[i])#dwelling time constraints

            for i in range(stationNum+1,2*stationNum):
                m.addConstr(timetableArr[i]==timetableDep[i-1]+linkrun20[i-N-1]-(1-x[i-1])*accExtra-(1-x[i])*decExtra)#running time constraints

            for i in range(stationNum,2*stationNum):
                if downxk0[i-N]==0 or downxk1[i-N]==0:
                    m.addConstr(x[i]==1)
                if i==N+N/2:
                    m.addConstr(timetableDep[i]-PretimetableDep[k,i]<=Delta1)
                if i==2*N-1:
                    m.addConstr(timetableDep[i]-PretimetableDep[k,i]<=Delta2)


#------------------------------------------ passenger flow part---------------------------

            for i in range(stationNum,2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    m.addConstr(wb3[i,j]<=w3[i,j])#no. of psgers from i to j
                    m.addConstr(wb3[i,j]<=M*x[i])
                    m.addConstr(wb3[i,j]<=M*x[j])
                    m.addConstr(wb3[i,j]>=w3[i,j]-M*(2-x[i]-x[j]))

            for i in range(stationNum,2*stationNum-1):
                m.addConstr(quicksum(w3[i,j] for j in range(i+1,2*stationNum))==w2[i])#number of passengers on platform
                m.addConstr(quicksum(wb3[i,j] for j in range(i+1,2*stationNum))==wb2[i])#number of psgers wanting to get on
                m.addConstr(quicksum(nb3[i,j] for j in range(i+1,2*stationNum))==nb2[i])#no. of psgers can get on

            remainPsn=0
            for i in range(stationNum,2*stationNum-1):
                if k==upOpt0 and int(startStatus[k-1])>i:
                    temp=START/W+H2
                    for j in range(i+1,2*stationNum):
                        Standardw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(StandardtimetableDep[k,i]-START)*OD2[i,j]

                        m.addConstr(w3[i,j]==nWantOn2[i,j]+
                                                PsnVol0[2*N-1-i,temp]*(timetableDep[i]-START)*OD2[i,j])#number of passengers on platform
                else:
                    temp=FinaltimetableDep[k-1,i]/W+H2
                    temp=int(temp)
                    if temp>12:
                        temp=12
                    for j in range(i+1,2*stationNum):
                        Standardw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[2*N-1-i,temp]*(StandardtimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]

                        m.addConstr(w3[i,j]==SFinalwr3[k-1,i,j]+
                                                PsnVol0[2*N-1-i,temp]*(timetableDep[i]-FinaltimetableDep[k-1,i])*OD2[i,j])#number of passengers on platform

            for i in range(stationNum+1,2*stationNum):
                m.addConstr(na2[i]==quicksum(nb3[j,i] for j in range(stationNum,i)))
                #number of passenger getting off

            i=stationNum
            m.addConstr(n2[i]==nb2[i])
            m.addConstr(nb2[i]<=wb2[i])
            m.addConstr(nb2[i]<=Capacity)
            m.addConstr(nb2[i]>=wb2[i]-M2*(1-indicator1[i]))
            m.addConstr(nb2[i]>=Capacity-M2*(indicator1[i]))
            m.addConstr(indicator1[i]<=1+(-wb2[i]+Capacity)/M1)
            m.addConstr(indicator1[i]>=(-wb2[i]+Capacity)/M1)

            for i in range(stationNum+1,2*stationNum-1):
                m.addConstr(n2[i]==n2[i-1]-na2[i]+nb2[i])
                m.addConstr(nb2[i]<=wb2[i])
                m.addConstr(nb2[i]<=Capacity-n2[i-1]+na2[i])
                m.addConstr(nb2[i]>=wb2[i]-M2*(1-indicator1[i]))
                m.addConstr(nb2[i]>=Capacity-n2[i-1]+na2[i]-M2*(indicator1[i]))
                m.addConstr(indicator1[i]<=1+(-wb2[i]+Capacity-n2[i-1]+na2[i])/M1)
                m.addConstr(indicator1[i]>=(-wb2[i]+Capacity-n2[i-1]+na2[i])/M1)
                # no. of psnger in train when departing from station i

            for i in range(stationNum,2*stationNum-1):
                temp=sum(Standardw3[k,i])
                for j in range(i+1,2*stationNum):
                    m.addConstr(wr3[i,j]==w3[i,j]-nb3[i,j])# remaining passengers in platform
                    m.addConstr(nb3[i,j]<=wb3[i,j])
                    m.addConstr(nb3[i,j]>=wb3[i,j]-M*indicator2[i,j])
                    m.addConstr(nb3[i,j]>=nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01))-M*(1-indicator2[i,j]))
                    m.addConstr(indicator2[i,j]>=(wb3[i,j]-nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01)))/M)
                    m.addConstr(indicator2[i,j]<=1+(wb3[i,j]-nbx2[i,j]*(Standardw3[k,i,j]/float(temp+0.01)))/M)
                    m.addConstr(nb3[i,j]<=nbx2[i,j])
                    m.addConstr(nbx2[i,j]<=nb2[i])
                    m.addConstr(nbx2[i,j]<=M*x[j])
                    m.addConstr(nbx2[i,j]>=nb2[i]-M*(1-x[j]))

#-----------------------------------optimization part---------------------------------------

            obj = LinExpr()

            for i in range(stationNum,2*stationNum-1):
                for j in range(i+1,2*stationNum):
                    obj+=wr3[i,j]

            m.update()
            m.setObjective(obj, GRB.MINIMIZE)
            m.params.MIPGap=Gap

            startuptime=clock()
            m.optimize()
            enduptime=clock()

            for i in range(N):
                downxk0[i]=downxk1[i]
                downxk1[i]=1

            for i in range(stationNum,2*stationNum):
                FinaltimetableArr[downOpt,i]=int(timetableArr[i].X)
                FinaltimetableDep[downOpt,i]=int(timetableDep[i].X)
                FinalX[downOpt,i]=int(x[i].X)
                downxk1[i-N]=int(x[i].X)
                Finalw2[downOpt,i]=int(w2[i].X)
                Finaln2[downOpt,i]=int(n2[i].X)
                Finalnb2[downOpt,i]=int(nb2[i].X)
                Finalwb2[downOpt,i]=int(wb2[i].X)
                if i<2*stationNum-1:
                    for j in range(i+1,2*stationNum):
                        Finalw3[downOpt,i,j]=round(w3[i,j].X)
                        Finalwr3[downOpt,i,j]=round(wr3[i,j].X)
                        Finalwb3[downOpt,i,j]=round(wb3[i,j].X)
                        Finalnb3[downOpt,i,j]=round(nb3[i,j].X)

#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#----------------------------Passenger flow simulation in the downstream direction----------------------------
    if downOpt<serviceNum-10:#downstream
        temp=int(START)/W+H2
        k=downOpt

#--------------------------------------------downstream existing services-------------------------------

        if k<upOpt0 and 0<int(startStatus[k]):#
            if 2*stationNum>int(startStatus[k])>=1.5*stationNum:
                PsnInTraink=np.zeros(2*N)
                for i in range(2*N):
                    PsnInTraink[i]=PsnInTrain0[k,i]#initialize passenger number in existing service
            elif 1.5*stationNum>int(startStatus[k]):
                PsnInTraink=np.zeros(2*N)
                if int(startStatus[k])>stationNum:
                    for i in range(2*N):
                        PsnInTraink[i]=PsnInTrain0[k,i]#initialize passenger number in existing service


            if 2*stationNum-1>int(startStatus[k])>=1.5*stationNum:
                i=int(startStatus[k])
                SFinalna2[k,i]=FinalX[k,i]*PsnInTraink[i]+(1-FinalX[k,i+1])*PsnInTraink[i+1]#some passenger get off at int(startStatus[k])
                for j in range(i+1,2*stationNum):
                    SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]
                    SFinalw2[k,i]+=SFinalw3[k,i,j]
                    SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                    SFinalwb2[k,i]+=SFinalwb3[k,i,j]

                SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-sum(PsnInTraink)+SFinalna2[k,i])
                SFinaln2[k,i]=sum(PsnInTraink)-SFinalna2[k,i]+SFinalnb2[k,i]

                for j in range(i+1,2*stationNum):
                    SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                    if j==i+1:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+(1-FinalX[k,j])*PsnInTraink[j]
                    else:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]
                    SFinalwr2[k,i]+=SFinalwr3[k,i,j]

                for i in range(int(startStatus[k])+1,2*stationNum-1):
                    if k==upOpt0 or int(startStatus[k-1])>i:
                        temp=START/W+H2
                        for j in range(i+1,2*stationNum):
                            SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]#number of passengers on platform
                            SFinalw2[k,i]+=SFinalw3[k,i,j]
                            SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                            SFinalwb2[k,i]+=SFinalwb3[k,i,j]
                    else:
                        temp=FinaltimetableDep[k-1,i]/W+H2
                        temp=int(temp)
                        if temp>12:
                            temp=12
                        for j in range(i+1,2*stationNum):
                            SFinalw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[2*N-1-i,temp]*(FinaltimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]#number of passengers on platform
                            SFinalw2[k,i]+=SFinalw3[k,i,j]
                            SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                            SFinalwb2[k,i]+=SFinalwb3[k,i,j]

                        SFinalna2[k,i]=FinalX[k,i]*PsnInTraink[i]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[i+1]
                        for j in range(int(startStatus[k]),i):
                            SFinalna2[k,i]+=SFinalnb3[k,j,i]

                    SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-SFinaln2[k,i-1]+SFinalna2[k,i])
                    SFinaln2[k,i]=SFinaln2[k,i-1]-SFinalna2[k,i]+SFinalnb2[k,i]

                    for j in range(i+1,2*stationNum):
                        SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                        if j==i+1:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,j])*PsnInTraink[j]
                        else:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]
                        SFinalwr2[k,i]+=SFinalwr3[k,i,j]

            elif 1.5*stationNum>int(startStatus[k]):
                i=int(startStatus[k])
                SFinalna2[k,i]=PsnInTraink[i]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[i+1]+FinalX[k,i]*(1-FinalX[k,i+1])*(1-FinalX[k,i+2])*PsnInTraink[i+2]
                for j in range(i+1,2*stationNum):
                    SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]#some passenger get off at int(startStatus[k])
                    SFinalw2[k,i]+=SFinalw3[k,i,j]
                    SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                    SFinalwb2[k,i]+=SFinalwb3[k,i,j]

                SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-sum(PsnInTraink)+SFinalna2[k,i])
                SFinaln2[k,i]=sum(PsnInTraink)-SFinalna2[k,i]+SFinalnb2[k,i]
                for j in range(i+1,2*stationNum):
                    SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                    if j==i+1:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[j]
                    elif j==i+2:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,i+1])*(1-FinalX[k,i+2])*PsnInTraink[i+2]
                    else:
                        SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]
                    SFinalwr2[k,i]+=SFinalwr3[k,i,j]

                for i in range(int(startStatus[k])+1,2*stationNum-1):
                    if k==upOpt0 or int(startStatus[k-1])>i:
                        temp=START/W+H2
                        for j in range(i+1,2*stationNum):
                            SFinalw3[k,i,j]=nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(FinaltimetableDep[k,i]-START)*OD2[i,j]#some passenger get off at int(startStatus[k])
                            SFinalw2[k,i]+=SFinalw3[k,i,j]
                            SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                            SFinalwb2[k,i]+=SFinalwb3[k,i,j]
                    else:
                        temp=FinaltimetableDep[k-1,i]/W+H2
                        temp=int(temp)
                        if temp>12:
                            temp=12
                        for j in range(i+1,2*stationNum):
                            SFinalw3[k,i,j]=SFinalwr3[k-1,i,j]+PsnVol0[2*N-1-i,temp]*(FinaltimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]#number of passengers on platform
                            SFinalw2[k,i]+=SFinalw3[k,i,j]
                            SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                            SFinalwb2[k,i]+=SFinalwb3[k,i,j]

                    if i<stationNum-2:
                        SFinalna2[k,i]=FinalX[k,i]*PsnInTraink[i]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[i+1]+FinalX[k,i]*(1-FinalX[k,i+1])*(1-FinalX[k,i+2])*PsnInTraink[i+2]
                    elif i==stationNum-2:
                        SFinalna2[k,i]=FinalX[k,i]*PsnInTraink[i]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[i+1]
                    for j in range(int(startStatus[k]),i):
                        SFinalna2[k,i]+=SFinalnb3[k,j,i]


                    SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-SFinaln2[k,i-1]+SFinalna2[k,i])
                    SFinaln2[k,i]=SFinaln2[k,i-1]-SFinalna2[k,i]+SFinalnb2[k,i]
                    for j in range(i+1,2*stationNum):
                        SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                        if j==i+1:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,i+1])*PsnInTraink[j]
                        elif j==i+2:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]+FinalX[k,i]*(1-FinalX[k,i+1])*(1-FinalX[k,i+2])*PsnInTraink[i+2]
                        else:
                            SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]
                        SFinalwr2[k,i]+=SFinalwr3[k,i,j]

#--------------------------------------------downstream future services-------------------------------
        elif k>=upOpt0:

            for i in range(stationNum,2*stationNum-1):
                if k==upOpt0 and int(startStatus[k-1])>i:
                    temp=START/W+H2
                    for j in range(i+1,2*stationNum):
                        SFinalw3[k,i,j]==nWantOn2[i,j]+PsnVol0[2*N-1-i,temp]*(FinaltimetableDep[i]-START)*OD2[i,j]#number of passengers on platform
                        SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                        SFinalw2[k,i]+=SFinalw3[k,i,j]
                        SFinalwb2[k,i]+=SFinalwb3[k,i,j]
                else:
                    temp=FinaltimetableDep[k-1,i]/W+H2
                    temp=int(temp)
                    if temp>12:
                        temp=12
                    for j in range(i+1,2*stationNum):
                        SFinalw3[k,i,j]=Finalwr3[k-1,i,j]+PsnVol0[2*N-1-i,temp]*(FinaltimetableDep[k,i]-FinaltimetableDep[k-1,i])*OD2[i,j]
                        SFinalwb3[k,i,j]=FinalX[k,i]*FinalX[k,j]*SFinalw3[k,i,j]
                        SFinalw2[k,i]+=SFinalw3[k,i,j]
                        SFinalwb2[k,i]+=SFinalwb3[k,i,j]

            i=stationNum
            SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity)
            SFinaln2[k,i]=SFinalnb2[k,i]
            for j in range(i+1,2*stationNum):
                SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]# remaining passengers in platform


            for i in range(stationNum+1,2*stationNum-1):
                for j in range(i):
                    SFinalna2[k,i]+=SFinalnb3[k,j,i]
                SFinalnb2[k,i]=min(SFinalwb2[k,i],Capacity-SFinaln2[k,i-1]+SFinalna2[k,i])
                SFinaln2[k,i]=SFinaln2[k,i-1]-SFinalna2[k,i]+SFinalnb2[k,i]
                for j in range(i+1,2*stationNum):
                    SFinalnb3[k,i,j]=(float(SFinalnb2[k,i])/(SFinalwb2[k,i]+0.01))*SFinalwb3[k,i,j]
                    SFinalwr3[k,i,j]=SFinalw3[k,i,j]-SFinalnb3[k,i,j]# remaining passengers in platform
                    SFinalwr2[k,i]+=SFinalwr3[k,i,j]


        print 'downOpt: ',downOpt
        print 'is finished'
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------


    downOpt=downOpt+1
    if downOpt==serviceNum-10:
        break



#-------------------------------------------------------------------------------------------------------------
#---------------------------------------------------Plot figures-------------------------------------------------
#-------------------------------------------------------------------------------------------------------------

#------------------------------------Timetable only-----------------------------------
distance0=np.zeros((stationNum))

for i in range(1,stationNum):
    distance0[i]=distance0[i-1]+sectionlen[i-1]


fig=plt.figure()
ax=fig.add_subplot(111)

for i in range(0,stationNum):
    plt.plot([-50,xlength],[distance0[i],distance0[i]],'b:')

space1=np.zeros((2*stationNum))
time1=np.zeros((2*stationNum))


for k in range(existingService):
    if (int(startStatus[k])==0 and timeInstance[k]>0):
        continue;
    else:
        if int(startStatus[k])<stationNum:
            for i in range(int(startStatus[k]),stationNum):
                space1[2*i]=distance0[i]
                space1[2*i+1]=distance0[i]
                time1[2*i]=FinaltimetableArr[k,i]
                time1[2*i+1]=FinaltimetableDep[k,i]
            plt.plot(time1[2*int(startStatus[k]):2*stationNum],space1[2*int(startStatus[k]):2*stationNum],'b')
            for j in range(stationNum):
                space1[2*j]=distance0[stationNum-1-j]
                space1[2*j+1]=distance0[stationNum-1-j]
                time1[2*j]=FinaltimetableArr[k,j+stationNum]
                time1[2*j+1]=FinaltimetableDep[k,j+stationNum]
            plt.plot(time1,space1,'r')
            a=FinaltimetableDep[k,stationNum-1]
            b=FinaltimetableArr[k,stationNum]
            plt.plot([a,a,(a+b)/2.0,b,b],[distance0[stationNum-1],distance0[stationNum-1]+300,distance0[stationNum-1]+800,
            distance0[stationNum-1]+300,distance0[stationNum-1]],'g')
        else:
            for j in range(2*stationNum-int(startStatus[k])):
                space1[2*j]=distance0[2*stationNum-int(startStatus[k])-1-j]
                space1[2*j+1]=distance0[2*stationNum-int(startStatus[k])-1-j]
                time1[2*j]=FinaltimetableArr[k,j+int(startStatus[k])]
                time1[2*j+1]=FinaltimetableDep[k,j+int(startStatus[k])]
            plt.plot(time1[0:4*stationNum-2*int(startStatus[k])],space1[0:4*stationNum-2*int(startStatus[k])],'r')



for k in range(existingService,serviceNum):
    for i in range(stationNum):
        space1[2*i]=distance0[i]
        space1[2*i+1]=distance0[i]
        time1[2*i]=FinaltimetableArr[k,i]
        time1[2*i+1]=FinaltimetableDep[k,i]
    plt.plot(time1,space1,'b')
    a=FinaltimetableDep[k-existingService,2*stationNum-1]
    b=FinaltimetableArr[k,0]
    plt.plot([a,a,(a+b)/2.0,b,b],[0,-300,-800,-300,0],'g')

for k in range(existingService,serviceNum-10):
    for j in range(stationNum):
        space1[2*j]=distance0[stationNum-1-j]
        space1[2*j+1]=distance0[stationNum-1-j]
        time1[2*j]=FinaltimetableArr[k,j+stationNum]
        time1[2*j+1]=FinaltimetableDep[k,j+stationNum]
    plt.plot(time1,space1,'r')
    a=FinaltimetableDep[k,stationNum-1]
    b=FinaltimetableArr[k,stationNum]
    plt.plot([a,a,(a+b)/2.0,b,b],[distance0[stationNum-1],distance0[stationNum-1]+300,distance0[stationNum-1]+800,
    distance0[stationNum-1]+300,distance0[stationNum-1]],'g')

AfterTimetableArrUp=np.zeros([2,N]);
AfterTimetableDepUp=np.zeros([2,N]);
AfterTimetableArrDown=np.zeros([2,N]);
AfterTimetableDepDown=np.zeros([2,N]);

AfterTimetableArrUp[0,0]=FinaltimetableDep[serviceNum-1-19,2*N-1]+backturning
AfterTimetableDepUp[0,0]=AfterTimetableArrUp[0,0]+stationdt[0]
AfterTimetableArrUp[1,0]=FinaltimetableDep[serviceNum-1-18,2*N-1]+backturning+60
AfterTimetableDepUp[1,0]=AfterTimetableArrUp[1,0]+stationdt[0]

for i in range(1,N):
    AfterTimetableArrUp[0,i]=AfterTimetableDepUp[0,i-1]+linkrun10[i-1]
    AfterTimetableDepUp[0,i]=AfterTimetableArrUp[0,i]+stationdt[i]
    AfterTimetableArrUp[1,i]=AfterTimetableDepUp[1,i-1]+linkrun10[i-1]
    AfterTimetableDepUp[1,i]=AfterTimetableArrUp[1,i]+stationdt[i]


for i in range(stationNum):
    space1[2*i]=distance0[i]
    space1[2*i+1]=distance0[i]
    time1[2*i]=AfterTimetableArrUp[0,i]
    time1[2*i+1]=AfterTimetableDepUp[0,i]
plt.plot(time1,space1,'b--')
a=FinaltimetableDep[serviceNum-1-19,2*stationNum-1]
b=AfterTimetableArrUp[0,0]
plt.plot([a,a,(a+b)/2.0,b,b],[0,-300,-800,-300,0],'g')

for i in range(stationNum):
    space1[2*i]=distance0[i]
    space1[2*i+1]=distance0[i]
    time1[2*i]=AfterTimetableArrUp[1,i]
    time1[2*i+1]=AfterTimetableDepUp[1,i]
plt.plot(time1,space1,'b--')
a=FinaltimetableDep[serviceNum-1-18,2*stationNum-1]
b=AfterTimetableArrUp[1,0]
plt.plot([a,a,(a+b)/2.0,b,b],[0,-300,-800,-300,0],'g')

AfterTimetableArrDown[0,0]=FinaltimetableDep[serviceNum-1-9,N-1]+backturning
AfterTimetableDepDown[0,0]=AfterTimetableArrDown[0,0]+stationdt[N+0]
AfterTimetableArrDown[1,0]=FinaltimetableDep[serviceNum-1-8,N-1]+backturning+60
AfterTimetableDepDown[1,0]=AfterTimetableArrDown[1,0]+stationdt[N+0]
for i in range(1,N):
    AfterTimetableArrDown[0,i]=AfterTimetableDepDown[0,i-1]+linkrun20[i-1]
    AfterTimetableDepDown[0,i]=AfterTimetableArrDown[0,i]+stationdt[N+i]
    AfterTimetableArrDown[1,i]=AfterTimetableDepDown[1,i-1]+linkrun20[i-1]
    AfterTimetableDepDown[1,i]=AfterTimetableArrDown[1,i]+stationdt[N+i]

for i in range(stationNum):
    space1[2*i]=distance0[stationNum-1-i]
    space1[2*i+1]=distance0[stationNum-1-i]
    time1[2*i]=AfterTimetableArrDown[0,i]
    time1[2*i+1]=AfterTimetableDepDown[0,i]
plt.plot(time1,space1,'r--')
a=FinaltimetableDep[serviceNum-1-9,stationNum-1]
b=AfterTimetableArrDown[0,0]
plt.plot([a,a,(a+b)/2.0,b,b],[distance0[stationNum-1],distance0[stationNum-1]+300,distance0[stationNum-1]+800,distance0[stationNum-1]+300,distance0[stationNum-1]],'g')

for i in range(stationNum):
    space1[2*i]=distance0[stationNum-1-i]
    space1[2*i+1]=distance0[stationNum-1-i]
    time1[2*i]=AfterTimetableArrDown[1,i]
    time1[2*i+1]=AfterTimetableDepDown[1,i]
plt.plot(time1,space1,'r--')
a=FinaltimetableDep[serviceNum-1-8,stationNum-1]
b=AfterTimetableArrDown[1,0]
plt.plot([a,a,(a+b)/2.0,b,b],[distance0[stationNum-1],distance0[stationNum-1]+300,distance0[stationNum-1]+800,distance0[stationNum-1]+300,distance0[stationNum-1]],'g')


plt.xlim(-40,xlength)
plt.ylim(-1000,23000)

plt.xticks(range(0,xlength,600),range(0,xlength/60,10))
plt.yticks(distance0,range(1,1+stationNum))

plt.xlabel('time (min)')
plt.ylabel('station')


#------------------------------------Passenger flow evolution-----------------------------------

fig4,ax4=plt.subplots()


for i in range(0,stationNum):
    plt.plot([-40,xlength],[distance0[i],distance0[i]],'b:')

space1=np.zeros((2*stationNum))
time1=np.zeros((2*stationNum))

SFinalw1=np.zeros((2*stationNum))
SFinalwr1=np.zeros((2*stationNum))

for i in range(2*stationNum):
    for k in range(serviceNum):
        SFinalwr2[k,i]=SFinalwr3[k,i].sum()
        if SFinalwr2[k,i]>SFinalwr1[i]:
            SFinalwr1[i]=SFinalwr2[k,i]
        if SFinalw2[k,i]>SFinalw1[i]:
            SFinalw1[i]=SFinalw2[k,i]


for k in range(upOpt0,existingService):
    if (int(startStatus[k])==0 and timeInstance[k]>0):
        continue;
    else:
        for i in range(int(startStatus[k]),stationNum):
            space1[2*i]=distance0[i]
            space1[2*i+1]=distance0[i]
            time1[2*i]=FinaltimetableArr[k,i]
            time1[2*i+1]=FinaltimetableDep[k,i]
            if i<stationNum-1:
                if k<existingService-1 or (k==existingService-1 and existingService<serviceNum):
                    if i==int(startStatus[k]) or (k>0 and int(startStatus[k-1])>i):
                        xa=START
                        xb=FinaltimetableDep[k,i]
                        ya=distance0[i]+0.6*sectionlen[i]*(PsnAtStationStart[i])/(0.01+SFinalw1[i])
                        yb=distance0[i]+0.6*sectionlen[i]*SFinalw2[k,i]/(0.01+SFinalw1[i])
                        xab=np.linspace(xa,xb,num=100)
                        if ya==yb:
                            yab=ya*np.ones(100)
                        else:
                            yab=np.linspace(ya,yb,100)
                        zab=distance0[i]*np.ones(100)
                        ax4.fill_between(xab, zab, yab,color='brown')
                    xa=FinaltimetableDep[k,i]
                    xb=FinaltimetableDep[k+1,i]
                    ya=distance0[i]+0.6*sectionlen[i]*SFinalwr2[k,i]/(0.01+SFinalw1[i])
                    yb=distance0[i]+0.6*sectionlen[i]*SFinalw2[k+1,i]/(0.01+SFinalw1[i])
                    xab=np.linspace(xa,xb,num=100)
                    if ya==yb:
                        yab=ya*np.ones(100)
                    else:
                        yab=np.linspace(ya,yb,100)
                    zab=distance0[i]*np.ones(100)
                    ax4.fill_between(xab, zab, yab,color='brown')
        plt.plot(time1[2*int(startStatus[k]):2*stationNum],space1[2*int(startStatus[k]):2*stationNum],'b')


for k in range(existingService,serviceNum):
    for i in range(stationNum):
        space1[2*i]=distance0[i]
        space1[2*i+1]=distance0[i]
        time1[2*i]=FinaltimetableArr[k,i]
        time1[2*i+1]=FinaltimetableDep[k,i]
        if i<stationNum-1:
            if k==existingService and int(startStatus[k-1])>i:
                xa=START
                xb=FinaltimetableDep[k,i]
                ya=distance0[i]+0.6*sectionlen[i]*(PsnAtStationStart[i])/(0.01+SFinalw1[i])
                yb=distance0[i]+0.6*sectionlen[i]*SFinalw2[k,i]/(0.01+SFinalw1[i])
                xab=np.linspace(xa,xb,num=100)
                if ya==yb:
                    yab=ya*np.ones(100)
                else:
                    yab=np.linspace(ya,yb,100)
                zab=distance0[i]*np.ones(100)
                ax4.fill_between(xab, zab, yab,color='brown')
            if k<serviceNum-1:
                xa=FinaltimetableDep[k,i]
                xb=FinaltimetableDep[k+1,i]
                ya=distance0[i]+0.6*sectionlen[i]*SFinalwr2[k,i]/(0.01+SFinalw1[i])
                yb=distance0[i]+0.6*sectionlen[i]*SFinalw2[k+1,i]/(0.01+SFinalw1[i])
                xab=np.linspace(xa,xb,num=100)
                if ya==yb:
                    yab=ya*np.ones(100)
                else:
                    yab=np.linspace(ya,yb,100)
                zab=distance0[i]*np.ones(100)
                ax4.fill_between(xab, zab, yab,color='brown')
    plt.plot(time1,space1,'b')
    a=FinaltimetableDep[k-existingService,2*stationNum-1]
    b=FinaltimetableArr[k,0]


plt.xlim(-40,xlength)
plt.ylim(-1000,23000)

plt.xticks(range(0,xlength,600),range(0,xlength/60,10))
plt.yticks(distance0,range(1,1+stationNum))

plt.xlabel('time (min)')
plt.ylabel('station')

plt.show()

#-------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------