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

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

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

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

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) % tTest 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'*trips

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

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

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

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;

Ограничения

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

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

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

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

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

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

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

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

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

tourIdxs = conncomp(Gsol);
numtours = max(tourIdxs); % Number of subtours
fprintf('# of subtours: %d\n',numtours);

# of subtours: 27

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

% 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

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

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

disp(output.absolutegap)

     0

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