How to implement a multiobjective TSP?

[grafkevin]

Hi everyone,

This question is a rather long one. It is referring to this example code: https://github.com/Project-Platypus/Platypus/blob/master/examples/tsp.py

I am working on a multiobjective TSP model where the salesman can choose between bus or flight to go from city i to j. The two conflicting objectives are to minimise travelling time and carbon emissions. This would be applicable to a tourist doing a citytrip, the tourist has the choice whether to go by plane (fast but carbon intensive) or by bus (slow but green).

I created my own "distance matrices" for carbon emissions and traveling time for 6 cities. This dataset I would like to feed the NSGA-II algorithm to find out the non-dominated solutions. Below is just the carbon emission matrix:

carbon_bus = {'Berlin': [999999.0, 4.58, 8.31, 20.71, 35.65, 13.89], 
              'Dresden': [4.58, 999999.0, 3.54, 15.96, 32.68, 10.92], 
              'Prague': [8.31, 3.54, 999999.0, 12.47, 30.88, 9.1], 
              'Budapest': [20.71, 15.96, 12.47, 999999.0, 28.83, 15.48], 
              'Rome': [35.65, 32.68, 30.88, 28.83, 999999.0, 21.73], 
              'Munich': [13.89, 10.92, 9.1, 15.48, 21.73, 999999.0]}

carbon_bus = pd.DataFrame(carbon_bus)
carbon_bus["carbon_bus(kg CO2)"] = ["Berlin","Dresden","Prague","Budapest","Rome","Munich"]
carbon_bus = carbon_bus.set_index("carbon_bus(kg CO2)")
print(carbon_bus)

I tried writing my emissions function based on the tsp()-function from the example code (https://github.com/Project-Platypus/Platypus/blob/master/examples/tsp.py):

def emissions(x):
    em = 0
    tour = x[0]
    for i in range(len(tour)-1):
        em = em + carbon_bus.iloc[tour[i],tour[i+1%len(tour)]]
    return em

I tried to feed this function into the solver but failed to get a reasonable solution:

problem = Problem(1, 1)
problem.types[0] = Permutation(range(len(6))) 
problem.directions[0] = Problem.MINIMIZE
problem.function = emissions
algorithm = GeneticAlgorithm(problem)
algorithm.run(1000, callback = lambda a : print(a.nfe, unique(nondominated(algorithm.result))[0]))

Questions:

1. How can I minimise the carbon objective by referring to the matrix given above?
2. Also, I don't know how to understand the line with tour=x[0] which I have taken from the tsp()-function. Is this referring to the decision variable x_ij?
3. How would you propose to implement a multiobjective TSP in playtypus?

Thanks in advance
Kevin

For more information on my model you can refer to my question in OR-Stackexchange: https://or.stackexchange.com/q/3933/3385

[dhadka]

1. There's a bug on this line: em = em + carbon_bus.iloc[tour[i],tour[i+1%len(tour)]]. Put parenthesis around (i+1).

2. x is an array of decision variables. In this case, you have a single decision variable, a permutation. x[0] is reading the value of that permutation, which is an array with the permutation values, e.g., [1, 5, 4, 2, 0, 3].

3. I've found it extremely helpful when solving TSP problems to include use the 2-opt local search heuristic. I expect you'd see substantially better results using this. It's not included in this example for simplicity.

[grafkevin]

Hi David, thank for responding so quickly on my three questions!

About Q1: How can I minimise the carbon objective by referring to the matrix given above?

I still have trouble obtaining a result. Here is my code:

carbon_bus = {'Berlin': [999999.0, 4.58, 8.31, 20.71, 35.65, 13.89], 
              'Dresden': [4.58, 999999.0, 3.54, 15.96, 32.68, 10.92], 
              'Prague': [8.31, 3.54, 999999.0, 12.47, 30.88, 9.1], 
              'Budapest': [20.71, 15.96, 12.47, 999999.0, 28.83, 15.48], 
              'Rome': [35.65, 32.68, 30.88, 28.83, 999999.0, 21.73], 
              'Munich': [13.89, 10.92, 9.1, 15.48, 21.73, 999999.0]}

carbon_bus = pd.DataFrame(carbon_bus)
carbon_bus["carbon_bus(kg CO2)"] = ["Berlin","Dresden","Prague","Budapest","Rome","Munich"]
carbon_bus = carbon_bus.set_index("carbon_bus(kg CO2)")

def emissions(x):
    em = 0
    tour = x[0]
    for i in range(len(tour)-1):
        em = em + carbon_bus.iloc[tour[i],tour[(i+1)%len(tour)]]
    return em

problem = Problem(1, 1)
problem.types[0] = Permutation(range(6)) 
problem.directions[0] = Problem.MINIMIZE
problem.function = emissions

algorithm = GeneticAlgorithm(problem)
algorithm.run(100000, callback = lambda a : print(a.nfe, unique(nondominated(algorithm.result))[0].objectives[0]))

And I get this error message:

c:\users\kevin\appdata\local\programs\python\python37\lib\site-packages\platypus\core.py in evaluate(self, solution)
    198             constrs = [constrs]
    199 
--> 200         if len(objs) != self.nobjs:
    201             raise PlatypusError("incorrect number of objectives: expected %d, received %d" % (self.nobjs, len(objs)))
    202 

TypeError: object of type 'numpy.float64' has no len()

About Q3: How would you propose to implement a multiobjective TSP in Platypus?

Thanks for the recommendation on the 2-opt operator, will try it out definitely. But I was more asking on how write the multiobjective TSP in Platypus. I am not a very experienced programmer as you can probably tell. Therefore, I would be very thankful if you could give me some inspiration by providing a rough outline on how to setup the code for a multiobjective TSP in Platypus.

I have another matrix with the traveling times between the 6 exemplary cities. I suppose I need to write another function which calculates the total traveling time for a given tour. The problem object is then changed as follows:

problem = Problem(nvars=1, nobjs=2)
problem.types[0] = Permutation(range(6)) 
problem.directions[0] = Problem.MINIMIZE
problem.function = emissions, time

Am I thinking in the right direction here? In the documentation I didn't find much about implementing multiple objective functions unfortunately.

[dhadka]

1. Try return em.item(). This will convert from the NumPy float64 value it's currently returning to a Python float.

2. For the multi-objective case, you'll need to make a single function that returns the two values:

def fitness(x):
    em = 0
    time = 0
    tour = x[0]
    for i in range(len(tour)-1):
        em = em + carbon_bus.iloc[tour[i],tour[(i+1)%len(tour)]]
        time = time + ...
    return [em, time]

and just set

problem.function = fitness

[grafkevin]

Hi David, thanks for all the help so far! I managed to implement the multiobjective function.

Now, I am trying to use local search operators like two_opt which you recommended. I guess having to rely only on random permutatations is not quite efficient. I have implemented the two_opt function in the normal TSP first. It works for me but it's taking quite some time to complete the function evaluations. Could you please take a look and tell me what could be improved?

def two_opt(x):
    route = x[0]
    bestdist = distance(route)
    bestr = route.copy()
    for i in range(0, len(route)-3):
        for j in range(i+2, len(route)-1):
            newr = route[:i+1] + route[j:i:-1] + route[j+1:] 
            newdist = distance(newr)
            if newdist < bestdist:
                bestr = newr.copy()
                bestdist = newdist
    return bestdist

problem = Problem(nvars=1, nobjs=1)
problem.types[0] = Permutation(range(len(cities)))
problem.directions[0] = Problem.MINIMIZE
problem.function = two_opt

progress = []
pop_size = 50

algorithm = GeneticAlgorithm(problem, population_size=pop_size, offspring_size=50)
algorithm.run(500, callback = lambda a: (print("Function evaluations:",a.nfe,algorithm.result[0].objectives[0]),
                                          progress.append(unique(nondominated(algorithm.result))[0])))

Also, I was wondering: Did you ever achieve to get the optimal route for PR76 with Platypus?

Thanks again
Kevin

[dhadka]

This is the version of 2-opt that I wrote in Java - https://github.com/MOEAFramework/MOEAFramework/blob/master/examples/org/moeaframework/examples/ga/tsplib/TSP2OptHeuristic.java. It's been a few years since I wrote that so I forget the details, but one difference I see is that you don't need to compute the full route each time, you can just compute the distance between the four cities (i, i+1, j, j+1).

I don't think I've ever seen the optimal result for PR76 without 2-opt, but that's not surprising. There are 76! = 1.88 x 10^111 different permutations. Without some local search, it's really, really hard to find the optimum. With 2-opt, you can find the optimum in 10-20 iterations. The MOEAFramework, the Java library I wrote that Platypus is based on, has all this implemented.

[dhadka]

Adding to ☝️ , it also can plot the search, showing the best route found so far (red), and the current population (gray lines) - http://moeaframework.org/images/screenshot-5.png

[grafkevin]

Hi David, I implemented your version of two_opt which really seems to work better. Now, I have some problem with the distance matrix which is needed for the two_opt operator. The get_total_distance function is the objective function to be optimized. Here's my code:

def two_opt(x):
    #route = x[0]
    d = distMatrix
    for i in range(0,len(route)-3):
        for j in range(i+2,len(route)-1):
            dist1 = d[route[i]][route[i+1]] + d[route[j]][route[j+1]]
            dist2 = d[route[i]][route[j]] + d[route[i+1]][route[j+1]]

            if dist2 < dist1:
                route = route[:i+1] + route[j:i:-1] + route[j+1:]
                print("route",route)
                dist = get_total_distance(route)
                print("dist",dist)
    return dist

def read_coordinates(inputfile):
    coord = [] 
    iFile = open(inputfile, "r")
    for i in range(6):           # skip first 6 lines
        iFile.readline()
    line = iFile.readline().strip()
    while line != "EOF":
        values = line.split()
        coord.append([float(values[1]), float(values[2])])
        line = iFile.readline().strip() 
    iFile.close()
    return coord

def get_distance_matrix(coord):
    distMatrix = squareform(pdist(coord,"euclidean"))   # returns an array
    return distMatrix 

def get_total_distance(x):
    route = x[0]
    total = 0
    d = distMatrix
    for i in range(len(route)-1):
        total += d[route[i]][route[i+1]]
    return total

Error:

c:\users\kevin\appdata\local\programs\python\python37\lib\site-packages\platypus\core.py in evaluate(self, solution)
    198             constrs = [constrs]
    199 
--> 200         if len(objs) != self.nobjs:
    201             raise PlatypusError("incorrect number of objectives: expected %d, received %d" % (self.nobjs, len(objs)))
    202 

TypeError: object of type 'numpy.float64' has no len()

The code runs fine when I feed the function a predefined route tough. Where is this error coming from?

Thanks again for the help!
Kevin

[dhadka]

Have it return total.item(). total is a numpy float, and needs to be converted to a Python float.

[grafkevin]

It works, thanks! The two_opt is quite fast and is getting me good results. But the algorithm.result[0] object does not give out the corresponding route to the given distance. I tried several ways to solve it but failed.

def read_coordinates(inputfile):
    coord = [] 
    iFile = open(inputfile, "r")
    for i in range(6):           # skip first 6 lines
        iFile.readline()
    line = iFile.readline().strip()
    while line != "EOF":
        values = line.split()
        coord.append([float(values[1]), float(values[2])])
        line = iFile.readline().strip() 
    iFile.close()
    return coord

def get_distance_matrix(coord):
    distMatrix = squareform(pdist(coord,"euclidean"))  # returns an array
    return distMatrix 

def get_total_distance(x):
    route = x
    total = 0
    d = distMatrix
    for i in range(len(route)-1):
        total += d[route[i]][route[i+1]]
    return total.item()

def two_opt(x):
    route = x[0]
    dist = get_total_distance(route)
    d = distMatrix
    for i in range(0,len(route)-3):
        for j in range(i+2,len(route)-1):
            dist1 = d[route[i]][route[i+1]] + d[route[j]][route[j+1]]
            dist2 = d[route[i]][route[j]] + d[route[i+1]][route[j+1]]
    
            if dist2 < dist1:
                route = route[:i+1] + route[j:i:-1] + route[j+1:]
                dist = get_total_distance(route)
    return dist

coord = read_coordinates("instances\pr76.tsp.txt")
distMatrix = get_distance_matrix(coord)
progress = []

problem = Problem(nvars=1, nobjs=1)
problem.types[0] = Permutation(range(len(cities)))
problem.directions[0] = Problem.MINIMIZE
problem.function = two_opt

algorithm = GeneticAlgorithm(problem)
algorithm.run(10000, callback = lambda a: (print(a.nfe, unique(nondominated(algorithm.result[0])),
                                           progress.append(unique(nondominated(algorithm.result))[0]))))

I get a result like this:
800 [Solution[[71, 47, 15, 55, 33, 69, 63, 35, 64, 72, 43, 1, 31, 56, 13, 50, 3, 46, 60, 61, 57, 26, 58, 22, 74, 2, 11, 9, 29, 10, 42, 37, 40, 59, 68, 48, 4, 21, 49, 5, 62, 27, 39, 45, 7, 44, 36, 70, 12, 30, 23, 14, 34, 24, 52, 6, 28, 20, 65, 18, 17, 16, 25, 54, 67, 0, 73, 53, 8, 32, 41, 51, 19, 38, 66, 75]|117040.59065613162|0]] None

When I plug this route into the get_total_distance() function shows me a much longer distance. I know I can print the routes in the two_opt() function to give me right routes but that's a workaround I try to avoid.

PR76 file

[github-actions]

This issue is stale and will be closed soon. If you feel this issue is still relevant, please comment to keep it active. Please also consider working on a fix and submitting a PR.