Main Content

여행하는 외판원 문제: 문제 기반

이 예제에서는 이진 정수 계획법을 사용하여 전형적인 여행하는 외판원 문제를 푸는 방법을 보여줍니다. 이 문제는 여러 경유지(도시)를 거쳐 원래 도시로 다시 돌아오는 최단 순회 경로를 구하는 것이 목적입니다. 이 경우에는 200개의 경유지가 있지만, 사용자는 nStops 변수를 변경하여 손쉽게 문제의 규모를 정할 수 있습니다. 초기 문제를 풀면 해에 하위 경로가 있음을 확인할 수 있습니다. 즉, 문제에서 구한 최적해가 모든 점을 통과하는 하나의 연속적인 경로를 제공하는 것이 아니라, 연속되지 않은 여러 개의 루프를 가지고 있음을 의미합니다. 그런 다음, 하위 경로를 확인하고 제약 조건을 추가한 후 하위 경로가 제거될 때까지 최적화를 다시 실행하는 반복적인 작업을 하게 될 것입니다.

이 문제에 대한 솔버 기반 접근법은 여행하는 외판원 문제: 솔버 기반 항목을 참조하십시오.

문제 정식화

다음과 같이 정수 선형 계획법을 사용하는 여행하는 외판원 문제를 정식화합니다.

  • 가능한 모든 여행(trip), 즉 서로 다른 모든 경유지 쌍을 생성합니다.

  • 각 여행의 거리를 계산합니다.

  • 최소화할 비용 함수는 경로에 포함된 각 여행 거리의 합입니다.

  • 결정 변수는 이진 변수이고 각 여행과 연결되어 있습니다. 여기서 1은 경로에 존재하는 여행을 나타내고, 0은 경로에 존재하지 않는 여행을 나타냅니다.

  • 경로에 모든 경유지가 포함되도록 하기 위해 각 경유지가 정확히 두 개의 여행에 존재하도록 규정하는 선형 제약 조건을 포함시켜야 합니다. 이는 각 경유지에서 도착과 출발이 각각 한 번씩 일어남을 의미합니다.

경유지 생성하기

대략적으로 미대륙을 나타내는 다각형 내부에 임의 경유지를 생성합니다.

load('usborder.mat','x','y','xx','yy');
rng(3,'twister') % Makes stops in Maine & Florida, and is reproducible
nStops = 200; % You can use any number, but the problem size scales as N^2
stopsLon = zeros(nStops,1); % Allocate x-coordinates of nStops
stopsLat = stopsLon; % Allocate y-coordinates
n = 1;
while (n <= nStops)
    xp = rand*1.5;
    yp = rand;
    if inpolygon(xp,yp,x,y) % Test if inside the border
        stopsLon(n) = xp;
        stopsLat(n) = yp;
        n = n+1;
    end
end

점 간의 거리 계산하기

200개의 경유지가 있으므로 여행 횟수는 19,900개입니다. 즉, 이진 변수가 19,900개임을 의미합니다(변수 개수 = 200개 중에서 2개를 선택하는 경우의 수).

모든 여행, 즉 모든 경유지 쌍을 생성합니다.

idxs = nchoosek(1:nStops,2);

피타고라스 규칙을 사용할 수 있도록 지구가 평평하다고 간주하고 모든 여행 거리를 계산합니다.

dist = hypot(stopsLat(idxs(:,1)) - stopsLat(idxs(:,2)), ...
             stopsLon(idxs(:,1)) - stopsLon(idxs(:,2)));
lendist = length(dist);

위와 같이 dist 벡터에 대한 정의를 사용하면 경로 길이는 다음과 같습니다.

dist'*trips

여기서 trips는 해에서 구한 여행 횟수를 나타내는 이진 벡터입니다. 이 경로의 길이가 바로 최소화하고자 하는 경로의 거리입니다.

그래프 생성 및 지도 그리기

문제를 그래프로 나타냅니다. 경유지를 노드로 나타내고 여행을 간선으로 나타내는 그래프를 생성합니다.

G = graph(idxs(:,1),idxs(:,2));

그래프 플롯을 사용하여 경유지를 표시합니다. 그래프 간선 없이 노드를 플로팅합니다.

figure
hGraph = plot(G,'XData',stopsLon,'YData',stopsLat,'LineStyle','none','NodeLabel',{});
hold on
% Draw the outside border
plot(x,y,'r-')
hold off

변수와 문제 생성하기

잠재적 여행을 나타내는 이진 최적화 변수를 갖는 최적화 문제를 생성합니다.

tsp = optimproblem;
trips = optimvar('trips',lendist,1,'Type','integer','LowerBound',0,'UpperBound',1);

문제에 목적 함수를 포함시킵니다.

tsp.Objective = dist'*trips;

제약 조건

각 경유지마다 도착하는 여행과 출발하는 여행이 각각 하나씩 있어야 하므로, 경유지마다 두 개의 여행이 연결되어 있는 선형 제약 조건을 생성합니다.

그래프 표현을 사용해 경유지에 연결된 모든 간선을 찾아 해당 경유지에서 시작하거나 끝나는 모든 여행을 식별합니다. 경유지마다 해당 경유지에 대한 여행의 합을 2로 규정하는 제약 조건을 생성합니다.

constr2trips = optimconstr(nStops,1);
for stop = 1:nStops
    whichIdxs = outedges(G,stop); % Identify trips associated with the stop
    constr2trips(stop) = sum(trips(whichIdxs)) == 2;
end
tsp.Constraints.constr2trips = constr2trips;

초기 문제 풀기

이제, 문제를 풀 준비가 되었습니다. 반복 출력값을 표시하지 않으려면 디폴트 표시를 끄십시오.

opts = optimoptions('intlinprog','Display','off');
tspsol = solve(tsp,'options',opts)
tspsol = struct with fields:
    trips: [19900×1 double]

해 시각화하기

해에서 구한 여행을 간선으로 사용하여 새 그래프를 생성합니다. 이를 위해, 값이 정확히 정수가 아닌 경우에는 해를 반올림한 후 결과 값을 logical형으로 변환합니다.

tspsol.trips = logical(round(tspsol.trips));
Gsol = graph(idxs(tspsol.trips,1),idxs(tspsol.trips,2),[],numnodes(G));
% Gsol = graph(idxs(tspsol.trips,1),idxs(tspsol.trips,2)); % Also works in most cases

기존 플롯에 새 그래프를 중첩하고 간선을 강조 표시합니다.

hold on
highlight(hGraph,Gsol,'LineStyle','-')
title('Solution with Subtours')

지도에서 알 수 있듯이 해에는 여러 하위 경로가 있습니다. 지금까지 지정한 제약 조건은 이러한 하위 경로의 발생을 막지 않습니다. 어떤 하위 경로의 발생도 허용하지 않으려면 엄청나게 많은 수의 부등식 제약 조건이 필요할 것입니다.

하위 경로 제약 조건

모든 하위 경로 제약 조건을 추가할 수는 없으므로 반복법을 사용하십시오. 현재 해에서 하위 경로를 검출한 후 부등식 제약 조건을 추가하여 특정 하위 경로가 발생하지 않도록 합니다. 이렇게 하면 몇 번의 반복으로 적합한 경로를 찾을 수 있습니다.

부등식 제약 조건을 사용하여 하위 경로를 제거합니다. 이에 대한 예로, 하위 경로에 5개의 점이 있다고 하면 이 점을 연결하는 5개의 직선이 하위 경로를 생성하게 됩니다. 이 5개의 점 사이의 직선의 개수가 4개보다 작거나 같아야 함을 나타내는 부등식 제약 조건을 만들어 이 하위 경로를 제거하십시오.

이 5개 점 사이의 모든 직선을 찾은 후 존재하는 직선의 개수가 4개를 넘지 않도록 해에 제약 조건을 적용합니다. 해에 5개 이상의 직선이 존재하면 해에 하위 경로가 생기기 때문에 이는 올바른 제약 조건입니다(n개 노드와 n개 간선이 포함된 그래프는 항상 순환을 포함함).

현재 해에서 간선을 사용하여 생성한 그래프인 Gsol의 연결성분을 식별하여 하위 경로를 검출합니다. conncomp는 각 간선이 속한 하위 경로 개수를 가진 벡터를 반환합니다.

tourIdxs = conncomp(Gsol);
numtours = max(tourIdxs); % Number of subtours
fprintf('# of subtours: %d\n',numtours);
# of subtours: 27

선형 부등식 제약 조건을 포함시켜 하위 경로를 제거하고 단 하나의 하위 경로만 남을 때까지 솔버를 반복해서 호출합니다.

% Index of added constraints for subtours
k = 1;
while numtours > 1 % Repeat until there is just one subtour
    % Add the subtour constraints
    for ii = 1:numtours
        inSubTour = (tourIdxs == ii); % Edges in current subtour
        a = all(inSubTour(idxs),2); % Complete graph indices with both ends in subtour
        constrname = "subtourconstr" + num2str(k);
        tsp.Constraints.(constrname) = sum(trips(a)) <= (nnz(inSubTour) - 1);
        k = k + 1;        
    end
    
    % Try to optimize again
    [tspsol,fval,exitflag,output] = solve(tsp,'options',opts);
    tspsol.trips = logical(round(tspsol.trips));
    Gsol = graph(idxs(tspsol.trips,1),idxs(tspsol.trips,2),[],numnodes(G));
    % Gsol = graph(idxs(tspsol.trips,1),idxs(tspsol.trips,2)); % Also works in most cases
    
    % Plot new solution
    hGraph.LineStyle = 'none'; % Remove the previous highlighted path
    highlight(hGraph,Gsol,'LineStyle','-')
    drawnow

    % How many subtours this time?
    tourIdxs = conncomp(Gsol);
    numtours = max(tourIdxs); % Number of subtours
    fprintf('# of subtours: %d\n',numtours)    
end
# of subtours: 20
# of subtours: 7
# of subtours: 9
# of subtours: 9
# of subtours: 3
# of subtours: 2
# of subtours: 7
# of subtours: 2
# of subtours: 1
title('Solution with Subtours Eliminated');
hold off

해의 품질

해는 하나의 폐루프이므로 하나의 실현 가능한 경로를 나타냅니다. 하지만 이것이 최소 비용 경로일까요? 이를 확인하는 한 가지 방법은 출력 구조체를 검사하는 것입니다.

disp(output.absolutegap)
   1.8300e-04

절대 간격이 작다는 사실은 해가 최적이거나 해의 총 길이가 최적에 가깝다는 것을 의미합니다.

관련 항목