Python Gurobi not able to get optimize solution

Question

I'm working on a Gurobi optimize question but can't not find the solution I want. This question is trying to find the best path for picking materials in warehouse.

Here is the current solution.

x(0,5,33) 1.0
q(1,33) 2.0
x(1,2,33) 1.0
q(2,33) 2.0
x(2,3,33) 1.0
q(3,33) 2.0
x(3,4,33) 1.0
q(4,33) 2.0
x(4,1,33) 1.0
q(5,33) 2.0
x(5,0,33) 1.0

Here is my code.

from gurobipy import*

N = [0, 1, 2, 3, 4, 5] #Material location in warehosue, 0 is depot
S = [1] #The configuration need to be pick up

#Material list with location n
I, K = multidict({1:11, 2:12, 3:13, 4:42, 5:41})

#Distance between each location
R1 = {(0,1): 3, (0,2): 6, (0,3): 5, (0,4):4, (0,5): 2,\
      (1,2): 3, (1,3): 4, (1,4): 3, (1,5): 3,\
     (2,3): 3, (2,4): 6, (2,5): 4, (3,4): 3, (3,5): 7, (4,5): 6}
R = {}
for i,j in R1:
    R[j,i]=R1[i,j]
    R[i,j]=R1[i,j] #把0,1 -> 1,0

#Capacity
C = 10
#the batch of picking
L = [l for l in range(1,60)]

#BOM list (S:K1,K2,K3,K4,K5) 
BOM = {1:[11,12,13,42,41]}

#Quantity demand of configuration
S, D = multidict({1:2})

#Decision variable
model = Model("Picking path")
x = {}
q = {}
u = {}
for l in L:
    for i in range(len(N)+1):
        q[i,l]=model.addVar(vtype="C", name="q(%s,%s)" %(i,l))
        for j in range(len(N)+1):
            x[i,j,l]=model.addVar(vtype="B", name="x(%s,%s,%s)" %(i,j,l))
for l in L:
    for i in range(len(N)+1):
        u[i,l]=model.addVar(vtype="C", name="u(%s)" %i)

model.update()

#Constraints
BigM = 10**6
for l in L:
    model.addConstr(quicksum(x[0,j,l] for j in N if j!=0)*BigM >= quicksum(q[j,l] for j in N), name="A.3")
    model.addConstr(quicksum(x[i,0,l] for i in N if i!=0)*BigM >= quicksum(q[j,l] for j in N), name="A.4")
    for j in N:
        model.addConstr(quicksum(x[i,j,l] for i in N if i!=j)\
                       -quicksum(x[j,v,l] for v in N if v!=j)==0, name="A.2")
        model.addConstr(quicksum(x[i,j,l] for i in N if i!=j)*BigM >= q[j,l], name="A.5")
        model.addConstr(q[j,l] <= C, name="A.6")
            
for i in range(1,len(N)):
    for s in S:
        for k in K:
            model.addConstr(quicksum(q[i,l] for l in L)\
                            ==quicksum(D[s] for s in BOM if K[i] in BOM[s]), name="A.7")
                
#subtour
for l in L:
    model.addConstr(u[0,l]==0)
    for i in range(len(N)+1):
        for j in N:
            if i!=j:
                model.addConstr(u[i,l]-u[j,l] + (len(N)+1)*x[i,j,l] <= (len(N)+1), name="A.8")


#Objective
model.setObjective(quicksum(R[i,j]*x[i,j,l] for i in N for j in N for l in L if i!=j), GRB.MINIMIZE)

model.optimize()

I would like to have the solution go as 0-5 locations-0 as 1 big circle. But now it's two circel (0-5-0 / 1-2-3-4-1), as below.

How can I modify the code to solve this issue?