time-dependent TSP as a ILP

[rafaelmarques7]

Hello,

this is not a question regarding a program stack error, but regarding an error about the solution value. Because of this, if this is the wrong forum to place this question, I will remove it.

I am using Google's Operation Research Tools to provide the solvers necessary to solve the time-dependent Traveling Salesman Problem. In order to do so, i follow the Integer Linear formulation for the TDTSP.

Basically, in this formulation, we have a 3-d cost matrix, with the indexes (i,j,t) and a list of decision variables, X, also with three indexes, X(i,j,t) that determine if the solution component x_ijt (the arc connecting cities i and j at the time t) belongs to the solution.

The objective is to minimize the total cost of the trip. The trip has to respect some constraints, that is: each city must be visited exactly once, and we return to the original city.

Below, is a figure representing the ILP forulation and the constraints. (Note: in constraint 5, the negative part of the summ, the index is XJit, and not XIjt)

Below I present my complete code to solve the problem:

    from __future__ import print_function
    from ortools.constraint_solver import pywrapcp
    import numpy
    import random_tdtsp

    def main(num_cities = 5):
        # Solver 
        parameters = pywrapcp.Solver.DefaultSolverParameters()
        solver = pywrapcp.Solver('simple_CP', parameters)

        # TDTSP object we will work with
        tdtsp_object = random_tdtsp.TDTSP_generator(num_cities)
        opt_sol, opt_path =  tdtsp_object['best_cost'], tdtsp_object['solution_path']

        # Cost     
        cost = tdtsp_object['cost_matrix']

        # Variables
        x = {}
        for i in range(0, num_cities):
            for j in range(0, num_cities):
                for t in range(0, num_cities):
                    x[i, j, t] = solver.IntVar(0, 1, 'x[%i, %i, %i]' % (i, j, t))

        # Transform x into an array
        x_array = [x[i,j,t] for i in range(0, num_cities) for j in range(0, num_cities) for t in range(0, num_cities)]

        # Constraints
        #2)
        for i in range(0, num_cities):
            solver.Add(solver.Sum(x[i,j,t] for j in range(0, num_cities) for t in range(0, num_cities)) == 1)

        #3)
        for j in range(0, num_cities):
            solver.Add(solver.Sum(x[i,j,t] for i in range(0, num_cities) for t in range(0, num_cities)) == 1)

        #4)
        for t in range(0, num_cities):
            solver.Add(solver.Sum(x[i,j,t] for i in range(0, num_cities) for j in range(0, num_cities)) == 1)

        #5)
        for i in range(1, num_cities):
            solver.Add(
                (
                    solver.Sum((t+1)*x[i,j,t] for j in range(0, num_cities) for t in range(1, num_cities) ) 
                    - solver.Sum((t+1)*x[j,i,t] for j in range(0, num_cities) for t in range(0, num_cities) )
                )
             == 1
            )

        # Objective
        obj_expr = solver.IntVar(0, 100000, "obj_expr")
        solver.Add(obj_expr == solver.Sum(x[i,j,t] * int(cost[i,j,t]) for i in range(0, num_cities) for j in range(0, num_cities) for t in range(0, num_cities)))
        objective = solver.Minimize(obj_expr, 1)

        # Collector 
        collector = solver.LastSolutionCollector()
        # add variables
        for var in x_array:
            collector.Add(var)
        # add objective
        collector.AddObjective(obj_expr)

        # Decision builder
        db = solver.Phase(x_array, solver.CHOOSE_FIRST_UNBOUND, solver.ASSIGN_MIN_VALUE)

        # Solve 
        solver.Solve(db, [objective, collector])

        if collector.SolutionCount() > 0:
            best_solution = collector.SolutionCount() -1
            print('objective value: ', collector.ObjectiveValue(best_solution))
            print('optimal obj val: ', opt_sol)
            print('path: ')
            for var in x_array:
                if collector.Value(best_solution, var) == 1:
                    print(var)
            print('optimal path: ', opt_path)
        else:
            print('no solution')

    if __name__ == "__main__":
        main()

furthermore, in order to compare the solution, I used a function that I previously had, which generated a TDTSP instance, with a known optimal solution the code of this function is below.

    from random import shuffle, randint
    import numpy as np

    def TDTSP_generator(num_cities = 10):
        "returns a TDTSP object with: 'cost_matrix' (np.matrix), 'solution_path', 'best_cost'"
        #print("generating tdtsp instance with " + str(num_cities) + " cities")
        #1) specify an optimum tour
        opt_sol = [0]
        part_sol = [i for i in range(1, num_cities)]
        shuffle(part_sol) #holds a shuffled list
        opt_sol += part_sol
        opt_sol.append(opt_sol[0])

        #2) randomly assign values to the dual variables, alfa, beta, gama
        min_int = 0
        max_int = 100
        alfa = [randint(min_int, max_int) for i in range(num_cities)]
        beta = [randint(min_int, max_int) for i in range(num_cities)]
        gama = [randint(min_int, max_int) for i in range(num_cities)]

        #allocate memory for the cost matrix
        cmat = np.zeros([num_cities, num_cities, num_cities])

        #3) select Cijt that satisfies the CS condition for the basic variables
        #i.e. Cit,jt+1,t = alfa_it + beta_it+1 + gama_t
        for t in range(0, num_cities):
            i = opt_sol[t]
            j = opt_sol[t+1]
            cmat[t][i,j] = alfa[i] + beta[j] + gama[t]

        #4) select Cijt that satisfies the dual feasibility for non basic bariables
        #ii.e. Cijt >= alfa_i + beta_j + gama_t
        for t in range(0, num_cities):
            for i in range(0, num_cities):
                for j in range(0, num_cities):
                    if cmat[t][i,j] == 0:
                        min_val = alfa[i] + beta[j] + gama[t]
                        cmat[t][i,j] = randint(min_val, min_val + 100)

        #solution cost
        cost = 0
        for t in range(0, num_cities):
            i = opt_sol[t]
            j = opt_sol[t+1]
            cost += cmat[t][i,j]

        #return
        return_data = {
            "cost_matrix": cmat,
            "solution_path": opt_sol,
            "best_cost": cost,
        }
        return return_data

    if __name__ == "__main__":
        data = TDTSP_generator()
        print(data)

Note: the solutions generated by the ILP formulation presented, always start and end at node 0. The code I had to generate a random TDTSP instance with a known solution was for a random depot, and I made a slight alteration in order to start always at 0. but i dont believe it influences the results)

The problems I am experiencing are:

1) The solution presented by the ILP solver does not coincide with the optimal solution (the error is of (normaly) one wrong solution component )

2) How do I print the route of the solution? What I am doing now is printing the IntVar which have a 1 value (belonging thus to the solution, which displays something like x0,5,0, where the [0,5,0] means, from city 0 to city 5 at time 0, and the (0..1) are the allowed values for the variable.

I appreciate your help, and I will be available to provide any further feedback.

[lperron]

The title says solve with ILP, and then you use the CP solver. Therefore, either you use a MIP solver (preferably a good one), or you use the CP constraint (NoCycle in case).

You can also look at the CP-SAT solver (ortools/sat/doc/index.md) with the circuit constraint.

Anyway, the recommended way is to use the routing library that is tailored for this.

[BRaabe99]

Hi rafaelmarques7! Did you manage to solve this issue efficiently? I used your model, converted it to a CP-SAT Model and solved a 42 node Problem, but it took 2h and i need it to be able to solve ~150 nodes. Did you manage to use the circuit contraint correctly or did you find a solution in the routing library?