Gurobiを用いたTSP問題の解法

久保ら「あたらしい数理最適化　Python言語とGurobiで解く」（近代科学社）に掲載されているプログラムは，次のサイトに掲載されている．

あたらしい数理最適化 -GurobiとPython言語で解く- サポートページ

の「ソースコードzipファイル」にある．

しかし，ソースコードは，python 2.Xを想定している．しかし，現在のGurobiはpython 3系列に移行したので，print文などの変更が必要になる．

特にzipファイルのtsp.py（あたらしい数理最適化p.80-86）については，修正が必要だったので，以下に掲載している．なお，プログラム中，使っていないaddcut2の関数については省略した．

"""
tsp.py:  solve the traveling salesman problem 

minimize the travel cost for visiting n customers exactly once
approach:
	- start with assignment model
	- add cuts until there are no sub-cycles

Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
"""
import math
import random
import networkx
from gurobipy import *

def solve_tsp(V,c):
	"""solve_tsp -- solve the traveling salesman problem 
		- start with assignment model
		- add cuts until there are no sub-cycles
	Parameters:
		- V: set/list of nodes in the graph
		- c[i,j]: cost for traversing edge (i,j)
	Returns the optimum objective value and the list of edges used.
	"""
	
	def addcut(cut_edges):
		G = networkx.Graph()
		G.add_edges_from(cut_edges)
		Components = list(networkx.connected_components(G))	#listが必要

		if len(Components) == 1:
			return False
		for S in Components:
			model.addConstr(quicksum(x[i,j] for i in S for j in S if j>i) <= len(S)-1)
			print("cut: len(%s) <= %s" % (S,len(S)-1))
		return True


	# main part of the solution process:
	model = Model("tsp")

	model.Params.OutputFlag = 0 # silent/verbose mode
	x = {}
	for i in V:
		for j in V:
			if j > i:
				x[i,j] = model.addVar(ub=1, name="x(%s,%s)"%(i,j))
	model.update()
	
	for i in V:
		model.addConstr(quicksum(x[j,i] for j in V if j < i) + \
						quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i)

	model.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j > i), GRB.MINIMIZE)

	EPS = 1.e-6
	while True:
		model.optimize()
		edges = []
		for (i,j) in x:
			if x[i,j].X > EPS:
				edges.append( (i,j) )

		if addcut(edges) == False:
			if model.IsMIP:	 # integer variables, components connected: solution found
				break
			for (i,j) in x:	 # all components connected, switch to integer model
				x[i,j].VType = "B"
			model.update()

	return model.ObjVal,edges


def distance(x1,y1,x2,y2):
	"""distance: euclidean distance between (x1,y1) and (x2,y2)"""
	return math.sqrt((x2-x1)**2 + (y2-y1)**2)

def make_data(n):
	"""make_data: compute matrix distance based on euclidean distance"""
	V = range(1,n+1)
	x = dict([(i,random.random()) for i in V])
	y = dict([(i,random.random()) for i in V])
	c = {}
	for i in V:
		for j in V:
			if j > i:
				c[i,j] = distance(x[i],y[i],x[j],y[j])
	return V,c
	
			
if __name__ == "__main__":
	import sys

	# Parse argument
	if len(sys.argv) < 2:
		print("Usage: %s instance" % sys.argv[0])
		print("Using randomized example instead")
		n = 200
		seed = 1
		random.seed(seed)
		V,c = make_data(n)
	else:
		from read_tsplib import read_tsplib
		try:
			V,c,x,y = read_tsplib(sys.argv[1])
		except:
			print("Cannot read TSPLIB file",sys.argv[1])
			exit(1)

	obj,edges = solve_tsp(V,c)

	print("\nOptimal tour:",edges)
	print("Optimal cost:",obj)

Gurobiではなく，PySCIPOptを使って，似たようなことを行っている人がいた．

GitHubに次のサイトがある．

examplesのフォルダーのなかのfinishedに，「あたらしい数理最適化」のプログラムがある．その中のtsp.pyは，上記のプログラムに対して，PySCIPOptを使うようにしたプログラムになっている．上記のプログラムとPySCIPOptのtsp.pyも同じ出力結果であることを確認した．