code: def mtztw(n,c,e,l):
      
      model = Model("tsptw - mtz")
  
      x,u = {},{}
      for i in range(1,n+1):
          u[i] = model.addVar(lb=e[i], ub=l[i], vtype="C", name="u(%s)"%i)
          for j in range(1,n+1):
              if i != j:
                  x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))
  
      for i in range(1,n+1):
          model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i)
          model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i)
  
      for i in range(1,n+1):
          for j in range(2,n+1):
              if i != j:
                  M = max(l[i] + c[i,j] - e[j], 0)
                  model.addCons(u[i] - u[j] + M*x[i,j] <= M-c[i,j], "MTZ(%s,%s)"%(i,j))
  
      model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize")
  
      model.data = x,u
      return model
docs: mtzts: model for the traveling salesman problem with time windows
      (based on Miller-Tucker-Zemlin's one-index potential formulation)
      Parameters:
          - n: number of nodes
          - c[i,j]: cost for traversing arc (i,j)
          - e[i]: earliest date for visiting node i
          - l[i]: latest date for visiting node i
      Returns a model, ready to be solved.