Problem with current time

Hey,

Currently i try to Model a Orienteering Problem with time windows in a combination with a Pickup and Delivery Problem. Instead of gaining Score by visiting vertices, the score is gained by pickups and deliveries of vehicles in the path. 

Orienteering Problem with time windows:

Pick -Up and Delivery: 

https://link.springer.com/article/10.1007/s00291-016-0454-y#Sec2

So far the Orienteering Problem Works good, same as PuD. The Problem i have is with my "timenow". It declares the time at the current vertex i. I have no idea, how i can declare it in a nother way. Here i just made it as a grb.var with an index in range n. 

my Code:

import random
import gurobipy as grb
import math
rnd = random

n = 15#Vertices
Distance = 25#Max Distance
npv = 4 #number of pickup vertices and number of delivery vertices 
P= [rnd.randint(1,n-1) for i in range(npv)] #set of pickup vertices
D = [rnd.randint(1,n-1) for i in range(npv)]#set of delivery vertices 
for i in range (npv): #start != ende 
    if (D[i] == P[i]):
        D[i] -= 1
# recv= []
# for i in range(len(P)):
#     recv.append(rnd.randint(10,20))
recv={i: rnd.randint(8,11) for i in D}

def distance(points, i, j):
    dx = points[i][0] - points[j][0]
    dy = points[i][1] - points[j][1]
    return math.sqrt(dx*dx + dy*dy)

random.seed(1)
points = []
for i in range(n):
    points.append((random.randint(0,25),random.randint(0,25)))
opt_model = grb.Model(name="Platooning Model")
times={i: rnd.randint(8,11) for i in range(n)}
timee={i: rnd.randint(14,21) for i in range(n)}
timestart = 8
M = 1000000000




# <= Variables
x_vars = {}
for i in range(n):
    for j in range(n):
        x_vars[i,j] = opt_model.addVar(vtype=grb.GRB.BINARY,name='e'+str(i)+'_'+str(j))
u={}
for i in range(1,n):
    u[i]=opt_model.addVar(vtype=grb.GRB.CONTINUOUS,
                          name='e'+str(i))

timenow = opt_model.addVars(n,vtype = grb.GRB.CONTINUOUS,name='s'+str(i) )
opt_model.addConstr(timenow[0] == timestart)
# print(P,D)
print(recv)
# print(times)
# print(timee)
# print(timenow)

# <= Constraint (Mandatory Edges and excluding vertexes) Eq(1) hier 1,n zu n
opt_model.addConstr((grb.quicksum(x_vars[0,j] for j in range(n))) == 1)
opt_model.addConstr((grb.quicksum(x_vars[i,0] for i in range(n))) == 1)
opt_model.addConstr((grb.quicksum(x_vars[i,i] for i in range(1,n))) == 0)
# <= Constraint (Distance) Eq(3) Hier statt n-1 -> n
for i in range(n):
opt_model.addConstr(grb.quicksum(x_vars[i,j]*distance(points, i, j) for j in range(1,n) ) <= Distance)
# <= Constraint (Equality & Single edge in and out) Eq(2)
for k in range(1, n):
opt_model.addConstr(grb.quicksum(x_vars[i,k] for i in range(n))== grb.quicksum(x_vars[k,j] for j in range(n)))
opt_model.addConstr(grb.quicksum(x_vars[i,k] for i in range(n))<=1) 
opt_model.addConstr(grb.quicksum(x_vars[k,j] for j in range(n))<=1)
#<= Constraint (Subtour elimination) Eq(4) Eq(5)
for i in range(1,n):
opt_model.addConstr(2 <= u[i])
opt_model.addConstr(u[i] <= n)
for i in range(1,n):
for j in range(1,n):
opt_model.addConstr((u[i] - u[j] +1 <= (n-1)*(1-x_vars[i,j])))
#objective = grb.quicksum(recv[i]*x_vars[P[i],j] for i in range(len(P)) for j in range(n))
#Objective nach Delivery 
objective = grb.quicksum(recv[j]*x_vars[i,j] for i in range(n) for j in D)

for i in range(len(P)):
opt_model.addConstr(( grb.quicksum(x_vars[P[i],j] for j in range(n)) - grb.quicksum(x_vars[D[i], j] for j in range(n))==0))
for i in range(n):
opt_model.addConstr(times[i] <= timenow[i]) 
opt_model.addConstr(timenow[i]<=timee[i])

opt_model.addConstrs( timenow[P[i]] <= timenow[D[i]] for i in range(len(P)))
#13 Pud
opt_model.addConstrs( timenow[j] >= x_vars[i,j] * ( timenow[i] + distance(points, i, j)) for j in range (1,n) for i in range (1,n))
#BIG M
#opt_model.addConstrs(timenow[i] + distance(points, i, j) - timenow[j]<= M * (1-x_vars[i,j]) for j in range (1,n) for i in range (1,n))

opt_model.ModelSense = grb.GRB.MAXIMIZE
opt_model.setObjective(objective)
# opt_model.update()
# opt_model.Params.MIPGap = 0.15 #stop 15% gap
# opt_model.Params.TimeLimit = 60 #seconds timelimit
opt_model.optimize()


import matplotlib.pyplot as plt
select = grb.tuplelist((i,j) for i,j in x_vars.keys() if x_vars[i,j].X > 0.5)
x_val = [x[0] for x in points]
y_val = [x[1] for x in points]
for i,j in select:
plt.annotate('$Standort_%d$'%(i),(x_val[i]+2,y_val[i],))#0-9 
plt.scatter(x_val[1:],y_val[1:], c='b') #locations for 1 till n+1 customer/clients
for i,j in select:
plt.plot([x_val[i],x_val[j]],[y_val[i],y_val[j]], c = 'g', zorder=0,alpha=0.3)
plt.plot(x_val[0], y_val[0], c = 'r', marker='s')
platoonsize = 0
for i,j in select:
for p in P:
if p == i:
platoonsize = platoonsize + 1
for d in D:
if d == j:
if j>i:
platoonsize = platoonsize - 1
print("Maximale Platoongröße = ",platoonsize)
plt.show()
print(select)
print(recv)
print(P,D)