巡回セールスマン問題: 問題ベース

ライブスクリプトを開く

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

この問題に対するソルバーベースのアプローチについては、巡回セールスマン問題: ソルバーベースを参照してください。

問題の定式化

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

可能なすべての経路 (つまり、訪問先のすべての明確なペア) を生成します。
各経路の距離を計算します。
最小化するコスト関数はツアー内の各経路の経路間の距離の合計です。
決定変数は 2 進数であり、各経路に関連付けられています。ここで、1 はツアー上に存在する経路を表し、各 0 はツアー上に存在しない経路を表します。
訪問先がすべてツアーに含まれていることを確認するには、各訪問先が 2 つの経路上に必ず配置されるように線形制約を含めます。これは訪問先に 1 回到着し、訪問先から 1 回出発することを意味します。

訪問先の生成

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

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 の 2 値変数が存在します (# 変数 = 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 は解の経路を表す 2 値ベクトルです。これが最小化するツアーの距離です。

グラフの作成と地図の描画

問題をグラフとして表します。訪問先をノード、経路をエッジとしてグラフを作成します。

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 つの経路が結合されているという線形制約を作成します。つまり、各訪問先に到着する経路と各訪問先から出発する経路が存在しなければならないということです。

グラフ表現を使用して、訪問先と接続されているすべてのエッジを見つけることにより、訪問先で開始または終了する経路をすべて特定します。各訪問先について、その訪問先の経路の合計は 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 つの点の間に存在する線は 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

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

% 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

解の質

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

disp(output.absolutegap)

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