ドキュメンテーション

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

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

この例では、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,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')
hold off

問題の定式化

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

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

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

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

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

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

点間の距離の計算

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

変数および問題の作成

問題および 2 値変数を作成します。

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

目的関数を問題に含めます。

tsp.Objective = dist'*trips;

等式制約

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

最初の種類の等式制約を指定します。nStops 経路を指定し、問題に含めなければなりません。

constrips = sum(trips) == nStops;
tsp.Constraints.constrips = constrips;

2 番目の種類の等式制約を指定するには、各訪問先に 2 つの経路を結合し、訪問先ごとの経路を見つけ、その訪問先の経路数を合計する必要があります。その訪問先が開始点と終了点の両方である経路を見てみます。

constr2trips = optimconstr(nStops,1);
for stops = 1:nStops
    whichIdxs = (idxs == stops);
    whichIdxs = any(whichIdxs,2); % start or end at stops
    constr2trips(stops) = sum(trips(whichIdxs)) == 2;
end
tsp.Constraints.constr2trips = constr2trips;

初期問題を解く

問題を解く準備が整いました。解を迅速に得るために役立つ 'round-diving' ヒューリスティックを使用してソルバーを呼び出し、整数前処理をオフにしてさらに高速化します。

opts = optimoptions('intlinprog','Display','off','Heuristics','round-diving',...
    'IPPreprocess','none');
tspsol = solve(tsp,'options',opts)
tspsol = struct with fields:
    trips: [19900×1 double]

解の可視化

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

地図上で確認できるように、解には複数のサブツアーが存在します。これまでに指定された制約では、これらのサブツアーの発生を回避することはできません。サブツアーの発生を回避するには、途方もない数の不等式制約が必要になります。

サブツアーの制約

サブツアーの制約をすべて追加することはできないため、反復法を使用します。現在の解にあるサブツアーを検出し、不等式制約を追加してこれらの特定のサブツアーの発生を回避します。このようにすることで、数回の反復で適切なツアーを見つけることができます。

不等式制約によってサブツアーを消去します。サブツアー内に 5 つの点がある場合、これらの点を接続してサブツアーを作成し、これが実際にどのように機能するかの例を示します。これらの 5 つの点の間に存在する線は 4 つ以下でなければならないという不等式制約を実装することによってサブツアーを消去します。

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

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

tours = detectSubtours(tspsol.trips,idxs);
numtours = length(tours); % 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
        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);
        a = false(length(idxs),1);
        for jj = 1:length(variations)
            whichVar = (sum(idxs==subTourIdx(variations(jj,1)),2)) & ...
                       (sum(idxs==subTourIdx(variations(jj,2)),2));
            a = a | whichVar;
        end
        tsp.Constraints.(sprintf('subtourconstr%i',k)) = sum(trips(a)) <= length(subTourIdx)-1;
        k = k + 1;
    end
    % Try to optimize again
    [tspsol,fval,exitflag,output] = solve(tsp,'options',opts);

    % Visualize result
    lh = updateSalesmanPlot(lh,tspsol.trips,idxs,stopsLon,stopsLat);

    % How many subtours this time?
    tours = detectSubtours(tspsol.trips,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

解の質

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

disp(output.absolutegap)
     0

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