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

ライブスクリプトを開く

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

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

問題の定式化

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

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

訪問先の生成

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

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

点間の距離の計算

訪問先が 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'*x_tsp

ここで、x_tsp は 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

制約

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

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,:) = whichIdxs'; % Include in the constraint matrix
end
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);

解の経路をエッジとして使用して、新しいグラフを作成します。これを行うために、厳密に整数でない値が含まれる場合に備えて解を丸め、結果の値を logical に変換します。

x_tsp = logical(round(x_tsp));
Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));
% Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,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 つのサブツアーのみが残るまでソルバーを繰り返し呼び出します。

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 = find(tourIdxs == 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);
    x_tsp = logical(round(x_tsp));
    Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));
    % Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2)); % Also works in most cases
    
    % Visualize result
    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)

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