op-Traveling Salesman Problem

TSP(Traveling Salesman Problem)

TSP 개념 및 응용분야

+) tsp 테스트 DATA 링크
https://www.math.uwaterloo.ca/tsp/world/index.html

TSP 모델링

TSP 코드(예시)

from ortools.linear_solver import pywraplp

nCity = 6

DIST = [ [0, 702, 454, 842, 2396, 1196],
[702, 0, 324, 1093, 2136, 764],
[454, 324, 0, 1137, 2180, 798],
[842, 1093, 1137, 0, 1616, 1857],
[2396, 2136, 2180, 1616, 0, 2900],
[1196, 764, 798, 1857, 2900, 0]
]

solver = pywraplp.Solver.CreateSolver("SAT")
X = {}
for i in range(nCity):
    for j in range(nCity):
        if i != j:
            X[i, j] = solver.IntVar(0, 1, "X"+str(i)+str(j))

U = {}
for i in range(1, nCity):
    U[i] = solver.IntVar(1, nCity-1, "U[%i]" %i)


# 제약조건
# 도시 j로 한 번은 들어와야 한다.
for j in range(nCity):
    solver.Add(sum([X[i,j] for i in range(nCity) if i!=j]) == 1)

# 도시 j로 부터 한 번은 나와야 한다.
for i in range(nCity):
    solver.Add(sum([X[i,j] for j in range(nCity) if i!=j]) == 1)

# 방문 순서 제약

for i in range(1, nCity):
    for j in range(1, nCity):
        if i!=j:
            solver.Add(U[i]-U[j] + nCity*X[i,j] <= nCity-1)


# Objective
objective_terms = []
for i in range(nCity):
    for j in range(nCity):
        if i != j:
            objective_terms.append(DIST[i][j] * X[i, j])

solver.Minimize(solver.Sum(objective_terms))


if 1:
    with open('or9-1.lp', "w") as out_f:
        lp_text = solver.ExportModelAsLpFormat(False)
        out_f.write(lp_text)
# Solve
status = solver.Solve()

# Print solution.
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print(f"Total cost = {solver.Objective().Value():.1f}\n", )
    for i in range(nCity):
        for j in range(nCity):
            if i != j:
                if X[i, j].solution_value() > 0.5:
                    print(f"X{i} --> X{j}")
    for i in range(1, nCity):
        print(f"{i} 도시 방문순서: ", U[i].solution_value())
else:
    print("No solution found.")

TSP 상황별 해법

1. TSP – book salesperson

# 2번만 풀이
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

def create_data_model():
    data = {}
    data["distance_matrix"] = [
    [0, 125, 225, 155, 215],
    [125, 0, 85, 115, 135],
    [225, 85, 0, 165, 190],
    [155, 115, 165, 0, 195],
    [215, 135, 190, 195, 0]
    ]
    data["num_vehicles"] = 1
    data["depot"] = 0
    return data

data = create_data_model()

# 색인 관리자, 일종의 정보관리자
manager = pywrapcp.RoutingIndexManager(len(data["distance_matrix"]), data["num_vehicles"], data["depot"])

# Routing 객체 생성
routing = pywrapcp.RoutingModel(manager)

# 두 점 사이의 거리를 반환한다.
def distance_callback(from_index, to_index):
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    return data["distance_matrix"][from_node][to_node]

# 거리 계산 함수를 callback 함수로 등록하고 활용함.
transit_callback_index = routing.RegisterTransitCallback(distance_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

# 시작해를 찾기 위한 휴리스틱 메소드를 등록
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

def print_solution(manager, routing, solution):
    print(f"Objective: {solution.ObjectiveValue()} miles")
    index = routing.Start(0)
    plan_output = "Route for vehicle 0:\n"
    route_distance = 0
    while not routing.IsEnd(index):
        plan_output += f" {manager.IndexToNode(index)} ->"
        previous_index = index
        index = solution.Value(routing.NextVar(index))
        route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
    plan_output += f" {manager.IndexToNode(index)}\n"
    print(plan_output)
    plan_output += f"Route distance: {route_distance}miles\n"


# 문제 풀이
solution = routing.SolveWithParameters(search_parameters)
if solution:
    print_solution(manager, routing, solution)

# 해 경로를 리스트에 저장하기 (Optional)
def get_routes(solution, routing, manager):
    routes = []
    for route_nbr in range(routing.vehicles()):
        index = routing.Start(route_nbr)
        route = [manager.IndexToNode(index)]
        while not routing.IsEnd(index):
            index = solution.Value(routing.NextVar(index))
            route.append(manager.IndexToNode(index))
        routes.append(route)
    return routes

routes = get_routes(solution, routing, manager)

# Display the routes.
for i, route in enumerate(routes):
    print('Route', i, route)

+) 3번 질문에 대한 답 : nearest neighbor algorithm
경로 : basin -> wald -> bon -> mena -> kiln -> basin / 거리 = 785miles

"nearest neighbor algorithm" 방식이 거리가 더 크게 나온다.(최적해가 아님)

이유 : basin에서 kiln까지 가는과정의 거리는 최소일지 모르나 다시 basin으로 돌아와야 하는 과정에서 비효율적인 거리가 발생하기 때문인것으로 보인다

2. TSP – Open Tour

1) TSP – Open Tour 개념정리

2) TSP – Open Tour 실습문제

from ortools.linear_solver import pywraplp
# 2번 풀이
# 추가적인 9행과 9열을 만들어서 의미없는 0을 집어넣어야 open-tsp문제를 해결할 수 있다
DIST =[
    [0, 20, 30, 25, 12, 33, 44, 57, 0],
    [22, 0, 19, 20, 20, 29, 43, 45, 0],
    [28, 19, 0, 17, 38, 48, 55, 60, 0],
    [25, 20, 19, 0, 28, 35, 40, 55, 0],
    [12, 18, 34, 25, 0, 21, 30, 40, 0],
    [35, 25, 45, 30, 20, 0, 25, 39, 0],
    [47, 39, 50, 35, 28, 20, 0, 28, 0],
    [60, 38, 54, 50, 33, 40, 25, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0]]

nCity = len(DIST)

solver = pywraplp.Solver.CreateSolver("SAT")
X = {}
for i in range(nCity):
    for j in range(nCity):
        if i != j:
            X[i, j] = solver.IntVar(0, 1, "X"+str(i)+str(j))

U = {}
for i in range(1, nCity):
    U[i] = solver.IntVar(1, nCity-1, "U[%i]" %i)


# 제약조건
# 도시 j로 한 번은 들어와야 한다.
for j in range(nCity):
    solver.Add(sum([X[i,j] for i in range(nCity) if i!=j]) == 1)

# 도시 j로 부터 한 번은 나와야 한다.
for i in range(nCity):
    solver.Add(sum([X[i,j] for j in range(nCity) if i!=j]) == 1)

# 방문 순서 제약

for i in range(1, nCity):
    for j in range(1, nCity):
        if i!=j:
            solver.Add(U[i]-U[j] + nCity*X[i,j] <= nCity-1)
solver.Add(X[5,8] == 1)

# Objective
objective_terms = []
for i in range(nCity):
    for j in range(nCity):
        if i != j:
            objective_terms.append(DIST[i][j] * X[i, j])

solver.Minimize(solver.Sum(objective_terms))


if 1:
    with open('과제4', "w") as out_f:
        lp_text = solver.ExportModelAsLpFormat(False)
        out_f.write(lp_text)
# Solve
status = solver.Solve()

# Print solution.
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print(f"Total dollar = {solver.Objective().Value():.1f}\n", )
    for i in range(nCity):
        for j in range(nCity):
            if i != j:
                if X[i, j].solution_value() > 0.5:
                    print(f"X{i} --> X{j}")
    for i in range(1, nCity):
        print(f"{i} 도시 방문순서: ", U[i].solution_value())
else:
    print("No solution found.")

3. TSP – 서비스 센터 : 추가 조건부여

from ortools.linear_solver import pywraplp

# DIST : 고객사이의 이동시간을 나타냄(단위 : 분)
DIST =[
[0, 20, 15, 19, 24, 14, 21, 11],
[20, 0, 18, 22, 23, 22, 9, 10],
[15, 18, 0, 11, 21, 14, 32, 12],
[19, 22, 11, 0, 20, 27, 18, 15],
[24, 23, 21, 20, 0, 14, 25, 20],
[14, 22, 14, 27, 14, 0, 26, 17],
[21, 9, 32, 18, 25, 26, 0, 20],
[11, 10, 12, 15, 20, 17, 20, 0]
]

nCity = len(DIST)

solver = pywraplp.Solver.CreateSolver("SAT")

# 고객 i -> 고객 j로 이동하는것을 표현하는 변수
X = {}
for i in range(nCity):
    for j in range(nCity):
        if i != j:
            X[i, j] = solver.IntVar(0, 1, "X"+str(i)+str(j))

# 
U = {}
for i in range(1, nCity):
    U[i] = solver.IntVar(1, nCity-1, "U[%i]" %i)


# 제약조건
# 도시 j로 한 번은 들어와야 한다.
for j in range(nCity):
    solver.Add(sum([X[i,j] for i in range(nCity) if i!=j]) == 1)

# 도시 j로 부터 한 번은 나와야 한다.
for i in range(nCity):
    solver.Add(sum([X[i,j] for j in range(nCity) if i!=j]) == 1)

# 방문 순서 제약

for i in range(1, nCity):
    for j in range(1, nCity):
        if i!=j:
            solver.Add(U[i]-U[j] + nCity*X[i,j] <= nCity-1)

# 고객방문순서 추가 제약조건
# 2,3번 고객이 먼저 해결되어야하나, 2번고객다음에 무조건 6번고객이 와야하므로 3번고객을 첫번째로 방문해야함 
solver.Add(U[3] == 1) 
solver.Add(U[2] == 2)
solver.Add(U[6] == 3)

# Objective
objective_terms = []
for i in range(nCity):
    for j in range(nCity):
        if i != j:
            objective_terms.append(DIST[i][j] * X[i, j])

solver.Minimize(solver.Sum(objective_terms))


if 1:
    with open('과제4-3번', "w") as out_f:
        lp_text = solver.ExportModelAsLpFormat(False)
        out_f.write(lp_text)
# Solve
status = solver.Solve()

# Print solution.
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print(f"Total time = {solver.Objective().Value():.1f}\n", )
    for i in range(nCity):
        for j in range(nCity):
            if i != j:
                if X[i, j].solution_value() > 0.5:
                    print(f"X{i} --> X{j}")
    for i in range(1, nCity):
        print(f"{i}번 고객 방문순서: ", U[i].solution_value())
else:
    print("No solution found.")

4. TSP – 미팅순서 : 행렬 스스로 작성

#%% 1번 : 행렬 스스로 작성
# 회의실에 들어오고 나가는 직원숫자로 행렬 DIST를 만들었다
# DIST[0][1] : project1 -> project2 수행 시 "회의실로 들어가는 직원수 + 나오는 직원수"
DIST =[
    [0, 4, 4, 6, 6, 5],
    [4, 0, 6, 4, 6, 3],
    [4, 6, 0, 4, 8, 7],
    [6, 4, 4, 0, 6, 5],
    [6, 6, 8, 6, 0, 5],
    [5, 3, 7, 5, 5, 0]]

#%% 2번 : 풀이
from ortools.linear_solver import pywraplp

# open-tsp문제이기 때문에 의미없는 7번째행과 7번째열을 생성하여 0값을 집어넣는다 
DIST =[
    [0, 4, 4, 6, 6, 5, 0],
    [4, 0, 6, 4, 6, 3, 0],
    [4, 6, 0, 4, 8, 7, 0],
    [6, 4, 4, 0, 6, 5, 0],
    [6, 6, 8, 6, 0, 5, 0],
    [5, 3, 7, 5, 5, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]]

nCity = len(DIST)

solver = pywraplp.Solver.CreateSolver("SAT")
X = {}
for i in range(nCity):
    for j in range(nCity):
        if i != j:
            X[i, j] = solver.IntVar(0, 1, "X"+str(i)+str(j))

U = {}
for i in range(1, nCity):
    U[i] = solver.IntVar(1, nCity-1, "U[%i]" %i)


# 제약조건
# 프로젝트 j로 한 번은 들어와야 한다.
for j in range(nCity):
    solver.Add(sum([X[i,j] for i in range(nCity) if i!=j]) == 1)

# 프로젝트 j로 부터 한 번은 나와야 한다.
for i in range(nCity):
    solver.Add(sum([X[i,j] for j in range(nCity) if i!=j]) == 1)


# 방문 순서 제약
for i in range(1, nCity):
    for j in range(1, nCity):
        if i!=j:
            solver.Add(U[i]-U[j] + nCity*X[i,j] <= nCity-1)

# Objective
objective_terms = []
for i in range(nCity):
    for j in range(nCity):
        if i != j:
            objective_terms.append(DIST[i][j] * X[i, j])

# 직원수 최소화하기
solver.Minimize(solver.Sum(objective_terms))


if 1:
    with open('과제4- 문제 4', "w") as out_f:
        lp_text = solver.ExportModelAsLpFormat(False)
        out_f.write(lp_text)
# Solve
status = solver.Solve()

# Print solution.
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print(f"총 직원수 = {solver.Objective().Value():.1f}\n", )
    for i in range(nCity):
        for j in range(nCity):
            if i != j:
                if X[i, j].solution_value() > 0.5:
                    print(f"X{i} --> X{j}")
    for i in range(1, nCity):
        print(f"프로젝트 {i} 미팅순서: ", U[i].solution_value())
else:
    print("No solution found.")

5. TSP – 기판드릴링(아직 이해못함;;)

# 기판드릴링 코드 / 뭐에쓰는지 아직 이해못함;;
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
import math
import pandas as pd

def makeDIST(nP):
    DIST = list()
    nCity = len(nP)
    for i in range(nCity):
        DIST.append([])
        for j in range(nCity):
            if j != i:
                temp = math.hypot(nP['xc'][i] - nP['xc'][j], nP['yc'][i] - nP['yc'][j])
                DIST[i].append(int(temp))
            else:
                DIST[i].append(0)
    return DIST


def create_data_model(DIST):
    data = {}
    data["distance_matrix"] = DIST
    data["num_vehicles"] = 1
    data["depot"] = 0
    return data

def print_solution(manager, routing, solution):
    print(f"Objective: {solution.ObjectiveValue()} miles")
    index = routing.Start(0)
    plan_output = "Route for vehicle 0:\n"
    route_distance = 0
    while not routing.IsEnd(index):
        plan_output += f" {manager.IndexToNode(index)} ->"
        previous_index = index
        index = solution.Value(routing.NextVar(index))
        route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
    plan_output += f" {manager.IndexToNode(index)}\n"
    print(plan_output)
    plan_output += f"Route distance: {route_distance}miles\n"

# Main Flow
#f = open(".\\TSP_data\\a280.tsp", 'r')
f = open("C:/Users/jjb69/OneDrive/바탕 화면/데이터기반최적화/TSP_data/a280.tsp", 'r') # 기판

nPos = pd.DataFrame(columns=[ 'xc', 'yc'])

flag = 0
while True:
    line = f.readline().strip()
    if line == "EOF": break
    if line == "NODE_COORD_SECTION":
        flag = 1
        continue
    if flag:
        ss = line.split()
        nPos.loc[len(nPos)] = [float(ss[1]), float((ss[2]))]
f.close()

# Draw a picture
import matplotlib.pyplot as plt
import numpy as np
plt.scatter(nPos['xc'], nPos['yc'])
plt.show()

DIST = makeDIST(nPos)

data = create_data_model(DIST)

manager = pywrapcp.RoutingIndexManager(len(data["distance_matrix"]), data["num_vehicles"], data["depot"])

# Create Routing Model.
routing = pywrapcp.RoutingModel(manager)

def distance_callback(from_index, to_index):
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    return data["distance_matrix"][from_node][to_node]

transit_callback_index = routing.RegisterTransitCallback(distance_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
solution = routing.SolveWithParameters(search_parameters)

# Print solution on console.
if solution:
    print_solution(manager, routing, solution)