pyomo: plotting the optimal route for tsp in pyomo

Question

I am still learning pyomo, and so far i have made some progress:

This link! gives an example of tsp in pyomo. I reproduced the code as below. And everything is working just fine. However, I am not able to print the optimal route, could someone help or give me a fair idea of how to print and plot the optimal route?

The code:

from pyomo.environ import * 
from pyomo.opt import SolverFactory
import pyomo.environ

n=13
distanceMatrix=[[0,8,4,10,12,9,15,8,11,5,9,4,10],
[8,0,27,6,8,6,17,10,12,9,8,7,5],
[4,7,0,7,9,5,8,5,4,8,6  ,10,8],
[10,6   ,7,0,6,11,5 ,9,8,12,11,6,9],
[12,8   ,19,6,   0,7,9,6,9,8,4,11,10],
[9,6,5,11,7,0,10,4,3,10,6,5,7],
[15,7   ,8,15,19,10,0,10,9,8,5,9,10],
[8,10   ,5,9,6,4,10,0,11,5,9,6,7],
[11,12,4,8, 19,13,9,11,0, 9,11,11,6],
[5,9,8,12,8,10,8,5,9,0,6,7,5],
   [9,8,6,11,14,6,5,9,11,6,0,10,7],
   [4,7,10,6,31,5,9,6,11,7,10,0,9],
   [10,5,8,9,10,7,10,7,6,5,7,9,0]] 
startCity = 0

model = ConcreteModel()
model.M = Set(initialize=range(1, n+1))
model.N = Set(initialize=range(1, n+1))
model.c = Param(model.N, model.N, initialize=lambda model, i, j:    distanceMatrix[i-1][j-1])
model.x = Var(model.N, model.N, within=Binary)


def obj_rule(model):  

    return sum(model.c[n,j]*model.x[n,j] for n in model.N for j in model.N)
model.obj = Objective(rule=obj_rule,sense=minimize)

def con_rule(model, n):
    return sum(model.x[j,n] for j in model.N if j < n) + sum(model.x[n,j]   for j in model.N if j > n) == 2

model.con = Constraint(model.N, rule=con_rule,doc='constraint1')
opt = SolverFactory("glpk")
results = opt.solve(model)
results.write()
print('Printing Values')
print(value(model.obj))

Answer 1

First you have to think about what form the solution takes. If we look at this imlementation, we see that we have the matrix x with dimensions NxN and the domain binary (0 or 1). Think about what it means for an element x[j][k] to be equal to 1. And what does it mean if j < k or j > k?

To extract the values for x, a simple way might be

import numpy as np
N = len(model.N)
x = np.zeros((N,N))
for (i,j), val in model.x.get_values().items():
    x[i-1,j-1] = val

Then you can use the values for x however you like.