Problema del viajante: Basado en problemas
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.
Si desea ver el enfoque basado en solvers para este problema, consulte Problema del viajante: Basado en solvers.
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 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'*trips
donde trips
es el vector binario que representa los viajes que realiza la solución. 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
Crear variables y problemas
Cree un problema de optimización con variables de optimización binarias que representan los viajes potenciales.
tsp = optimproblem; trips = optimvar('trips',lendist,1,'Type','integer','LowerBound',0,'UpperBound',1);
Incluya la función objetivo en el problema.
tsp.Objective = dist'*trips;
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.
Utilice la representación de la gráfica para identificar todos los viajes que comienzan y terminan en una parada encontrando todos los bordes que están conectados a esa parada. Para cada parada, cree la restricción de que la suma de los viajes para esa parada es igual a dos.
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;
Resolver el problema inicial
El problema está listo para resolverlo. Para suprimir la salida iterativa, desactive la visualización predeterminada.
opts = optimoptions('intlinprog','Display','off'); tspsol = solve(tsp,'options',opts)
tspsol = struct with fields:
trips: [19900×1 double]
Visualizar la solución
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
.
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
Superponga la nueva gráfica sobre la gráfica existente y destaque sus bordes.
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 nodos y 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.
% 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
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)
1.8300e-04
El tamaño pequeño del intervalo absoluto implica que la solución es óptima o tiene una longitud total casi óptima.