Формулировка задачи

Сформулируйте задачу коммивояжера для целочисленного линейного программирования следующим образом:

Сгенерируйте упоры

Сгенерируйте случайные упоры внутри неочищенного полигонального представления континентальной части США

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 двоичных переменных (# переменные = 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 - вектор двоичного решения. Это расстояние тура, которое вы стараетесь минимизировать.

Создайте график и рисуйте карту

Представьте задачу как график. Создайте график, где упоры являются узлами, а отключения - ребрами.

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

Ограничения

Создайте линейные ограничения, чтобы каждая остановка имела два связанных отключения, потому что должно быть отключение каждой остановки и отключение, отходящее от каждой остановки.

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);

Двоичные границы

Все переменные принятия решений являются двоичными. Теперь установите 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')

Как видно на карте, решение имеет несколько субтуров. Указанные до настоящего времени ограничения не препятствуют возникновению этих субтур. В порядок, чтобы предотвратить любое возможное подтекст, вам понадобится невероятно большое количество ограничений неравенства.

Субтактные ограничения

Поскольку вы не можете добавить все субтактные ограничения, примите итеративный подход. Обнаружите субтуры в текущем решении, затем добавьте ограничения неравенства, чтобы предотвратить эти конкретные субтуры. Делая это, вы находите подходящий тур в нескольких итерациях.

Исключить подтуры с ограничениями неравенства. Пример того, как это работает, если у вас есть пять точек в подтексте, тогда у вас есть пять линий, соединяющих эти точки, чтобы создать подтекст. Устраните этот подтур, реализуя ограничение неравенства, чтобы сказать, что между этими пятью точками должно быть меньше или равно четырем линиям.

Еще больше, найти все линии между этими пятью точками и ограничить решение, чтобы не иметь более четырех из этих линий. Это является правильным ограничением, потому что, если бы пять или более линий существовали в решении, то решение имело бы подтату (a графика с $$n$$ узлы и $$n$$ ребра всегда содержат цикл).

Обнаружение субтругов путем идентификации связанных компонентов в Gsol, графиков, построенный с ребрами в текущем решении. conncomp возвращает вектор с номером подтура, которому принадлежит каждое ребро .

tourIdxs = conncomp(Gsol);
numtours = max(tourIdxs); % 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 = 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

Качество решения

Решение представляет допустимое турне, потому что это один замкнутый цикл. Но это тур с минимальными затратами? Один из способов выяснить это - изучить структуру output.

disp(output.absolutegap)

     0

Малость абсолютной погрешности подразумевает, что решение либо оптимально, либо имеет общую длину, которая близка к оптимальной.