Got Stuck in Constructing Routes in a periodic vehicle routing model

Hi, 

I am trying to construct a periodic vehicle routing problem with a slight difference, my model should decide on the visiting schedules of each node based on their required visiting frequencies. I have currently 40 nodes and a depot, 3 vehicles, and 40 days in my planning horizon. As I said, my model basically should do two things, i) it should satisfy the visiting times of each node (regarding their frequency of visit)  and İİ) a vehicle completes its tour (i.e, depart from depot (node 0), go to its assigned nodes for that period, and return to the depot at the end of the same period.) But each vehicle does not suppose to depart from depot at a given period as in the classical VRP, neither they should do that in every period. Everything should suppose to happen just when it is needed (and the visiting frequencies should determine this condition), and in the remaining times vehicles can remain idle in the depot, as this is what we essentially want when we try to minimize distance. similarly, not all the vehicles are supposed to be utilized if I can handle the tours with the same vehicle. I plan to solve this problem in a very large network consisting of 175 nodes, and I know when I do that I'm gonna need 3 vehicles or more. But since I cannot run my code in a simple pilot network, I cannot move forward. Here is my current work:

import random
import pandas as pd
import gurobipy as gp
from gurobipy import*
import matplotlib.pyplot as plt
mdl = Model("periodic vrp")
'''Sets and parameters'''

N = []
for i in range(31):
N.append(i)

I = []
for i in range(1,31):
I.append(i)

I2 = [3,9,12,15,16,21,23,27,28] #set of nodes to be visited twice in 40 days

I4 = [11,14] #set of nodes to be visited four times in 40 days

K = [1, 2, 3]

T = []
for i in range(1,41):
T.append(i)

T2 = []
for i in range(1,21):
T2.append(i)

C=5

random.seed(1)

def distance(p1,p2):
x_dist=(p1[0]-p2[0])
y_dist=(p1[1]-p2[1])
euclidian=(x_dist**2+y_dist**2)**0.5
return euclidian

D = [] # set of main distribution centers
D_loc = {}
for i in range(31):
D.append(i)
D_loc[i] = ((random.randint(0,100)),random.randint(0,100))

DIST = {}
for i in D:
for j in D: 
DIST[(i,j)] = distance(D_loc[i],D_loc[j])

q = {1:0.72,2:0.58,3:0.93,4:0.68,5:0.52,6:0.38,7:0.55,8:0.28,9:1.25,10:0.2,11:1.77,12:0.92,13:0.43,14:1.73,15:1.2,16:1.1,17:0.45,18:0.53,19:0.43,20:0.5,21:0.92,22:0.57,23:0.77,24:0.35,25:0.71,26:0.63,27:0.77,28:0.87,29:0.38,30:0.61}

FREQ = {1:1,2:1,3:2,4:1,5:1,6:1,7:1,8:1,9:2,10:1,11:4,12:2,13:1,14:4,15:2,16:2,17:1,18:1,19:1,20:1,21:2,22:1,23:2,24:1,25:1,26:1,27:2,28:2,29:1,30:1}


'''Decision Variables'''

pairs= [(n,m,k,t) for n in N for m in N for k in K for t in T]
x = mdl.addVars(pairs, vtype=GRB.BINARY, name ="x") #this is 1 if vehicle-k goes from node-n to node-m at period-t

bins=[(i,k,t) for i in I for k in K for t in T]
y = mdl.addVars(bins, vtype=GRB.BINARY, name ="y") #this is 1 if node-i is visited by vehicle-k at period-t

visits=[(i,t) for i in I for t in T]
v = mdl.addVars(visits, vtype=GRB.BINARY, name ="v") #this is 1 if node-i is visited at period-t

mtzs=[(i,k,t) for i in I for k in K for t in T]
u = mdl.addVars(mtzs, name ="u") #this is for the MTZ constraints of the capacitated VRP

'''Constraints'''

#----------flow balance constraints 
# for i in I:
# for k in K:
# for t in T:
# mdl.addConstr(quicksum(x[n,i,k,t] for n in N if n !=i) == quicksum(x[i,n,k,t] for n in N if n !=i))

for i in I:
for k in K:
for t in T:
mdl.addConstr(quicksum(x[n,i,k,t] for n in N if n !=i) == y[i,k,t])
for i in I:
for k in K:
for t in T:
if t !=40:
mdl.addConstr(quicksum(x[i,n,k,t+1] for n in N if n !=i) == y[i,k,t])


#----------connection of decision variables to each other
mdl.addConstrs(quicksum(y[i,k,t] for k in K) == v[i,t] for i in I for t in T)

for n in I:
for i in I:
for t in T:
if n !=i:
mdl.addConstr(quicksum(x[n,i,k,t] for k in K) <= (v[n,t]+v[i,t])/2)

for t in T:
for j in I:
mdl.addConstr(quicksum(x[n,i,k,t] for n in N for k in K) == v[i,t])

#------------
mdl.addConstrs(quicksum(y[i,k,t] for k in K for t in T) == FREQ[i] for i in I) #every node should be visited as much as its required frequency

mdl.addConstrs(quicksum(x[0,i,k,40] for i in I) == 0 for k in K)
mdl.addConstrs(quicksum(x[0,i,k,t] for i in I) <= 1 for t in T for k in K)

# mdl.addConstrs(quicksum(x[i,0,k,t] for i in I) == 1 for t in T for k in K)

mdl.addConstrs(u[j,k,t] - u[i,k,t] >= q[j]/FREQ[j]-C*(1-x[i,j,k,t]) for i in I for j in I for k in K for t in T if i!=j)
mdl.addConstrs(q[i] <= u[i,k,t] for i in I for k in K for t in T) 
mdl.addConstrs(u[i,k,t] <= C for i in I for k in K for t in T) 


start_time, end_time = T[0], T[-1] + 1
mdl.addConstrs(
(v[i, t] == 1)
>> (
quicksum(
v[i, j] for j in range(max(t + 1, start_time), min(t + 20, end_time))
) == 0
)
for i in I2
for t in T
)

start_time, end_time = T[0], T[-1] + 1
mdl.addConstrs(
(v[i, t] == 1)
>> (
quicksum(
v[i, j] for j in range(max(t + 1, start_time), min(t + 10, end_time))
) == 0
)
for i in I4
for t in T
)


'''Objective Function'''
obj=quicksum(x[n,m,k,t]*DIST[(n,m)] for n in N for m in N for k in K for t in T)
mdl.setObjective(obj, GRB.MINIMIZE)

'''Solve'''

mdl.optimize()

print('Obj: %g' % mdl.objVal)

for v in mdl.getVars():
if v.x==1:
print(v)

I am supposed to get complete tours in each period regarding the visiting frequencies, but this code makes partial visits to some nodes at different periods, although I wrote all my constraints with respect to the logic '\forall t \in T'. I hope I could explain myself well, and you can respond to this post in a short while, I have a similar post like this in the programming section of the community as well as an extension to my coding problem. You can reply to whichever you like. 

Thank you in advance for your help, I really appreciate your support.

Best regards, Benay