Main Content

Problema del viajante: basado en solvers

Abrir script en vivo

Este ejemplo muestra cómo utilizar programación de enteros binaria para resolver el problema clásico del viajante. Este problema implica encontrar el trayecto (ruta) cerrado más corto pasando por un conjunto de paradas (ciudades). En este caso hay 200 paradas, pero puede cambiar fácilmente la variable nStops para obtener un tamaño de problema diferente. Resolverá el problema inicial y verá que la solución tiene subtrayectos. Esto significa que la solución óptima encontrada no proporciona una ruta continua por todos los puntos, sino que tiene varios bucles desconectados. Después utilizará un proceso iterativo para determinar los subtrayectos, añadiendo restricciones y volviendo a ejecutar la optimización hasta que los subtrayectos se eliminen.

Para el enfoque basado en problemas, consulte Traveling Salesman Problem: Problem-Based.

Formulación del problema

Formule el problema del viajante para programación lineal de enteros de la siguiente forma:

  • Genere todos los viajes posibles, es decir, todas las parejas de paradas diferentes.

  • Calcule la distancia para cada viaje.

  • La función de coste que se desea minimizar es la suma de las distancias del viaje para cada viaje del trayecto.

  • Las variables de decisión son binarias y están asociadas a cada viaje, donde cada 1 representa un viaje que existe en el trayecto y cada 0 representa un viaje que no está en el trayecto.

  • Para garantizar que el trayecto incluya todas las paradas, añada la restricción lineal de que cada parada esté exactamente entre dos viajes. Esto implica una llegada y una salida en cada parada.

Generar paradas

Genere paradas aleatorias dentro de una representación poligonal en bruto de los EE. UU. continentales.

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

Calcular distancias entre puntos

Dado que existen 200 paradas, hay 19.900 viajes, es decir, 19.900 variables binarias (# variables = 200 elegir 2).

Genere todos los viajes, es decir, todas las parejas de paradas.

idxs = nchoosek(1:nStops,2);

Calcule todas las distancias de los viajes, asumiendo que la tierra es plana para utilizar el teorema de Pitágoras.

dist = hypot(stopsLat(idxs(:,1)) - stopsLat(idxs(:,2)), ...
             stopsLon(idxs(:,1)) - stopsLon(idxs(:,2)));
lendist = length(dist);

Con esta definición del vector dist, la longitud de un trayecto es

dist'*x_tsp

donde x_tsp es el vector de solución binario. Esta es la distancia de un trayecto que se intenta minimizar.

Crear una gráfica y diseñar un mapa

Represente el problema como una gráfica. Cree una gráfica donde todas las paradas sean nodos y los viajes sean bordes.

G = graph(idxs(:,1),idxs(:,2));

Muestre las paradas utilizando una gráfica. Represente los nodos sin los bordes de la gráfica.

figure
hGraph = plot(G,'XData',stopsLon,'YData',stopsLat,'LineStyle','none','NodeLabel',{});
hold on
% Draw the outside border
plot(x,y,'r-')
hold off

Restricciones

Cree las restricciones lineales de que cada parada tenga dos viajes asociados, ya que debe haber un viaje hasta cada parada y un viaje desde cada parada.

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

Límites binarios

Todas las variables de decisión son binarias. Establezca el argumento intcon en el número de variables de decisión y ponga un límite inferior de 0 y un límite superior de 1 en cada una.

intcon = 1:lendist;
lb = zeros(lendist,1);
ub = ones(lendist,1);

Optimizar utilizando intlinprog

El problema está listo para su resolución. Para suprimir la salida iterativa, desactive la visualización predeterminada.

opts = optimoptions('intlinprog','Display','off');
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);

Cree una nueva gráfica con los viajes de la solución como bordes. Para hacerlo, redondee la solución en caso de que algunos valores no sean exactamente enteros, y convierta los valores resultantes a 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

Visualizar la solución

hold on
highlight(hGraph,Gsol,'LineStyle','-')
title('Solution with Subtours')

Como puede verse en el mapa, la solución tiene varios subtrayectos. Las restricciones especificadas hasta ahora no evitan estos subtrayectos. Para evitar cualquier posible subtrayecto, necesitaría un número increíblemente amplio de restricciones de desigualdad.

Restricciones de subtrayecto

Dado que no se pueden añadir todas las restricciones de subtrayecto, adopte un enfoque iterativo. Detecte los subtrayectos en la solución actual y, después, añada restricciones de desigualdad para evitar esos subtrayectos en particular. Al hacerlo, se encuentra un trayecto adecuado en unas pocas iteraciones.

Elimine los subtrayectos con restricciones de desigualdad. Un ejemplo de cómo funciona es que, si tiene cinco puntos en un subtrayecto, tendrá cinco líneas conectando estos puntos para crear el subtrayecto. Elimine este subtrayecto implementando una restricción de desigualdad que diga que debe haber cuatro líneas o menos entre estos cinco puntos.

Aún mejor, encuentre todas las líneas entre estos cinco puntos y restrinja la solución para que no tenga más de cuatro de estas líneas. Se trata de una restricción correcta porque, si en una solución existen cinco o más líneas, la solución tendría un subtrayecto (una gráfica con n nodos y n bordes siempre contiene un ciclo).

Detecte el subtrayecto identificando los componentes conectados en Gsol, la gráfica creada con los bordes en la solución actual. conncomp devuelve un vector con el número del subtrayecto al que pertenece cada borde.

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

Incluya las restricciones de desigualdad lineales para eliminar subtrayectos y llame repetidamente al solver hasta que solo quede un subtrayecto.

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

Calidad de la solución

La solución representa un trayecto factible porque se trata de un bucle cerrado único. Pero, ¿se trata de un trayecto con coste mínimo? Una forma de descubrirlo es examinar la estructura de salida.

disp(output.absolutegap)
     0

El tamaño pequeño del intervalo absoluto implica que la solución es óptima o tiene una longitud total casi óptima.

Temas relacionados