이 페이지의 최신 내용은 아직 번역되지 않았습니다. 최신 내용은 영문으로 볼 수 있습니다.

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

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

문제 기반 접근법에 대해서는 여행하는 외판원 문제: 문제 기반 항목을 참조하십시오.

지도와 경유지 그리기

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

figure;

load('usborder.mat','x','y','xx','yy');
rng(3,'twister') % makes a plot with 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
plot(x,y,'Color','red'); % draw the outside border
hold on
% Add the stops to the map
plot(stopsLon,stopsLat,'*b')

문제 정식화

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

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

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

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

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

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

점 간의 거리 계산하기

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'*x_tsp

여기서 x_tsp는 이진 해 벡터입니다. 이 경로의 길이가 바로 최소화하고자 하는 경로의 거리입니다.

등식 제약 조건

이 문제에는 두 가지 유형의 등식 제약 조건이 있습니다. 첫 번째는 총 여행 횟수가 200이 되어야 한다는 것입니다. 두 번째는 각 경유지마다 두 개의 여행이 있어야 한다는 것입니다(즉, 각 경유지마다 도착하는 여행과 출발하는 여행이 각각 하나씩 있어야 함).

Aeq*x_tsp = beq 형식으로 nStops개의 여행이 있어야 한다는 것을 규정하는 첫 번째 등식 제약 조건을 지정합니다.

Aeq = spones(1:length(idxs)); % Adds up the number of trips
beq = nStops;

각 경유지마다 두 개의 여행이 연결되어 있어야 한다는 것을 규정하는 두 번째 등식 제약 조건을 지정하려면 Aeq 행렬을 희소 행렬로 확장하십시오.

Aeq = [Aeq;spalloc(nStops,length(idxs),nStops*(nStops-1))]; % allocate a sparse matrix
for ii = 1:nStops
    whichIdxs = (idxs == ii); % find the trips that include stop ii
    whichIdxs = sparse(sum(whichIdxs,2)); % include trips where ii is at either end
    Aeq(ii+1,:) = whichIdxs'; % include in the constraint matrix
end
beq = [beq; 2*ones(nStops,1)];

이진 범위

모든 결정 변수가 이진 변수입니다. 이제, intcon 인수를 결정 변수의 개수로 설정하고, 각 변수의 하한과 상한을 각각 0과 1로 지정합니다.

intcon = 1:lendist;
lb = zeros(lendist,1);
ub = ones(lendist,1);

intlinprog를 사용하여 최적화하기

이제 문제를 풀 수 있습니다. 해를 신속하게 구할 수 있도록 'round-diving' 발견법(Heuristics)을 사용하여 솔버를 호출하고 속도가 더욱 빨라지도록 정수 전처리를 끕니다.

opts = optimoptions('intlinprog','Display','off','Heuristics','round-diving',...
    'IPPreprocess','none');
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);

해 시각화하기

segments = find(x_tsp); % Get indices of lines on optimal path
lh = zeros(nStops,1); % Use to store handles to lines on plot
lh = updateSalesmanPlot(lh,x_tsp,idxs,stopsLon,stopsLat);
title('Solution with Subtours');

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

하위 경로 제약 조건

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

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

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

detectSubtours 함수는 해를 분석한 후 벡터로 구성된 셀형 배열을 반환합니다. 셀형 배열에 있는 각 벡터는 해당 특정 하위 경로에 포함된 경유지를 포함합니다.

tours = detectSubtours(x_tsp,idxs);
numtours = length(tours); % number of subtours
fprintf('# of subtours: %d\n',numtours);
# of subtours: 27

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

A = spalloc(0,lendist,0); % Allocate a sparse linear inequality constraint matrix
b = [];
while numtours > 1 % repeat until there is just one subtour
    % Add the subtour constraints
    b = [b;zeros(numtours,1)]; % allocate b
    A = [A;spalloc(numtours,lendist,nStops)]; % a guess at how many nonzeros to allocate
    for ii = 1:numtours
        rowIdx = size(A,1)+1; % Counter for indexing
        subTourIdx = tours{ii}; % Extract the current subtour
%         The next lines find all of the variables associated with the
%         particular subtour, then add an inequality constraint to prohibit
%         that subtour and all subtours that use those stops.
        variations = nchoosek(1:length(subTourIdx),2);
        for jj = 1:length(variations)
            whichVar = (sum(idxs==subTourIdx(variations(jj,1)),2)) & ...
                       (sum(idxs==subTourIdx(variations(jj,2)),2));
            A(rowIdx,whichVar) = 1;
        end
        b(rowIdx) = length(subTourIdx)-1; % One less trip than subtour stops
    end

    % Try to optimize again
    [x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,A,b,Aeq,beq,lb,ub,opts);
    
    % Visualize result
    lh = updateSalesmanPlot(lh,x_tsp,idxs,stopsLon,stopsLat);
    
    % How many subtours this time?
    tours = detectSubtours(x_tsp,idxs);
    numtours = length(tours); % 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)
     0

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

관련 항목