mtz: model for the traveling salesman problem with time windows
    (based on Miller-Tucker-Zemlin's one-index potential formulation, stronger constraints)
    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.

def mtz2tw(n,c,e,l):
    """mtz: model for the traveling salesman problem with time windows
    (based on Miller-Tucker-Zemlin's one-index potential formulation, stronger constraints)
    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.
    """
    model = Model("tsptw - mtz-strong")
    
    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 j in range(2,n+1):
            if i != j:
                M1 = max(l[i] + c[i,j] - e[j], 0)
                M2 = max(l[i] + min(-c[j,i], e[j]-e[i]) - e[j], 0)
                model.addCons(u[i] + c[i,j] - M1*(1-x[i,j]) + M2*x[j,i] <= u[j], "LiftedMTZ(%s,%s)"%(i,j))

    for i in range(2,n+1):
        model.addCons(e[i] + quicksum(max(e[j]+c[j,i]-e[i],0) * x[j,i] for j in range(1,n+1) if i != j) \
                        <= u[i], "LiftedLB(%s)"%i)

        model.addCons(u[i] <= l[i] - \
                        quicksum(max(l[i]-l[j]+c[i,j],0) * x[i,j] for j in range(2,n+1) if i != j), \
                        "LiftedUB(%s)"%i)

    model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize")

    model.data = x,u
    return model