ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

巡回セールスマン問題

この例では、0-1 整数計画法を使用して古典的な巡回セールスマンの問題を解く方法を説明します。この問題では、一連の訪問先 (都市) を通過して最短の閉じられたツアー (パス) を見つけます。この場合、200 の訪問先がありますが、nStops 変数を容易に変更して異なるサイズの問題を取得できます。最初の問題を解くと、解にサブツアーがあることがわかります。これは、見つけた最適解では訪問先をすべて通過する 1 つの連続したパスが得られず、代わりに複数の切断されたループが存在することを示します。その後、反復処理を使用してサブツアーの特定、制約の追加を行い、サブツアーが消去されるまで最適化を再実行します。

地図と訪問先の描画

アメリカ合衆国本土の粗い多角形表現内にランダムな訪問先を生成します。

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,xx,yy) % 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')
hold off

問題の定式化

整数線形計画法の巡回セールスマンの問題を次のように定式化します。

  • 可能なすべての経路 (つまり、訪問先のすべての明確なペア) を生成します。

  • 各経路の距離を計算します。

  • 最小化するコスト関数はツアー内の各経路の経路間の距離の合計です。

  • 決定変数は 2 進数であり、各経路に関連付けられています。ここで、1 はツアー上に存在する経路を表し、各 0 はツアー上に存在しない経路を表します。

  • 訪問先がすべてツアーに含まれていることを確認するには、各訪問先が 2 つの経路上に必ず配置されるように線形制約を含めます。これは訪問先に 1 回到着し、訪問先から 1 回出発することを意味します。

点間の距離の計算

訪問先が 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 は 2 進数の解ベクトルです。これが最小化するツアーの距離です。

等式制約

問題には 2 種類の等式制約があります。1 番目の制約は合計 200 の経路が適用されなければならないことです。2 番目の制約は、各訪問先に 2 つの経路が結合されていなければならないことです (つまり、各訪問先に到着する経路と各訪問先から出発する経路が存在しなければならないことです)。

最初の種類の等式制約を指定します。nStops 経路を Aeq*x_tsp = beq 形式で指定しなければなりません。

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

2 番目の種類の等式制約を指定するには、各訪問先に 2 つの経路を結合し、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)];

2 進数範囲

決定変数はすべて 2 進数です。ここで、intcon 引数を決定変数の数に設定し、それぞれ下限を 0、上限を 1 に設定します。

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

intlinprog を使用した最適化

問題を解く準備が整いました。ソルバーを呼び出します。

opts = optimoptions('intlinprog','Display','off');
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);

解の可視化

hold on
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 つの点の間に存在する線は 4 つ以下でなければならないという不等式制約を実装することによってサブツアーを消去します。

さらには、これらの 5 つの点の間に存在するすべての線を見つけ、これらの線が 4 つを超えて存在しないように解を制約します。解に 5 つ以上の線が存在している場合、解にはサブツアーが含まれる ( $n$ ノードと $n$ エッジをもつグラフには必ず円が含まれる) ので、これは正しい制約です。

関数 detectSubtours は解を解析し、ベクトルのセル配列を返します。セル配列内の各ベクトルは、その特定のサブツアーに含まれる訪問先を含みます。

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

線形不等式制約を含めてサブツアーを消去し、1 つのサブツアーのみが残るまでソルバーを繰り返し呼び出します。

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

title('Solution with Subtours Eliminated');
hold off
# of subtours: 16
# of subtours: 10
# of subtours: 5
# of subtours: 3
# of subtours: 4
# of subtours: 1

解の質

解は単一の閉ループなので、解は実行可能なツアーを表します。しかし、これは最小コストのツアーでしょうか。これを確認する 1 つの方法は、出力構造体を調べることです。

disp(output.absolutegap)
    0.0021

ギャップの絶対値が小さいことは解が最適であるか、ツアーの長さの合計が最適に近いことを意味します。

この情報は役に立ちましたか?