В этом примере показано, как использовать бинарное целочисленное программирование, чтобы решить классическую задачу коммивояжера. Эта проблема включает нахождение самого короткого закрытого тура (путь) через набор остановок (города). В этом случае существует 200 остановок, но можно легко изменить nStops переменная, чтобы получить различный проблемный размер. Вы будете решать начальную задачу и видеть, что решение имеет подтуры. Это означает, что найденное оптимальное решение не дает один непрерывный путь через все точки, но вместо этого имеет несколько разъединенных циклов. Вы будете затем использовать итеративный процесс определения подтуров, добавления ограничений и повторного выполнения оптимизации, пока подтуры не будут устранены.
Для основанного на решателе подхода к этой проблеме смотрите проблему Коммивояжера: основанный на решателе.
Сформулируйте проблему коммивояжера для целочисленного линейного программирования можно следующим образом:
Сгенерируйте все возможные прохождения, означая все отличные пары остановок.
Вычислите расстояние для каждого прохождения.
Функция стоимости, чтобы минимизировать является суммой расстояний прохождения для каждого прохождения на туре.
Переменные решения являются двоичным файлом, и сопоставленный с каждым прохождением, где каждый 1 представляет прохождение, которое существует в туре, и каждый 0 представляет прохождение, которое не находится в туре.
Чтобы гарантировать, что тур включает каждую остановку, включайте линейное ограничение, что каждая остановка включена точно два прохождения. Это означает одно прибытие и одно отклонение от остановки.
Сгенерируйте случайные остановки в грубом многоугольном представлении континентальных США.
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));
Наложите новый график на существующем графике и подсветите его ребра.
hold on highlight(hGraph,Gsol,'LineStyle','-') title('Solution with Subtours')
Как видно на карте, решение имеет несколько подтуров. Ограничения, заданные до сих пор, не предотвращают эти подтуры. Для того, чтобы предотвратить любой возможный подтур, вам было бы нужно невероятно большое количество ограничений неравенства.
Поскольку вы не можете добавить все подтуристические ограничения, проявите итерационный подход. Обнаружьте подтуры в текущем решении, затем добавьте ограничения неравенства, чтобы предотвратить те конкретные подтуры. Путем выполнения этого вы находите подходящий тур в нескольких итерациях.
Устраните подтуры по-разному ограничения. Пример того, как это работы - то, если у вас есть пять точек на подтуре, затем у вас есть пять линий, соединяющих те точки, чтобы создать подтур. Устраните этот подтур путем реализации ограничения неравенства, чтобы сказать, что там должно быть меньше чем или равно четырем линиям между этими пятью точками.
Еще больше найдите все линии между этими пятью точками и ограничьте решение не иметь больше чем четыре из этих существующих линий. Это - правильное ограничение потому что, если бы пять или больше из линий существовали в решении, то решение имело бы подтур (график с узлы и ребра всегда содержат цикл).
Обнаружьте подтуры путем идентификации связанных компонентов в Gsol, график создается с ребрами в текущем решении. conncomp возвращает вектор с количеством подтура к который каждое ребро belongs.
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)); % 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
Малость абсолютного разрыва подразумевает, что решение или оптимально или имеет общую длину, которая является близко к оптимальному.