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


model = AbstractModel()

model.M = Set()
model.N = Set()

model.n = Param()
model.c = Param(model.M, model.N)

model.x = Var(model.M, model.N, domain=Binary)
model.u = Var(model.M, domain=NonNegativeIntegers)

def object(model):
    return sum(model.c[i,j]*model.x[i,j] for (i,j) in model.M*model.N if i!=j)

model.obj = Objective(rule=object)

def const1(model,j):
    return sum(model.x[i,j] for i in model.M if i!=j) == 1

model.cons = Constraint(model.N, rule= const1)

def const2(model,i):
    return sum(model.x[i,j] for j in model.N if j!=i) ==1

model.cons2 = Constraint(model.M, rule=const2)

def const3(model,i,j):
    if i==j or i <2 or j<2:
        return Constraint.Skip
    return model.u[i]-model.u[j]+model.n*model.x[i,j] <= model.n-1

model.cons3 = Constraint(model.M, model.N, rule=const3)

instance = model.create("salesman.dat")
instance.pprint()

opt = SolverFactory('glpk')

results = opt.solve(instance, tee=True)
results.write()
instance.solutions.load_from(results)

for v in instance.component_objects(Var, active=True):
    print ("Variable",v)
    varobject = getattr(instance, str(v))
    for index in varobject:
        print ("   ",index, varobject[index].value)
Solve Within PyCharm