Facultad

de
Ciencias

OPTIMIZACIÓN DEL PROBLEMA DEL

VIAJANTE: MÉTODOS DE BRANCH AND
BOUND Y LIN KERNIGHAN
(Optimisation of the Travelling Salesman Problem:
Branch and Bound and Lin-Kernighan Methods)

Trabajo de Fin de Grado

para acceder al

GRADO EN MATEMÁTICAS

Autora: Marı́a Pino Rey

Directora: Mónica Blanco Gómez
Septiembre - 2024
Resumen
El problema del Viajante (TSP) es formulado como un problema de programación lineal entera
que trata de encontrar un ciclo hamiltoniano que minimice los costes de las aristas utilizadas sobre
un grafo completo ponderado G = (V, A).

Al tratarse de un problema NP-completo tiene diversas técnicas de resolución tanto exactas

como heurı́sticas. Se planteará el método exacto de Branch and Bound con la regla de ramificación
Little y el método heurı́stico de mejora iterativa k − opt Lin Kernighan, además de exponer pequeños
problemas ilustrativos de ejemplo para ver el procedimiento de los métodos.

Finalmente se implementarán en Matlab los algoritmos mencionados para la resolución de dos

problemas reales comparando las soluciones y tiempos de ejecución.

Palabras clave: Problema del viajante, grafo, programación lineal, ramificación, intercambio

Abstract
The Traveling Salesman Problem (TSP) is formulated as an integer linear programming problem
what tries to find a Hamiltonian cycle that minimises the costs of the edges used on a weighted
graph G = (V, A).

As it is an NP-complete problem, it has several exact and heuristic solution techniques. The
Branch and Bound exact method with the Little branching rule and the Lin Kernighan heuristic
iterative improvement method will be presented, as well as small illustrative example problems to
show the procedure of the methods.

Finally, the aforementioned algorithms will be implemented in Matlab to solve two real problems,
comparing the solutions and execution times.

Key words: Travelling Salesman Problem, graph , linear programming, branching, exchange
Agradecimientos

Agradecer a mi tutora, Mónica, por aceptar guiarme en este proyecto, ayudarme cuando lo
necesité, toda su dedicación y todo su tiempo empleado.

A todos mis profesores, si no fuera por ellos no serı́a posible la realización de este trabajo.

A mis amigos, los que ya lo eran y los que tuve la suerte de que ahora lo sean. Por estar siempre
ahı́ y confı́ar en mı́ incluso cuando ni yo misma lo hacı́a, celebrando desde la distancia mis logros
como si fueran los suyos.
A Sergi, por hacerme sin saberlo, los malos momentos mucho más llevaderos.

A mi familia, que me dieron ánimo y fuerza para que nunca me diese por vencida. En especial a
mi tı́a, por estar siempre más que pendiente de su “tercera hija”.

A mis abuelos, que seguro que estarán muy orgullosos de mi.

¡Abu, lo conseguı́!

A mi madre, mi ejemplo a seguir, por ser mi gran apoyo incondicional y enseñarme que con
trabajo y esfuerzo conseguiré todo lo que me proponga.

Y por último como todo lo que hago, va dedicado a ti, Papá, por no dejar que nunca me rinda y
acompañarme siempre en cada paso que soy capaz de dar.
Índice general

1. Introducción 1
1.1. Preliminares de Teorı́a de Grafos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2. Preliminares de Programación lineal . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2. El problema 4
2.1. Planteamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2. Historia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3. Formulación TSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.4. Complejidad Temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.5. Técnicas de Resolución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.5.1. Vecino más próximo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3. Ramificación y acotación (Branch and Bound) 13

3.1. Planteamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.1. Relajación del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.2. Regla de Ramificación Little . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1.3. Exploración del árbol de decisión . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2. Implementación en Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3. Ejemplo ilustrativo para sTSP con Branch and Bound . . . . . . . . . . . . . . . . 21

4. Lin-Kernighan 24
4.1. Planteamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.2. Comentarios sobre el algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
4.3. Implementación en Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.4. Ejemplo ilustrativo para sTSP con Lin-Kernighan . . . . . . . . . . . . . . . . . . 33
4.5. Pruebas en Matlab con implementación de Lin-Kernighan . . . . . . . . . . . . . . 36

5. Aplicación y conclusiones 38
5.1. Ejemplo práctico TSP Euclı́deo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.1.1. Solver intlinprog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.1.2. Branch and Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.1.3. Vecino más cercano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.1.4. Lin-Kernighan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.2. Ejemplo práctico TSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.2.1. Solver intlinprog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.2.2. Branch and Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2.3. Lin-Kernighan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

A. Programas en Matlab 47
A. Vecino más cercano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
B. Branch and Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
B.1. BB TSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
B.2. BB base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
B.3. reducematriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
B.4. calcula penalización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3
 B.5. verificar exciclo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
C. Lin-Kernighan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
C.1. LK TSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
C.2. LK base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
C.3. encontrar t4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
C.4. encontrar t5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
C.5. aumentar camino . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
C.6. calcula coste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
C.7. es tour . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
C.8. soluc aristas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

B. Tablas de datos 64
A. Tabla de coordenadas de las 38 ciudades . . . . . . . . . . . . . . . . . . . . . . . . 64
B. Tabla con las distancias entre las aristas existentes de las 18 ciudades . . . . . . . . 65

4
Capı́tulo 1

Introducción

Este documento estudiará el Problema del Viajante, uno de los problemas más estudiados en el
campo de las Matemáticas, desde las perspectivas de la Teorı́a de Grafos y la Programación Lineal
Entera abordando su resolución mediante dos técnicas especı́ficas explicadas tras la recopilación
de diversas fuentes y autores. La motivación principal es desarollar una aplicación práctica que
demuestre su utilidad en situaciones cotidianas a través de una de sus múltiples aplicaciones como
es la logı́stica o el transporte.

La estructura de este documento es la siguiente:

En el este capı́tulo se expondrá una serie de preliminares para comprender de manera más clara
el contexto de lo que se tratará a continuación.

En el segundo capı́tulo, se planteará el problema más en profundidad ası́ como la explicación de

sus orı́genes y la formulación que será empleada en todo el documento.

El tercer capı́tulo, tratará de explicar el método exacto Branch and Bound ilustrandolo con un
breve ejemplo práctico además de la explicación de su implementación en Matlab. El cuarto, tomará
la misma estructura que el capı́tulo anterior pero en este caso será explicado el método heurı́stico
Lin Kernighan.

Por último, se desarrolla una de sus aplicaciones más conocida, el ámbito de la logı́stica/transporte.
Es planteada desde dos puntos de vista diferentes; teniendo en cuenta la propia red de carreteras
existente en España, concretamente en Galicia, y creando de manera artificial carreteras que
conectan de manera rectilı́nea todos los puntos escogidos del mapa. Para el análisis de los resultados
por los diferentes métodos planteatos en el documento como se dijo, se implementaron los métodos
en el programa Matlab.

1
 Se darán una serie de definiciones que serán de utilidad para la comprensión del documento.[7,
5, 6, 4, 20, 18]

1.1. Preliminares de Teorı́a de Grafos

Definición 1.1.1. Un grafo (finito) simple es una terna G = (V, A, φ) donde V = V (G) es
un conjunto finito no vacı́o (cuyos elementos se denominan vértices o nodos), A = A(G) es un
conjunto cuyos elementos se llaman aristas y φ : A −→ P (V ) es una aplicación que asocia a cada
arista un par de vértices. Es decir φ(x) = {i, j} con x ∈ A y i, j ∈ V . Se dice simple si no existen
lazos (i ̸= j ∀φ(x) = {i, j}) y dos vértices están unidos a lo sumo por una arista.
Definición 1.1.2. Un grafo es completo,denotado por Kn , si A = V × V ,es decir, para cualquier
2
pareja de vértices, existe una arista que los une. Siendo | A |= | V | .
Definición 1.1.3. Un grafo dirigido o digrafo es una terna G = (V, A, φ) donde V, A son como
en la definición de grafo y φ : A −→ V × V es una aplicación que asocia a cada arista un par de
vértices ordenados (una arista empieza en x y acaba en y).
Definición 1.1.4. Un grafo dirigido es simétrico si para toda arista (i, j) ∈ A también se cumple
que (j, i) ∈ A.
Definición 1.1.5. Un camino es una sucesión i0 x1 i1 x2 i2 · · · ik−1 xk ik donde is ∈ V y xs ∈ A una
arista incidente a is−1 y is donde todas las aristas son distintas y s ∈ {1, . . . k}.Se dice cerrado si
i0 = ik con is ∈ V . Se dice simple si no hay vértices repetidos salvo el primero y el último.
Definición 1.1.6. Un circuito es un camino cerrado.
Definición 1.1.7. Un ciclo es un camino cerrado simple no trivial (no es lazo).
Definición 1.1.8. Un grafo pesado o ponderado en aristas es un grafo G = (V, A) junto con
una aplicación w : A −→ R+ que asigna a cada arista un coste. El peso de un camino es la suma
de los pesos de las aristas que la componen.
Definición 1.1.9. Dados dos vértices i, j en un grafo pesado se llama distancia entre i e j, d(i, j),
al mı́nimo de los pesos de caminos que unen i con j.Un camino i0 e1 i1 · · · xk ik es de peso mı́nimo
si su peso es d(i0 , ik ).
Definición 1.1.10. Un subgrafo de un grafo G = (V, A) es un grafo H = (W, F ) donde W ⊆ V y
F ⊆A
Definición 1.1.11. Sea G = (V, A) un grafo simple .El subgrafo inducido por un subconjunto
W ⊂ V , denotado por G[W ], es el grafo H = (W, F ) donde el conjunto de aristas F contiene
una arista en A si y solo si ambos extremos de esta arista están en W. Es decir, si el subgrafo H
contiene todas las aristas que unen en G dos vértices de W.
Definición 1.1.12. Un árbol es un grafo conexo que no contiene ciclos y tiene n − 1 aristas.
Definición 1.1.13. Sea G un grafo y H ⊆ G un subgrafo. Se dice que H un árbol generador de
G si H es un árbol con V (H) = V (G).
Definición 1.1.14. Un árbol dirigido es un árbol que tiene un vértice distinguido al que se le
llama raı́z . Siendo a0 el nodo raı́z, se le llamará a a1 padre de a2 y a2 hijo de a1 .
Definición 1.1.15. Un árbol binario es aquel árbol que tiene exactamente dos hijos.
Definición 1.1.16. Un algoritmo de búsqueda de anchura sobre un árbol trata de construir un
árbol a partir de un vértice inicial i añadiendo los vértices sucesivos según su ditancia al vértice
original.
Definición 1.1.17. Un algoritmo de búsqueda en profundidad trata de construir un árbol a
partir de un vértice inicial i intentando volver atrás lo menos posible.

2
1.2. Preliminares de Programación lineal
Definición 1.2.1. Un Problema de Optimización (P) de dimensión finita puede definirse
como:
Dada f : Rn −→ R y K ⊆ Rn , el problema consistente en hallar x ∈ K tal que f (x) ≤ f (x)∀x ∈ K
formulado en la forma (
Minimizar f (x)
(P )
sujeto a x ∈ K
Se formula a partir de un conjunto de variables independientes,x, sobre las que pueden existir
condiciones (x ∈ K), llamadas restricciones del problema y una función dependiendo de las
variables que se quiere minimizar llamada función objetivo.

Definición 1.2.2. Un problema de optimización se considera un problema de Programación

Lineal(LP) si f es lineal, es decir, f (x) = cT x con c ∈ Rn y escrito es su forma estándar,
K = {x ∈ Rn : Ax = b, x ≥ 0} siendo A ∈ Rm×n y b ∈ Rm .
Definición 1.2.3. Un problema de Programación Lineal Entera Mixta(MILP) es un problema
lineal (LP) con algunas variables enteras.

Definición 1.2.4. Un problema de Programación Lineal Entera 0-1 es un problema MILP en

donde las variables unicamente toman valores 0 y 1.
Definición 1.2.5. Un problema lineal se dice relajado si se prescinde de alguna restricción del
problema original.

Definición 1.2.6. Un punto x que satisface todas las restricciones del problema, se considera
solución factible. El conjunto de todas las soluciones factibles se le llama región factible.
Definición 1.2.7. Una solución óptima ,x es el punto o puntos de la región factible que minimiza
la función objetivo, es decir, cuando se cumple que f (x) ≤ f (x) para toda x solución factible.

Definición 1.2.8. En un problema lineal 0-1 se le dice solución parcial a la solución en donde
algún valor está todavı́a sin fijar.
Definición 1.2.9. Se dice que x es un mı́nimo global de P si verfica f (x) ≤ f (x)∀x ∈ K
Definición 1.2.10. Se dice que x es un mı́nimo local de P si existe un entorno abierto de x, U,
tal que se verifica f (x) ≤ f (x)∀x ∈ K ∩ U

3
Capı́tulo 2

El problema

2.1. Planteamiento
El Problema del Viajante o también conocido como TSP por sus siglas del inglés Travelling
Salesman Problem, es un problema de optimización combinatoria cuyo número de posibles solu-
ciones crece exponencialmente con el tamaño del input. El problema,responde a la siguiente pregunta:

“Dada una lista de n ciudades y las distancias entre cada par de ciudades, ¿cuál es la ruta más
corta posible que visita cada ciudad exactamente una vez y vuelve a la ciudad de origen?”.
El TSP puede plantearse como un problema de Teorı́a de Grafos. Los nodos se corresponden
con las ciudades, las aristas con los caminos entre ciudades y los pesos de las aristas con los valores
de las distancias de ir de un lugar i a otro j representados como cij .

Partiendo de n lugares como una red de nodos, para cada arista (i, j) se tiene un valor cij
(Definición 1.1.8), tratando de encontrar en qué orden recorrer los nodos de la red de modo que
minimice la distancia total recorrida; dicho de otro modo, encontrar un circuito hamiltoniano
(Definición ??) con peso mı́nimo donde es indiferente el lugar que se designe como origen.

A su vez, puede ser planteado desde la Programación Lineal satisfaciendo lo que se denomina
función objetivo. Concretamente, atendiendo a la linealidad y naturaleza de las variables, se habla
de Programación Lineal Entera 0-1 (Definición 1.2.1) por tomar las variables xij , que se describirán
más adelante, valores 0 y 1.
Dentro de la Programación Lineal Entera, se encuentra la llamada Optimización Combinatoria,
donde hay un número finito aunque elevado de posibles decisiones. El TSP forma parte de esta
subclase de problemas.

Según su simetrı́a puede clasificarse en sTSP o TSP simétrico, donde es indiferente en qué
sentido recorrer las aristas, ya que tienen igual peso en ambos y en aTSP o TSP asimético donde sı́
importa el sentido de recorrido. De manera práctica, si una carretera es de doble sentido se hablará
de sTSP y si por la contra es de sentido único aTSP.

En general se hablará de TSP, ya que el sTSP es un caso particular del asimétrico.

2.2. Historia
El origen de este problema no está del todo claro, aunque se cree que se puede fechar por una
guı́a para vendendores ambulantes que incluye viajes de ejemplo por Alemania y Suiza, en el 1832.
Una versión de resolución mediante la Teorı́a de Grafos se puede datar en 1736 cuando Leonhard

4
Euler resolvió el Problema de los Siete Puentes de Königsberg o el Problema del caballo de ajedrez.

El primero de los problemas, el nombre, se debe a Königsberg, una ciudad de Prusia Oriental y
más tarde de Alemania que desde 1945 se convirtió en la ciudad rusa de Kaliningrado. Esta ciudad
está atravesada por el rı́o Pregolia que se bifurca y rodea con sus brazos a la isla Kneiphof quedando
ası́ dividido en cuatro regiones distintas que entonces estaban unidas mediante siete puentes.[21]

Figura 2.1: Grafo sobre mapa con los puentes de Königsberg [13]

Como se dijo, Euler resolvió el problema representando como se puede ver en la Figura 2.1 a la
ciudad de Königsberg en una especie de grafo donde los puentes a, b, c, d, e, f, g delimitados por dos
puntos (en rosa) determinan las cuatro partes del mapa {A, B, C, D} que se corresponderán con los
vértices. El objetivo era encontrar un recorrido para cruzar a pie toda la ciudad, pasando solo una
vez por cada uno de los puentes regresando al mismo punto de inicio. Años más tarde, él mismo,
concluyó que era imposible.

En el segundo de los problemas, se pedı́a que teniendo una cuadrı́cula de n × n casillas y un

caballo de ajedrez colocado en una posición cualquiera (a, b), el caballo pasase por todas las casillas
una sola vez con su movimiento en L y regresase a su punto de partida.
Este problema fue resuelto con éxito para el caso de un tablero de ajedrez completo (8 × 8) por
diversos expertos y métodos. Sin embargo, para otras dimensiones de tablero no es posible obtener
un resultado positivo. Es el caso de tableros en los que el número de casillas es impar o los tableros
rectangulares o en forma de cruz. Véase la demostración de la primera afirmación:
Ejemplo 2.2.1. Un tablero de ajedrez con un número impar de casillas no puede ser recorrido por
el caballo sin repetir ninguna casilla y volviendo al punto de partida.
Se resolverá atendiendo a la coloración de las casillas del tablero (blanco o negro). Suponga-
mos sin perder generalidad que el caballo comienza en una casilla blanca. Su movimiento en L
alterna tras cada movimiento el color de la casilla. Como se requiere que el punto final sea el de
partida según la suposición, se debe terminar en una casilla blanca. Esto es imposible de con-
seguir, ya que necesitarı́a realizar un número par de movimientos, pero el número de casillas es impar.

En 1857, se desarrolló el también conocido Icosian Game por Willian Rowan Hamilton que
consistı́a en un desafı́o matemático cuyo objetivo irradiaba en dar con el recorrido hamiltoniano por
las aristas de un dodecaedro para visitar una y solo una vez cada vértice coincidiendo el de llegada
con el de salida.
En la figura siguiente se representa una posible solución:

Desde otras perspectivas de resolución, en el siglo XIX el TSP fue formulado matemáticamente
por William Rowan Hamilton y Thomas Kirkman y desde entonces, ha sido estudiado por diversos
expertos en la materia. Tomó gran popularidad entre los años 50 y 60 debido a su estrecha relación
con los problemas de programación lineal de asignación y transporte. Ası́, en 1954 “Soluciones

5
 Figura 2.2: Posible solución al Icosian Game

de un problema del viajante de gran tamaño”de Dantzig, Fulkerson y Johnson en el Journal of

the Operations Research Society of America fue uno de los principales eventos en la historia de la
optimización[8]. El método utilizado parecı́a funcionar en pequeños problemas con unas 49 ciudades.

A continuación se muestra una tabla que contiene los avances en la resolución del TSP según el
número de nodos/ciudades a tener en cuenta [19]:

Año Autores Nodos

1954 G. Dantzig, R. Fulkerson, S. Johnson 49
1971 M. Held, R. Karp 64
1977 M. Grötschel, M. Padberg 120
1980 H. Crowder, M. Padberg 318
1987 M. Padberg, G. Rinaldi 532
1987 M. Grötschel, O. Holland 666
1987 M. Padberg, G. Rinaldi 2.392
1992 V. Chvátal, W. Cook (Concorde) 3.038
1994 D. Applegate, R. Bixby, V. Chvátal, W.Cook 7.397
1998 D. Applegate, R. Bixby, V. Chvátal, W. Cook (Concorde) 13.509
2001 D. Applegate, R. Bixby, V. Chvátal, W. Cook, K. Helsgaun (Concorde) 15.112
2004 D. Applegate, R. Bixby, V. Chvátal, W. Cook, K. Helsgaun (Concorde) 24.978
2006 V. Chvátal, W. Cook (Concorde) 85.900
2009 Y. Nagata 1.000.000
2010 S. Lin, B. Kernighan, K. Helsgaun 1.904.711

Tabla 2.1: Evolución del TSP según el número de ciudades a considerar.

* Concorde es un código informático en lenguaje de programación C para resolver el problema

del viajante simétrico u otros problemas de redes relacionados.

2.3. Formulación TSP

En la sección anterior, se expuso la primera formulación del TSP sobre un grafo (Figura 2.2)
tratando de encontrar un ciclo hamiltoniano, pero además, el problema puede ser formulado como
un problema de Programación Lineal Entera, como se verá a continuación.

El TSP puede ser formulado como un grafo completo no dirigido(sTSP) o dirigido (aTSP) en
donde las aristas A = {(i, j) ∀i, j ∈ V, i ̸= j} tienen distintos pesos o costes. Teniendo en cuenta
la Definición 1.1.8 puede construirse la matriz de costes C = [cij ] de dimensión n × n, por ser n

6
las ciudades distintas que existen, donde cada cij representa el peso de la arista (i, j). En forma
matricial se puede ver como:
 
∞ c12 . . . c1n
 c21 ∞ . . . c2n 
C= .
 
.. .. .. 
 .. . . . 
cn1 cn2 . . . ∞

En el caso del sTPS se cumplirá que cij = cji por tratarse de un grafo no dirigido. Los elementos
cij siempre tomarán valores positivos y en caso de que dos vértices no tengan una arista que los
conecte, cij tomará un valor arbitrariamente grande para que el algoritmo no tenga en cuenta dicha
arista. Cuando el grafo es completo, los únicos elementos que tomarán “valor ∞ ” serán los de la
diagonal, ya que no es posible ir de una ciudad en sı́ misma.
Teniendo un grafo, la construcción de la matriz de costes se realiza de manera muy sencilla.

Ejemplo 2.3.1. Supongamos que se tiene el grafo de 6 vértices que se muestra a continuación:

Figura 2.3: Grafo ponderado de 6 vértices

Las aristas punteadas no existı́an en el grafo inicial, por lo que se incluyen de manera artificial con
“coste infinito”para la creación de la matriz de costes asociada:
 
∞ 18 5 ∞ 9 4
18 ∞ 7 5 ∞ 6 
 
 5 7 ∞ 10 2 1 
C= ∞ 5 10 ∞ 7 3 

 
9 ∞ 2 7 ∞ 8
4 6 1 3 8 ∞

Por otra parte, formulándolo desde la Programación Lineal Entera, se plantea el problema de la
siguiente manera:

Se enumeran las ciudades de 1, . . . , n. Para cada par de ciudades (i, j) se consideran xij ,
variables booleanas, de modo que tomarán el valor 1 si la arista (i, j) está en la solución y 0
en otro caso.
Sea cij la distancia desde la ciudad i a la ciudad j.
El objetivo del problema y, por tanto, función a optimizar es conseguir, como se vino diciendo
hasta ahora, un ciclo Hamiltoniano de peso mı́nimo, por tanto, restricciones a tener en cuenta
serán que cada vértice tenga tan solo una arista de salida y una de entrada.

7
 Con las variables, restricciones y objetivo descritos se construye el modelo general de Programa-
ción Lineal en Enteros siguiente:

X
(P0 ) Minimizar f = cij xij (2.1)
i,j∈V

sujeto a
X
xij = 1 ∀i ∈ V \ {j} (2.2)
j∈V
X
xij = 1 ∀j ∈ V \ {i} (2.3)
i∈V
xij = {0, 1} ∀i, j ∈ V (2.4)

La ecuación (1) será la función objetivo, mientras que la (2) y la (3) serán las restricciones que
impondrán que se salga y se llegue una sola vez a cada ciudad respectivamente.

Notación 1. Dada una solución factible T del problema P0 se puede representar por el conjunto de
aristas que forman
P el ciclo, T = [(i1 , i2 ), (i2 , i3 ), . . . , (in−1 , in )(in , i1 )] .El coste total de dicho ciclo
será z(T ) := (i,j)∈T cij . Además, si T es una solución óptima, z(T ) tendrá el menor coste posible
y será, por tanto, una cota superior de v(P ) (valor óptimo de P ).
Observación 2.3.2. Por como se define z(T ), se escoge uno y solo un elemento (coste) en cada
fila y cada columna de la matriz C para contrubuı́r al coste final de la solución.
Por otro lado, para el caso simétrico, la solución T = [(in , in−1 ), (in−2 , in−3 ), . . . , (i2 , i1 )(i1 , in )]
tendrá coste, z(T ) = z(T ) aunque se traten de dos soluciones distintas.
Para ambos casos (simétrico o asimétrico) si xij = 1 automáticamente xji = 0 y viceversa

Aunque estas restricciones son necesarias, no son suficientes ya que se podrán formar ciclos
disjuntos como ocurre en esta posible solución de (P0 ) sobre la Figura 2.3:

Figura 2.4: Ciclos disjuntos con 6 nodos

En la Figura anterior, x12 = x23 = x31 = x45 = x56 = x64 = 1 por lo que se cumplen las restricciones
de (P0 ), sin embargo, no es un circuito hamiltoniano y, por tanto, no es una solución factible para la
resolución del problema que se quiere plantear. La ruptura de estos subcircuitos puede ser abordada
de varias formas introduciendo distintas restricciones.

Para cada subconjunto S ⊆ V con 2 ≤ |S| ≤ |V | − 1, una condición suficiente para evitar
subciclos entre los vértices de S es la propuesta por Dantzig, Fulkerson y Johnson (DFJ):

8
 X
xij ≤ |S| − 1 (S ⊂ V, 2 ≤ |S| ≤ |V | − 1)
i,j∈S

A continuación se demuestra que esta restricción serı́a suficiente para la eliminación de los ciclos
de menor longitud.
Proposición.
P 2.3.3. Dado un subconjunto S ⊆ V y una solución factible T del problema P0 . Si
i,j∈S x ij ≥ |S| =⇒ la solución T contiene al menos un ciclo entre los vértices de S .

Observación 2.3.4. Si |V (G)| = n y |A(G)| > n − 1 =⇒ G tiene al menos un ciclo.

Demostración de Proposición 2.3.3. Sea G el grafo, G[S] elPsubgrafo inducido por los vértices de S
(Definición 1.1.11) y sea TS = T ∩ A(G[S]). La condición i,j∈S xij ≥ |S| =⇒ TS tiene al menos
|S| aristas.
Por la Observación 2.3.4 como TS tiene más de |S| − 1 aristas =⇒ TS tiene al menos un ciclo.

Añadiendo esta restricción, el problema queda planteado como :

Problema 1 (Formulación del TSP).
X
(P) Minimizar f = cij xij (2.1)
i,j∈V

sujeto a
X
xij = 1 ∀i ∈ V \ {j} (2.2)
j∈V
X
xij = 1 ∀j ∈ V \ {i} (2.3)
i∈V
X
xij ≤ |S| − 1 (S ⊂ V, 2 ≤ |S| ≤ |V | − 1) (2.4)
i,j∈S

xij = {0, 1} ∀i, j ∈ V (2.5)

Observación 2.3.5. Los métodos de resolución que se exponen en las secciones sucesivas, serán
planteados basándose en esta formulación.
Definición 2.3.6. A lo largo del trabajo, a veces se llamará tour a una solución factible T y se
llamará cadena a los camino disjuntos más largos que forman parte de la solución parcial.

La mayorı́a de los documentos sobre este problema, hablan sobre su carácter NP-completo que
hace que esté en cotinuo estudio hasta la obtención de un método lo suficientemente eficiente. Para
entender este concepto es necesario introducir la teorı́a de complejidad.

2.4. Complejidad Temporal

La Teorı́a de la Complejidad, estudia la forma de clasificar problemas algorı́tmicos según la
dificultad para resolverlos. La clasificación tiene principalmente dos factores en cuenta: el tiempo que
tarda el algoritmo en ejecutarse y la memoria necesaria para almacenar los datos tras la ejecución
del problema.[23, 2, 24]

Definición 2.4.1. Se denomina función de complejidad a T : N −→ N donde T (n) es el tiempo

transcurrido en la ejecución del algoritmo. Un problema resuelto en tiempo polinómico equivale a
decir que dado un n ∈ N, T (n) está acotado por una función polinómica.

9
 Para poder clasificar un problema según su complejidad algorı́tmica, debe estar planteado como
un problema de decisión, es decir donde las únicas respuestas posibles sean “Sı́” o “No”.
El TSP planteado como un problema de decisión:

“¿Dado un número p ∈ N, existe un circuito Hamiltoniano de longitud menor o igual a p?”

La complejidad temporal puede ser usada para describir la cantidad de pasos necesarios para
resolver un problema o para describir cuanto tiempo lleva verificar si una solución es o no la óptima.
Definición 2.4.2. La clase P está formada por el conjunto de los problemas que pueden ser
resueltos con algoritmos en tiempo polinómico explorando una a una las distintas posibilidades.
Definición 2.4.3. La clase NP abarca el conjunto de todos los problemas de desición que pueden
ser resueltos con algoritmos de tiempo polinómico explorando todas las posibilidades a la vez hasta
encontrar la correcta dicho de otro modo, si la respuesta al problema es “Sı́” puede ser confirmado
en tiempo polinomial.
A dı́a de hoy, aún no se sabe si P = N P . Es uno de los problemas abiertos que el Instituto Clay
expone en los problemas del milenio en [3]. A raı́z de las distintas posibles soluciones que recibió
este problema se define:
Definición 2.4.4. Un problema es NP-completo si es NP y todo problema NP se puede reducir1
en tiempo polinomial en NP-completo.
Por lo que el TSP según su complejidad algorı́tmica temporal, Richard M. Karp en 1972 lo
definió como NP-completo.

Al categorizarse en este grupo, es complicado obtener la solución exacta en un tiempo compu-

tacionalmente bajo. Se utilizan métodos que pretenden obtener una solución óptima en el menor
tiempo posible.
Un breve resumen de la clasificación general de los métodos es la que se muestra a continuación:

2.5. Técnicas de Resolución

Algoritmos Exactos: Métodos que buscan obtener la solución óptima. Suelen aplicarse a
problemas sencillos ya que el tiempo de ejecución es muy elevado.
Algoritmos Heurı́sticos: Métodos que tratan de encontrar una buena aunque no siempre
óptima solución al problema con bajos tiempos de ejecución. Son empleados en problemas
donde la rapidez y la calidad de la solución priman de igual manera. Suelen ser algoritmos
ad hoc o lo que es lo mismo, se diseñan para un tipo de problema en especı́fico por lo que
consiguen acercarse a la solución en un tiempo razonable. Son utilizados cuando no se conoce
ningún algoritmo exacto que pueda resolver el problema o porque usandolo, es demasiado
costoso computacionalmente. A lo largo del documento, se empleará la palabra heurı́stica
para referirse a estos métodos.
Una breve clasificación dentro de este grupo es la siguiente:

• Constructivos: Construyen paso a paso una solución factible añadiendo en cada iteración
un nuevo elemento. Algunos de los más conocidos son:
◦ Heurı́sticos Voraces: Eligen en cada paso el mejor mı́nimo local. Toman decisiones
que parecen ser óptimas en el momento pero sin considerar las posibles consecuencias
a largo plazo.
1 Un problema de decisión R es polinomialmente reducible a un problema de decisión Q si hay una transformación

en tiempo polinomial de cada ejemplo IR del problema R a un ejemplo IQ del problema Q, donde los ejemplos de IR
y IQ tienen la misma respuesta

10
 ◦ Estrategia de descomposición: Divide el problema inicial en subproblemas más
sencillos de resolver individualmente. La solución final del problema original es la
combinación de las soluciones de los subproblemas.
• Métodos de búsqueda : Parten de una solución factible y la intentan mejorar proguesiva-
mente obteniendo en cada paso una nueva solución.
◦ Heurı́sticos de mejora iterativa: Utilizando una solución inicial se extraen k aristas
con el objetivo de añadirlas de modo que mejore la solución de partida.
◦ Heurı́sticos de multiarranque: Con ayuda de probabilidades y procesos aleatorios
mejoran la solución de partida.

El algortimo voraz más conocido e importante es el del Vecino más próximo que comienza en
una ciudad arbitraria y mientras haya ciudades no visitadas, se visita la ciudad más cercana que
todavı́a no haya aparecido en el recorrido (volviendo por último a la primera ciudad).

2.5.1. Vecino más próximo

Este método construtivo voraz trata de crear un ciclo Hamiltoniano de bajo coste. Dado un
vértice, se busca la arista con menor peso que sale de dicho vértice. Las primeras aristas escogidas
que aparentemente tienen poco peso pueden causar malas elecciones al final del proceso ya que sólo
quedarán aristas con peso muy elevado.

Algorithm 1 Algoritmo Vecino más cercano

INPUT : Matriz costes C

1. k=0 Seleccionar un vértice j ∈ V al azar para inicializar v0 = j y sea W = V .

2. k=k+1 Buscar el vértice más cercano a vk en W = V \ {vértices ya visitados}. Si el ciclo está
formado por los vértices v1 , v2 , . . . , vk , se busca el vértice que pasará a formar parte del ciclo
vk+1 tal que cvk vk+1 = min{cvk i con i ∈ W }.
3. Se sigue el proceso hasta que W = ∅. Una vez finalizado, vn será unido a v1 (considerando v1
el primer vértice y vn el último).

A continuación, para ver de manera gráfica este método, se representa un ejemplo con un grafo
de 5 vértices (partiendo del vértice A) en el que aparece resaltada la solución final:

Figura 2.5: Ciclo Hamiltoniano en grafo de 5 vértices creado mediante el método Vecino más cercano

11
 El coste total del recorrido por este método será 1 + 2 + 3 + 5 + 9 = 20 . Una mejor elección
del vértice E habrı́a sido la arista que sale hacia A con peso 2, sin embargo ese vértice ya formaba
parte del recorrido por lo que se tuvo que añadir una arista de mayor coste.
La solución óptima al problema, serı́a el ciclo [ABCED] cuyo coste es 1 + 2 + 7 + 5 + 2 = 17.

El auge en la resolución de este problema mediante diversas técnicas se debe a su importancia

práctica en la planificación de rutas, la secuenciación del genoma o la producción de chips entre otros.
Encontrar soluciones óptimas y eficientes puede tener un gran impacto en la eficiencia operativa. Por
lo tanto, la gran variedad de métodos desarrollados refleja la necesidad de abordar este problema
desde múltiples enfoques para satisfacer las necesidades especı́ficas de las distintas aplicaciones.

12
Capı́tulo 3

Ramificación y acotación (Branch

and Bound)

Los métodos exactos garantizan encontrar la solución en un número acotado de pasos. La manera
más intuitiva de resolver un problema de tales caracterı́sticas es calcular las n! posibles rutas con
sus costes asociados y concluir cuál tiene menor coste, es decir, mediante la búsqueda por fuerza
bruta. Sin embargo, para n = 9 ya se calcuları́an 9! = 362880 rutas, por lo que tendrı́a un coste
computacional demasiado elevado. Es por esto que se proponen técnicas más rápidas y menos
costosas como es el algoritmo de ramificación y acotación [25].

El término ramificación y acotación fue introducido en 1963 por John D. C. Little, Katta G.
Murty, Sweeney y Karel [16]. Como indica su nombre, el algoritmo consta de la realización de dos
pasos.

Primeramente, el método consiste en resolver el problema de optimización diviendo la región

factible en subconjuntos más pequeños determinados por soluciones parciales y representados
como nodos de un árbol de decisión. El proceso de división de los subconjuntos (de cada uno
de los nodos) radica en la ampliación de la solución parcial fijando una variable no fijada.
Esto es lo que se llama ramificación del nodo, pudiendo subdividirse cada uno en al menos
dos subconjuntos no vacı́os, puesto que cada variable toma dos posibles valores 0 o 1.

En este proceso de ramificación se van calculando cotas inferiores y superiores del valor de v(P )
en cada uno de los subproblemas, asociadas al coste de la solución parcial que corresponda.
Esto permite en última instancia decidir si una solución factible es óptima.

El algoritmo puede ser empleado para la resolución de problemas pequeños que contienen entre
40 y 60 ciudades. [12]

Este método puede plantearse de maneras muy diversas dependiendo de la relajación aplicada
sobre el problema (P ) o la regla de ramificación que determinará el orden de selección de los
subproblemas y procedimiento de cálculo de cotas inferiores.

A continuación se determinará el planteamiento y metodologı́a elegida.

13
3.1. Planteamiento

3.1.1. Relajación del problema

Considérese el problema (P ) con C = [cij ] la matriz de costes. La relajación del problema,
consiste en la eliminación de la restricción (4):
X
(AP) Minimizar cij xij (1)
i,j∈V

sujeto a
X
xij = 1 ∀i ∈ V \ {j} (2)
j∈V
X
xij = 1 ∀j ∈ V \ {i} (3)
i∈V
xij = {0, 1} ∀i, j ∈ V (5)

de modo con T factible las variables tomarı́an los valores

(
1 (i, j) ∈ E
xij =
0 (i, j) ∈ I

con E las aristas incluidas en la solución e I las que no lo están.

La elección de la eliminación de la restricción (4) se basa en la volumetrı́a de restricciones a añadir.
A priori, podrán formarse ciclos de menor longitud en la solución. Sin embargo, se evitarán de
manera manual como se verá más adelante.

Encontrar una solución óptima al problema relajado (AP ) proporcionará una cota superior para
v(P ) ,es decir,v(AP ) ≥ v(P ) ,por ser este menos restrictivo que el problema original.

Además, un modo de obtener una cota inferior de v(P ) en (AP ) será mediante la reducción de
la matriz de costes de la siguiente forma:

Definición 3.1.1. Una fila (columna) de una matriz se dice reducida si sus elementos son no
negativos y contiene al menos un cero.
Definición 3.1.2. Una matriz se dice reducida si cada fila y cada columna es reducida.
Sea C = [cij ] la matriz de costes asociada a un grafo ponderado. Se reducen las filas de modo
que para cada fila i, ui = min{cij : j ∈ V } luego ĉij := cij − ui ∀i, j ∈ V luego, Ĉ := [ĉij ].

De igual manera para la reducción de columnas. Para cada columna j, vj = min{ĉij : i ∈ V }

por lo que c̄ij := ĉij − vj . Ası́, se obtiene la matriz reducida por filas y columnas Cred := [c̄ij ].
Veamos un ejemplo sobre una matriz C:

     
∞ 8 25 3 ∞ 5 22 0 ∞ 5 14 0
8 ∞ 12 4  u=[3,4,12,3] 4 ∞ 8 0 v=[0,0,8,0] 4 ∞ 0 0
C=  → Ĉ =   → Cred = C̄ =  
25 12 ∞ 15 13 0 ∞ 3 13 0 ∞ 3
3 4 15 ∞ 0 1 12 ∞ 0 1 4 ∞

P P
El valor LB = i∈V ui + j∈V vj será una cota inferior de v(P ):

14
 P
Proposición.
P P 3.1.3. LB es una cota inferior del problema. O lo que es lo mismo, i∈V ui +
j∈V v j ≤ i,j∈V cij xij con la definición de ui , vj anterior.

Demostración. La obtención de una solución factible supone la elección de exactamente un elemento

por fila y por columna de la matriz de costes 2.3.2. Sea ui el mı́nimo de cada fila de C, vj el mı́nimo
de cada columna de Ĉ y z(T ) el coste total con T =[(i1 , i2 ), . . . , (ik−1 , ik )] una solución factible.
Si k = nPentonces ik i1 : P
n n Pn Pn
z(T
Pn ) = k=1 cik ik+1 = k=1 (cik iP k+1
− uk + uP
k) = k=1 ĉik ik+1 + Punk =
k=1
n n Pn
Pk=1 (ĉik ik+1 − v k+1 + v k+1 ) + u
k=1 kP = (c̄
k=1 P ik ik+1 + v k+1 ) + k=1 uk = k=1 c̄ik ik+1 +
n P n P n n n
u
k=1 k + v
k=1 k+1 = c̄
k=1 ik ik+1 + u
k=1 k + k=1 k v
Pn Pn
Como c̄ik ik+1 ≥ 0 =⇒ z(t) ≥ k=1 uk + k=1 vk ∀ T solución factible.

3.1.2. Regla de Ramificación Little

En este apartado se explicará la técnica de ramificación escogida ası́ como el orden del siguiente
subproblema a ramificar y la construción de la cota inferior de cada uno de ellos.
A lo largo de la historia del TSP se desarrollaron diversas técnicas de ramificación considerando
una buena regla la que genera pocos hijos de cualquiera de los nodos del árbol de decisión y por
tanto que deja excluı́das muchas soluciones no óptimas. La primera regla expuesta por John D. C.
Little, Katta G. Murty, Sweeney y Karel genera únicamente dos hijos por cada nodo del árbol de
decisión y sigue siendo la más efectiva hasta el momento.

Esta regla será por tanto la que se exponga a continuación:

Definición 3.1.4. Sea el subproblema relajado APk (iteración k) y su matriz de costes reducida
dada por Cred := [cij ], se define penalización sobre cada cij = 0 como:

Pij = min{cih : h ∈ V \ {j}} + min{chj : h ∈ V \ {i}} (3.1)

Una arista (r, s) ∈ A \ (Ek ∪ Ik ) será elegida si Prs = max{Pij : cij = 0} es decir, encontrar la
arista (r, s) consistirá en buscar una arista que no fue escogida todavı́a y que tiene el cero con la
penalización más alta.
Esta elección se realiza de modo que el algoritmo se centra en estudiar primero los subproblemas de
menor coste.

Definición 3.1.5. La penalización fila se define como Pi = min{cih : h ∈ V \ {j}} .

Definición 3.1.6. La penalización columna como Pj = min{chj : h ∈ V \ {i}}.
En caso de haber varios ceros con la misma Pij , se escogerá sin pérdida de generalidad el que
tenga mayor penalización fila.
Una vez elegida la arista, se ramificará generando dos hijos del nodo k; los nodos k + 1 y (k + 1)′
definidos como:

Ek+1 = Ek ∪ {(r, s)} Ik+1 = Ik

E(k+1)′ = Ek I(k+1)′ = Ik ∪ {(r, s)}

15
 Figura 3.1: Ramificación de un nodo

Cada uno de estos nuevos nodos, tendrán asociados una nueva matriz de costes reducida y una
cota inferior de la solución parcial. Dependiendo de si se esta en la rama de la derecha o en la rama
de la izquierda, estos elementos se calcularán de manera diferente:

Si xij = 1 =⇒ (i, j) ∈ Ek
• Matriz de costes

1. Eliminación de fila y columna Se sustrae la fila i y la columna j para asegurarse de

que no se elijan más aristas que salgan de i o lleguen a j.
2. Elemento cji = ∞ Para evitar la elección de la misma arista en el otro sentido
3. Evitar ciclos de menor longitud Si después de incluir la arista (i, j), incluir otra
arista (k, l) formarı́a un ciclo de menor longitud en la solución, se impone ckl = ∞.
Por ejemplo si se tiene la cadena [3 2 5 7] se impondrá c73 = ∞.
4. Reducción de la matriz Después de las transformaciones anteriores, se vuelve a
reducir la matriz en caso de que sea necesario que tendrá un coste de reducción
asociado costered .
• Cota inferior (LB) Será la cota inferior del nodo anterior más las reducciones de fila y
columna que resulten del último paso de la transformación en la matriz de costes, es
decir , LBk+1 = LBk + costered
′
Si xij = 0 =⇒ (i, j) ∈ Ik+1
• Matriz de costes Se impone cij = ∞ y se reduce la matriz cuyo coste de reducción
coincidirá con Pij .
• Cota inferior (LB) Será la LB del nodo anterior más la penalización del 0 correspondiente
por no haber escogido esa arista:
′
LBk+1 = LBk + Pij
.

Ahora, será necesario un criterio que determine que nodo de los dos nodos hijos generados será
el siguiente en ramificar.

3.1.3. Exploración del árbol de decisión

Se sigue un algoritmo de búsqueda en profundidad (Definición 1.1.17).

Definamos primero la terminologı́a empleada en los problemas de ramificación:

Definición 3.1.7. Un nodo vivo es un nodo del espacio de soluciones del que aún no se han
generado todos sus hijos.

16
Definición 3.1.8. Un nodo muerto es un nodo del que no se van a generar más hijos debido a
que:
No hay más hijos posibles
No producirá una solución mejor que la actual

Definición 3.1.9. Un nodo en expansión es un nodo del que se están generando hijos.
Inicialmente, se hará una búsqueda en profundidad “a ciegas” utilizando el método LIFO(Last
In, First Out). En este método, se van expandiendo los nodos siguiendo la rama en la que xij = 1
∀i, j ∈ V hasta que esta rama muera porque todos los vértices fueron visitados (|In | = ∅ y |En | = n),
obteniendo una primera solución factible T y su LB correspondiente, que será una cota superior de
v(P ). A esta cota superior se le denominará Z .
Para explorar los nodos de las ramas en las que xij = 0, se aplicará de nuevo el proceso anterior
para el nodo vivo que tenga menor LB asociada con la salvedad de que se parará de explorar
cuando el nodo a ramificar muera por tener LB ≥ Z. Por el contrario si en esta exploración se
encuentra una solución factible de menor coste, este será la nueva Z.

El proceso de ramificación terminará cuando todos los nodos vivos hayan sido explorados o
desechados. Obteniendo la solución óptima y su coste asociado.
Observación 3.1.10. Cuando se parte de un nodo en el que ya se han añadido n − 2 aristas, solo
se puede cerrar ciclo con las dos aristas restantes por lo que se escogerán ambas al tiempo.

17
 La subdivisión del conjunto de posibles soluciones de la región factible del problema P , se
representará como se ve en la figura siguiente:

Figura 3.2: Ramificación en árbol del problema

En la figura anterior se expanden en este orden los problemas relajados AP0 −→ AP1 −→ · · · −→
APn−2 y se crean los nodos vivos AP1′ , AP2′ , APn−2
′
.En las siguientes iteraciones el nodo de menor
′
LBk pasa a ser el nuevo nodo en expansión.

18
Algorithm 2 Algoritmo Branch and Bound

INPUT : Matriz de costes C y problema P

NODO RAIZ (PASO 0):

1. Reducir la matriz C por filas y columnas.
2. Calcular una cota inferior LB0 sumando las reducciones de filas y columnas.
3. Inicializar los conjuntos E = ∅ e I = ∅;

RESTO DE NODOS (PASO k):

1. Calcular la penalización Pij para cada cero de la matriz de costes.
2. Escoger la arista con el cero de mayor penalización. Si hay más de uno, escoger la de mayor
penalización fila.
3. Ramificar el árbol de búsqueda en el nodo calculando las cotas inferiores según corresponda:

1. Si la arista (i, j) no está en la solución, xij = 0:

a) Incrementar el coste del nodo anterior en la penalización del cero correspondiente,

LBk′ = LBk−1 + Pij .
b) Asignar un valor arbitrariamente grande (∞) en cij y reducir la matriz resultante.

c) Añadir (i, j) a Ik′ .

2. Si la arista (i, j) está en la solución, xij = 1:

a) Eliminar la fila i y la columna j de la matriz de costes. En la práctica, poner ∞ en

esas posiciones.
b) Asignar un valor arbitrariamente grande (∞) en cji y asegurarse de que no se formen
ciclos de menor longitud asignando valores ∞.
c) Reducir la matriz.
Pn Pn
d ) Calcular la cota inferior de dicha reducción: LBk = LBk−1 + i=1 ui + j=1 vj .
e) Añadir (i, j) a Ek .

4. Continuar la ramificación en la rama en la que xij = 1 repitiendo el proceso 1-4 hasta obtener
una matriz 2 × 2.
5. Reducir la matriz si fuese necesario e incrementar la cota. Ambas aristas con coste cero
formarán parte de la solución.
6. La última cota inferior obtenida será una cota superior Z del problema P .
7. Si LBk′ ≥ Z en todas las ramas, se acaba el proceso y el ciclo obtenido en el paso 6 será la
solución óptima.
8. En otro caso, mientras LB < Z repetir el proceso 1-6 para el nuevo nodo raı́z. Si en el proceso
se encuentra una nueva solución factible con menor coste, se actualizará al igual que el nuevo
coste.
El algoritmo termina cuando todos los nodos mueran por tener LB ≥ Z.
=0

19
3.2. Implementación en Matlab
Este método fue implementado en Matlab del modo que se muestra a continuación. Será expuesto
el orden de llamada a las distintas funciones, una breve descripción de cada una de las funciones asi
como de sus ’inputs‘ y ’outputs’:

Esquema de llamadas:

function [mejorsol, costesol, tiempo branch, tiempo] = BB TSP(matrizcostes)

• function [LB, E, Eord, LB I lista, matriz reducida lista, E lista, Eord lista]
= BB base(matrizred, LB, Z, E, Eord)
◦ function [max epenalizacion,fila,columna]=calcula epenalizacion(matriz)
◦ function [Eord, e1, e2] = verificar exciclo(Eord, arista nueva)
◦ function [matrizred, costered] = reducematriz(matrizcostes)
• function [matrizred, costered] = reducematriz(matrizcostes)
Breve resumen de cada función:
BB TSP: Función principal que implementa el método de branch and bound, tan solo es
necesaria una matriz de costes asociada a un grafo ponderado.
BB base: Función que a partir de un nodo raı́z con su matriz reducida y cota inferior que da
dicha solución hasta el momento, hace la ramificación a derecha hasta encontrar una solución
completa (aunque no necesariamente óptima).
En caso de tener más de dos input, ya existe una solución completa y se está ramificando a
derecha nodos vivos previamente almacenados.
reducematriz: Función que, a partir de una matriz de costes, obtiene la matriz reducida y el
coste de dicha reducción (suma de los valores para la reducción de filas y columnas).
calcula penalizacion: Función que a partir de una matriz, calcula la máxima penalización
además de las penalizaciones fila y columna.

verificar exciclo: Función que a partir de una serie de cadenas de aristas ya añadidas a la
solución, comprueba como poder “pegar” al añadir una arista nueva para evitar que se cree
un ciclo de menor longitud. Tiene como objetivo “colocar” esa nueva arista añadida en la
cadena correspondiente (uniendo dos cadenas si se diese el caso) además de reconocer los
extremos de los nuevos ciclos de longitud menor a n que se deben evitar.

Descripción de las variables empleadas:

matrizcostes : Matriz simétrica de tamaño n × n.

mejorsol : Vector fila con las aristas que están incluidas en la solución
costesol : Valor con el coste de la solución
tiempo branch : Valor en segundos del “branching a derecha”

tiempo : Valor en segundos correspondiente al tiempo de ejecución

matrizred : Matriz simétrica reducida de tamaño n × n
costered : Valor del coste de la reducción
LB : Cota inferior de v(P)

Z : Cota superior de v(P)

E : Matriz n × 2 que contiene las aristas incluidas en la solución

20
 LB I lista: Vector fila que almacena las cotas inferiores correspondientes a las soluciones
parciales de los nodos de la izquierda
matriz reducida lista : Array que tendrá almacenadas las matrices de los nodos de la
izquierda

E lista: Array con las aristas que se van incluyendo para ir creando la solución parcial, por
tanto, en el orden de inclusión
Eord lista: Array con las aristas de E ordenadas en cadenas disjuntas(en caso de no estar la
solución completa)
Eord : Array que contiene en las celdas una matriz k × 2, donde están concatenadas k aristas

e1,e2: Arista a evitar para que no forme ciclo

arista nueva : Vector 2 × 1 que se corresponde con la arista que es añadida a la solución en
cada paso.
max penalizacion: Valor de la penalización máxima y, por tanto, la elegida

fila: Índice de la fila de la matriz que contiene el 0 con mayor penalización

columna : Índice de la columna de la matriz que contiene el 0 con mayor penalización

3.3. Ejemplo ilustrativo para sTSP con Branch and Bound

A continuación se muestra un grafo con 4 nodos y pesos en sus aristas para ilustrar el
procedimiento de este método:

Figura 3.3: Grafo completo de ejemplo con 4 vértices

El grafo tiene asociada la matriz de costes siguiente:

 
∞ 8 3 15
0
8 ∞ 7 1
C =C= 3 7 ∞

2
15 1 2 ∞

Observación 3.3.1. La notación empleada para las matrices será, C k para las matrices asociadas
a los nodos de la derecha y C k′ para los de la izquierda. El superı́ndice k representa el nivel en el
que se encuentra asumiendo que el nodo raı́z está en el nivel 0. En caso de que la matriz sea
reducida tendrá el subı́ndice red.

21
El primer paso será obtener una matriz reducida:
Para ello, se calcula el mı́nimo de cada fila obteniéndose u = (3, 1, 2, 1) y restando en cada fila:
 
∞ 5 0 12
 7 ∞ 6 0
Ĉ 0 = 
1 5 ∞

0
14 0 1 ∞

Se hará lo mismo con las columnas de Ĉ 0 , siendo v = (1, 0, 0, 0) y restando en cada una de las
columnas:
 
∞ 5 0 12
0
6 ∞ 6 0
Cred = 0 5 ∞ 0


13 0 1 ∞
P P
LB0 = ui + vj = 3 + 1 + 2 + 1 + 1 + 0 + 0 + 0 = 8.
Ahora se calculará la penalización de cada uno de los ceros. Esta será la suma del valor mı́nimo de
la fila y la columna en la que se encuentre.

0(5+1) 0(6) 12
   
∞ 5 12 ∞ 5
0
 6 ∞ 6 0(6+0) 
= 6
 ∞ 6 0(6) 
Cred = (0+6) (0+0) (6)

0 5 ∞ 0  0 5 ∞ 0(0) 
13 0(1+5) 1 ∞ 13 0(6) 1 ∞

En la primera matriz se muestran los ceros con la penalización desglosada en penalización fila
y penalización columna, ya que como varios 0 tienen la misma penalización (en la matriz de la
derecha), se escogerá el 0 con mayor penalización fila.
En este caso, teniendo en cuenta todo lo anterior, la primera arista que ramificará el árbol será la
(2,4) y se seguirá el orden que se muestra:

Figura 3.4: Árbol correspondiente a la Figura 3.3

En la Figura 3.4 se representan las cotas inferiores en cada paso como LBk o LBk′ (el k representa
la profundidad en la que se encuentra en el árbol) según corresponda a ramas de la derecha o de la
izquierda respectivamente que son calculadas de distinta manera como se puede mostrar en el
algoritmo que se muestra en el apartado 3.1.2.

22
Ilustremos de forma práctica el procedimiento:

1 2 3 4
1 ∞ 5 0 12 1 2 3 
1 1 ∞ 5 0
C =  6 ∞ 6 0 =
2  
3 0 5 ∞
3  0 5 ∞ 0
4 13 ∞ 1
4 13 0 1 ∞

Como la arista (2, 4) ya está en el recorrido, la arista (4, 2) no podrá hacerlo y por ello se tomará
artificialmente un coste arbitrariamente grande (c42 = ∞).
Se comprueba que la matriz resultante sea reducida. La segunda columna y la tercera fila no tienen
ningún 0 por lo que se extraerá 5 y 1 respectivamente. Quedando ası́ la cota inferior
LB1 = 8 + 5 + 1 = 14 y
1 2 3 
1 1 ∞ 0 0
Cred =
3  0 0 ∞
4 12 ∞ 0
La cota inferior LB1′ , casualmente toma el mismo valor que LB1 pero no es obtenida de la misma
forma. Esta cota, es la suma de LB0 y 6 ( penalización de no escoger el 0 de la posición (2, 4) de la
matriz) y tiene la matriz asociada siguiente:

1 2 3 4
1 ∞ 5 0 12
C 1′ = 2 
6 ∞ 6 ∞ 
3 0 5 ∞ 0
4 13 0 1 ∞

Se sigue ramificando por el nodo de la derecha hasta obtener una primera solución sea o no la
1
óptima. Se calculan las penalizaciones de los ceros de Cred , se escogerá la siguiente arista donde el
0 tenga la mayor penalización y se repetirá el proceso anterior.

0(0+0) 0(0+0) 1 2 3 
 
∞ 1 2 
1 1 ∞ 0 0
Cred = 0(0+12) 0(0+0) ∞  =⇒ C 2 = = 1 ∞ 0
(12+0) 3 0 0 ∞
12 ∞ 0 3 0 ∞
4 12 ∞ 0

En este caso, el coste de la arista (3, 4) no aparece, ya que se sustrajo la columna 4 por lo que para
que en la solución parcial 2 − 4 − 3 se cumpla la condición (4) de P , será necesario que (3, 2) no
forme parte del ciclo y, por tanto, que tenga un coste arbitrariamente grande.
La matriz resultante es de tamaño, 2 × 2 luego ambas aristas con coste distinto de infinito
formarán parte de la solución, las aristas (3, 1) y (1, 2) . La cota superior será Z = 14.
Se cumple que LBk′ ≥ Z con k ∈ 1, 2 luego se habrá obtenido la solución óptima [1 2 4 3 1] cuyo
coste será 14.

Figura 3.5: Solución óptima del ejemplo

23
Capı́tulo 4

Lin-Kernighan

Teniendo en cuenta la clasificación explicada en la sección 2.5 sobre los métodos heurı́sticos, el
más destacable de búsqueda local es el algoritmo k − opt que será explicado en más detalle en lo
que sigue. Al tratarse de un algoritmo de mejora local, no puede demostrar la optimalidad de la
solución, pero si proporcionar una confianza estadı́stica adecuada para las aplicaciones prácticas
ya que, cuenta con tiempos de ejecución relativamente bajos. Esta técnica está planteada para
aplicarse sobre instancias del TSP simétrico.[9, 11, 14]

Una técnica heurı́stica de búsqueda local en un problema de Optimización Combinatoria sigue

las siguientes instrucciones:
1. Generar una solución factible al azar, es decir, un ciclo hamiltoniano T que satisface las
restricciones del problema P (Problema 2.3).
2. Haciendo transformaciones sobre T trata de buscar una solución factible T ′ que mejore el
coste de T .
3. Si se encuentra sustitución de T cumpliendo f (T ′ ) < f (T ) donde f es la función objetivo a
minimizar luego, T = T ′ y se repite el paso anterior.
4. Si no se encuentra ninguna mejora posible, T es un óptimo local y se repite el proceso desde
el principio.

La primera heurı́stica utilizada para la realización del paso 2, fue el intercambio 2 − opt propuesto
por Croes en 1958 que consiste en el intercambio de dos aristas por otra dos que mejoren la solu-
ción factible inicial. Años más tarde, se definió una generalización de la misma: el intercambio k −opt.

Definición 4.0.1. Un k-intercambio sobre una solución factible T es un intercambio de k aristas

de T por otras k que no pertenecen a T y tal que, se obtenga una nueva solución T ′ factible.
Se dice que T Y T ′ son ciclos vecinos.
Se dice k − óptimo(o k − opt) si T ′ tiene menor coste que T .
Cuanto mayor sea el valor de k, mayor probabilidad habrá de que el ciclo final T ′ sea óptimo. Sin
embargo, el esfuerzo computacional aumenta rápidamente al aumentar, k por lo que es difı́cil saber
de antemano qué k es el mejor para satisfacer que la solución tenga una buena calidad y a la vez que
tenga un buen tiempo de ejecución. En la práctica, se suele trabajar con 2 − opt y 3 − opt ya que
está demostrado que para valores de k superiores el tiempo extra de computación no compensa en
relación con la mejora de los resultados. Este último enunciado se cita en alguna de las referencias
empleadas, sin embargo, aunque parece lógico y en la implementación se evaluará dicha ´´obviedad”
realizando 2−opt iterativos, no hay documentos al alcance de todos con pruebas que ası́ lo determinen.

El algoritmo de Lin-Kernighan crea una nueva versión del algoritmo k − opt eliminando el
inconveniente de tener que imponer el k de antemano y, a pesar de ser planteado en el año 1973, se si-
gue considerando la mejor técnica heurı́stica disponible para la resolución del TSP simétrico.[9, 11, 14]

24
 Además, el algoritmo original tuvo diversas modificaciones para tratar de mejorar los resultados
y eficiencia del método. Alguna de ellas se pueden encontrar en los documentos que se referencian
más abajo ,pero este trabajo está únicamente centrado en el algoritmo base original (LK) ası́ como
algunos refinamientos que se añadieron posteriormente sobre este.

4.1. Planteamiento
Esta técnica va decidiendo cuál será el valor de k en cada iteración, intentando escoger el mayor
k posible para dar como resultado el ciclo de menor peso. Es decir, se inicializa con una solución
factible y se comprueba la existencia de ciclos vecinos de menor costo considerando un conjunto
creciente de intercambios potenciales, comenzando con k = 2 o k = 3.
Sea A el conjunto de todas las aristas y T el recorrido inicial. El algoritmo trata de encontrar
dos conjuntos de aristas, X = {x1 , x2 , . . . , xk } e Y = {y1 , y2 , . . . , yk∗ } tales que, si las aristas de X
se eliminan de T y se sustituyen por las aristas de Y pertenecientes a A − T , el recorrido T ′ tiene
menor peso que T . Los conjuntos X e Y comienzan vacı́os y se construyen elemento a elemento,
de modo que en el paso i se añaden xi e yi a los conjuntos respectivamente. Las aristas de X se
denominan aristas de salida y las de Y aristas de entrada.

Observación 4.1.1. No cualquier par de conjunto X,Y genera un k-intercambio. Esto se debe a
las diversas formas de reconexión de las partes disjuntas que crea la eliminación de las aristas de X.
Tomemos como ejemplo x1 = (i, j) y x2 = (k, l) aristas no adyacentes que forman parte de T . Si se
eliminan del recorrido inicial, quedarán dos partes disjuntas (A y B) que deberán ser reconectadas.
A continuación se muestran las posibles reconexiones:

Figura 4.1: Movimiento 2-opt

En la Figura anterior se muestran cuatro maneras de reconectar (elegir y1 , y2 ) las aristas

eliminadas, pero, solo las opciones (3)a y (3)b serán válidas (cambian el sentido de recorrido) para
un intercambio 2 − opt ya que la (3)c y (3)d violará la restricción (4) de no formación de ciclos
disjuntos.

A medida que aumentan las aristas a intercambiar, aumentan las posibilidades de reconexión
y, por tanto, la complejidad del algoritmo. Un ejemplo de posible intercambio 3 − opt se verı́a del
siguiente modo:

25
Figura 4.2: Movimiento 3-opt donde las aristas x1 , x2 , x3 son reemplazadas por las y1 , y2 , y3

Los criterios que se siguen para introducir aristas en los conjuntos X e Y son los siguientes:

Criterio de intercambio secuencial : xi e yi tienen que compartir un punto final, al igual

que yi y xi+1 . Si t1 denota uno de los extremos de x1 , se tiene en general xi = (t2i−1 , t2i ),
yi = (t2i , t2i+1 ) y xi+1 = (t2i+1 , t2i+2 ) ∀i ≥ 1.
La secuencia (x1 , y1 , x2 , y2 . . . xk , yk ) corresponde con un camino que alterna aristas de X e
Y denominado ciclo alterno (Figura 4.3) que será representado por la variable XY en lo
que sigue.

No siempre se cumple que un intercambio de aristas sea secuencial, como es el caso en la

Figura 4.4 donde no se pueden numerar las aristas para que el intercambio sea secuencial.

Figura 4.4: Intercambio no secuencial con

Figura 4.3: Ciclo alterno con 8 vértices k=4

Criterio de Factibilidad : Se requiere que xi = (t2i−1 , t2i ) ∀i ≥ 3 se elija de modo que, si t2i
se une a t1 , el camino resultante sea un ciclo. Es decir, que se pueda detener el proceso gene-
rando yi∗ = (t2i , t1 ) y obteniendo ası́ una solución factible al intercambiar X = {x1 , x2 , . . . , xi }
por Y = {y1 , y2 , . . . , yi∗ }.
Una condición necesaria, pero no suficiente para que el intercambio de lugar a un recorrido,
es que la secuencia creada sea cerrada; es decir, que exista yk∗ = (t2k , t1 ).

Criterio de Ganancia Positiva: Cada intercambio de aristas proporcionará una ganancia

gi = |xi | − |yi | donde |xi | y |yi | representan el coste de las aristas correspondientes.
La elección de yi se realiza de modo que la ganancia parcial Gi = g1 + g2 + · · · + gi sea positiva.
De esta forma, cuando se encuentra un gi negativo, por la aditividad de la ganancia, no será
necesario la detención del proceso (Proposición 4.1.2).
Criterio de disyuntividad : Se requiere que los conjuntos X, Y sean disjuntos. Ası́, xi no
puede ser una arista que se unió en un paso anterior e yi no puede ser una arista que fue

26
 eliminada previamente.

Criterio de parada: El criterio de parada decide que no se añadirán parejas en os conjuntos

X,Y cuando Gi ≤ G∗ .

Si el algoritmo termina con G∗ > 0 se habrá encontrado un k − opt que dará una solución factible
T con coste f (T ′ ) = f (T ) − G∗ . Luego, se podrá tomar T = T ′ para poder seguir mejorando la
′

solución.

En caso de terminar el algoritmo con G∗ = 0 se aplicará una técnica de backtracking para

garantizar el cumplimiento del criterio de ganancia positiva y/o el criterio de Factibilidad en caso
de no encontrar un xi o yi que ofrezcan una ganancia positiva y/o un yi∗ que cierre el ciclo, se
empleará una técnica de backtracking para estudiar si una elección (entre las válidas) distinta del
último vértice escogido, solventa el problema. Es importante destacar que esta técnica, aunque
puede resultar muy útil, incrementa el tiempo de ejecución y reduce la eficiencia, por lo que su uso
estará limitado a solo si no encuentra un k − opt a realizar.Esta técnica será detallada más adelante.

Proposición. 4.1.2. Si una secuencia de números tiene una suma positiva, existe una permutación
cı́clica de estos números tal que cada suma parcial es positiva.
Pn
Demostración. Sea la secuencia original g = (g1 , g2 , . . . , gn ) y sea G = i=1 gi > 0. Considerándose
Pi
las sumas parciales como Gi = j=1 gj ∀i = 1, 2, . . . , n se define s como el ı́ndice más grande
tal que Gs−1 es mı́nimo. Sea la permutación cı́clica de la secuencia comenzando desde el ı́ndice
s: (gs , gs+1 , . . . , gn , g1 , g2 , . . . , gs−1 ) = (h1 , h2 , . . . , hn ) = h para reducir la notación. Se pretende
Pi
demostrar que Hi = j=1 hj > 0 ∀i = 1, . . . , n
(
gi+s−1 − si si i = 1, . . . , n − s + 1
Representado, por tanto, hi =
gi−(n−s+1) − si si i = n − s + 2, . . . , n
Siendo Hi cada una de las sumas parciales de la nueva secuencia de modo que:

H1 = h1 = gs

H2 = h1 + h2 = gs + gs+1
.
hola que tal ..
Hn−s+1 = h1 + · · · + hn−s+1 = gs + gs+1 + · · · + gn

Hn−s+2 = h1 + · · · + hn−s+2 = gs + gs+1 + · · · + gn + g1

.
hola que tal ..
Hn = h1 + · · · + hn = gk + gs+1 + · · · + gn + g1 + g2 + · · · + gs−1 = G

Se pretende demostrar que cada Hi > 0 ∀i = 1, 2, . . . n

1≤i≤ Pn − s + 1P
i s+i−1
Hi = j=1 hj = j=s gj = Gs+i−1 − Gs−1 > 0.
Gs−1 es el mı́nimo y i + s − 1 ≥ s luego se cumplirá la desigualdad por la elección de s.
n−sP + 2 ≤ i ≤ nP
i n Pi−(n−s+1)
Hi = j=1 hj = j=s gj + j=1 gj = (Gn − Gs−1 ) + Gi−(n−s+1)
hola = Gn + (Gi−(n−s+1) − Gs−1 ) > 0.
Gn es la suma total que por definición es positiva. Además, se cumple que (Gi−(n−s+1) −
Gs−1 ) ≥ 0 ∀i ∈ {n − s + 2, . . . , n} por el mismo argumento que en el otro caso.

27
4.2. Comentarios sobre el algoritmo
Además de todos los criterios mencionados para mejorar la eficiencia, se añaden los siguientes
refinamientos sobre el backtracking:

El backtracking solo se realiza en el nivel i = 2 e inferior.

Cada t2i+1 (el que determina la arista yi ) está escogido con el menor peso posible.

Para reducir la búsqueda, para los vértices t3 y t5 en caso de ser necesaria la aplicación de
backtracking, este se realizará sobre los 5 vecinos más cercanos, es decir, se estudiarán los 5
y1 e y2 de menor peso respectivamente.
Se permite incumplir el criterio de intercambio secuencial en el nivel i = 2. Se permite que no
exista una arista y2∗ = (t4 , t1 ), con la condición de que exista la arista y3∗ = (t6 , t1 ). De este
modo, el proceso empezarı́a con un 3-intercambio en vez de con un 2-intercambio.
Este método, a pesar de ser como dije el mejor hasta el momento, tiene distintas interpretaciones
según los distintos artı́culos que empleé para este trabajo, luego fue complicado escoger qué formaba
parte del algoritmo original y qué no. El algoritmo básico inicial [15] deja demasiadas decisiones
sin determinar, por lo que se hace bastante difı́cil su interpretación total. Por esta controversia, he
interpretado y tomado las siguientes decisiones sobre los frentes que dejan abiertos:
La referencia citada anteriormente, expone que deben eliminarse las aristas que están ‘más
fuera de lugar’ por lo que se escogerá el yi que tenga menor peso y, por tanto, que proporcione
una mayor ganancia (gi = |xi | − |yi |).

Las aristas xi = (t2i−1 , t2i ) se crean de manera que t2i sea adyacente en T a t2i−1 por lo que
sin tener en cuenta ninguna restricción, tiene dos posibles opciones (a veces pudiendo estar
unı́vocamente determinada por el Criterio de Factibilidad anteriormente descrito). Para las
que ambas opciones son válidas, en la referencia [10] se explica que para x1 , x2 , x3 , el vértice
t2i se determinará aleatoriamente (en caso de mala elección se realizará backtracking sobre el
vértice no utilizado) pero a partir de x4 se escogerá la de peso máximo y, por tanto, de nuevo
que proporcione mayor ganancia.
Las aristas yi ∀i ≥ 2, se crearán con la certeza de que son capaces de “romper” un xi+1 válido,
∗
es decir, debe garantizarse la existencia de un yi+1 que cierre ciclo.
Siempre (a excepción en el backtracking del paso i = 2) debe cumplirse el criterio de
Factibilidad para encontrar intercambios secuenciales. Debe existir un yi∗ que cierre el ciclo en
caso de no poder incrementar más el valor de k. Sin embargo, ese cierre, en caso de no tener
ganancia positiva, no tiene que suponer una parada del algoritmo, simplemente no será una
solución válida.
En las dos referencias citadas, también aparece el enunciado “El algoritmo realiza un inter-
cambio 2 − opt o un 3 − opt seguido de una secuencia (vacı́a o no) de 2 − opt′ s” que de primera
mano sin justificación aparente puede resultar falsa.
El algoritmo obliga a tener la opción de cierre (con ganancia o no) en cada paso i con el ciclo
alterno x1 y1 x2 y2 . . . xi yi∗ , empezando desde i = 2 o i = 3.

Supongamos que el algoritmo termina con el cierre en el paso i + 1:

∗
ciclo alterno x1 y1 x2 y2 . . . xi yi xi+1 yi+1
aristas eliminar X = {x1 , . . . xi , xi+1 }
∗
aristas añadir Y = {y1 , . . . , yi , yi+1 }

Si se hubiera hecho el intercambio de aristas en el paso i anterior, el intercambio habrı́a sido:

28
 ciclo alterno x1 y1 x2 y2 . . . xi yi∗
aristas eliminar X = {x1 , . . . xi−1 , xi }
aristas añadir Y = {y1 , . . . , yi−1 , yi∗ }

Para llegar a la misma solución que habiendo hecho el intercambio realizado con el paso, i + 1
basta intercambiar {yi∗ , xi+1 } por {yi , yi+1
∗
}

{x1 , . . . , xi−1 , xi , xi+1 } −→ {y1 , . . . , yi−1 , yi∗ , xi+1 } −→ {y1 , . . . yi−1 , yi , yi+1
∗
}

Sin embargo, como esta planeando el algoritmo LK, no hay garantı́a de que cada uno de estos
2-intercambios (o el posible primer 3-intercambio) sean con ganancia.

29
 Ahora, se describe el algoritmo de Lin-Kernighan abordando todos los criterios y decisiones
anteriormente descritas:

Algorithm 3 Algoritmo Lin -Kernighan(LK)

INPUT : Matriz costes C , T solución factible construı́do mediante una heurı́stica constructiva
como puede ser el Vecino más próximo 2.5.1

1. i = 1 : G∗ = 0

a) Escoger el nodo t1 ∈ V de manera aleatoria.

b) x1 = (t1 , t2 ) es una de las dos arista de modo que se cumple que t2 sea adyacente en T a
t1 .
c) Ahora se debe escoger un nodo t3 (no escogido anteriormente), de modo que y1 =
(t2 , t3 ) ∈
/ T tenga peso mı́nimo y cumpliendo la condición G1 = g1 = |x1 | − |y1 | > 0. Si
no existe tal t3 , ir al paso 3 (d).

2. Sea i = i + 1. Elegir xi = (t2i−1 , t2i ) e yi de la siguiente forma:

(a) xi se escoge de manera que si t2i se une con t1 con (t2i , t1 ) ∈
/ T este camino alterno
genera un k-intercambio (se sabe que existe por el Criterio de factibilidad al haberlo
construido en el paso anterior).
Se compueba si este intercambio da un valor de ganancia positiva mejor que el que se
tiene: sea yi∗ la arista que conecta estos nodos y sea gi∗ = |xi | − |yi∗ |. Si Gi−1 + gi∗ >
G∗ =⇒ G∗ = Gi−1 + gi∗ y k = i produciendo, por tanto, un intercambio k − opt que será
utilizado en caso de no poder continuar con el proceso.
(b) Se escoge yi una arista con extremo t2i cumpliendo:
(a) Criterio de disyuntividad (no se habı́a escogido antes este vértice).
(b) Criterio de ganancia positiva (Gi > 0).
(c) El yi elegido sea capaz de romper una arista xi+1 (∃y∗i+1 ).
Se termina este paso si no existe yi o si Gi ≤ G∗ (Criterio de parada).
3. Si G∗ = 0 se aplica un backtracking limitado:

(a) Repetir el paso (2)(b) para i = 2 y buscar otro y2 tratando de incrementar su tamaño y
satisfaciendo G2 = g1 + g2 > 0.
(b) Si en ninguna opción del paso anterior se obtiene beneficio, volver al paso (2)(a) para
i = 2 y elegir otra posibilidad para x2 .
(c) Si el paso anterior también falla ir al paso (1)(c) para probar distintos y1 .
(d) Si el paso (c) no da mejora, se considera otra arista x1 en el paso (1)(b).
(e) Si (d) falla de nuevo, se selecciona otro t1 en el paso (1)(a).
El backtracking en el paso (3)(b) viola de forma temporal el Criterio de Factibilidad para
aumentar la efectividad del método.

El procedimiento termina cuando se probaron todos los n valores de t1 posibles y no se

encontró ninguna mejora.

30
4.3. Implementación en Matlab
La implementación en Matlab de este, fue lo que supuso mayor esfuerzo y gran parte del tiempo
de la realización de este trabajo. La mayorı́a de implementaciones existentes son de modificaciones
o no cumpliendo ı́ntegramente el propósito que expuso Kernighan de aumentar el k lo máximo
posible. Por ejemplo, en la referencia [10] en la implementación sustituye T por T ′ en cuanto se
encuentra ganancia sin tratar de aumentar el número de aristas a intercambiar. Aun ası́, por falta
de tiempo, quedaron puntualizaciones sin implementar.
Esquema de llamadas del programa:

function [T final,coste inicial ,coste final] = LK TSP(matrizcostes) La función an-

terior llamará de manera iterativa a lo siguiente:
function [T nuevo,G est,k,Xopt,Yopt] = LK base(T, matrizcostes)

• function E = soluc aristas(T)

• function [bool,sol] = es tour(E,X temp,Y temp)
• function coste = calcula coste(T, matrizcostes)
• function [T nuevo,X,XY,back 4,bool prov] = encontrar t4(X,E,Y,XY,M X,M Y,back 4)
◦ function [bool,sol] = es tour(E,X temp,Y temp)
⋄ function [Eord,e1,e2] = verificar exciclo(Eord,arista nueva)
• function [txi,sol ahead,Y,XY,back 5,contador 5]=
encontrar t5(X,E,Y,XY,M X,M Y,cont,back 5)
◦ function [bool,sol] = es tour(E,X temp,Y temp)
⋄ function [Eord,e1,e2] = verificar exciclo(Eord,arista nueva)
• function [Xnew,Ynew,XYnew,T nuevo,txi] =
aumentar camino(T,M X,M Y,X,Y,XY,txi,sol ahead)
◦ function [bool,sol] = es tour(E,X temp,Y temp)
⋄ function [Eord,e1,e2] = verificar exciclo(Eord,arista nueva)
Breve resumen de cada función:
LK TSP: Función que a partir de una matriz de costes de un grafo ponderado genera una
solución factible por un método constructivo para luego, tratar de mejorarla iterativamente
con el método de Lin-Kernighan.
LK base: Función principal que implementa el método de Lin-Kernighan con bracktracking,
es necesaria una solución inicial factible T y una matriz de costes asociada a un grafo.

soluc aristas: Función que a partir de la solución factible T n × 1 (donde n es el nº de

vértices/ciudades) que contiene los vértices ordenados del recorrido extrae una matriz n × 2
con las aristas que forman parte de dicho recorrido
verificar exciclo: Función que a partir de una serie de cadenas de aristas ya añadidas a la
solución, comprueba como poder “pegar” al añadir una arista nueva. Tiene como objetivo
“colocar” esa nueva arista añadida en la cadena correspondiente.

es tour: Función que comprueba la existencia de un ciclo cerrado tras la eliminación de k

aristas e introducción de otras k distintas, utilizando la función de pegado anteriormente
descrita.
encontrar t4: Función que determina el vértice t4 de x2 de modo que pueda cerrar ciclo en
caso de ser necesario y formar, por tanto y2∗ .
encontrar t5: Función que determina el vértice t5 de y2 de modo que rompa un x3 verificando
la existencia de y3∗ .

31
 aumentar camino: Función que intenta aumentar el camino alterno x1 y1 . . . xi yi en una
nueva pareja. Devuelve los conjuntos X,Y de aristas si es que ese nuevo intercambio genera
ganancia.
calcula coste: Función que a partir de una solución factible T y la matriz de costes asociada
al grafo ponderado calcula el coste del recorrido T .

Descripción de las variables empleadas:

matrizcostes : Matriz simétrica de tamaño n × n.

T final : Matriz n × 2 que contiene una solución válida
coste inicial: Valor del coste de T

coste final: Valor que se corresponde con el coste de T nuevo

T: Vector fila que es solución factible del problema, contiene los vértices en el orden en que se
visitan
T nuevo: Vector fila con la solución T mejorada tras haber intercambiado k aristas

G est: Valor de la mejora ganancia encontrada

Xopt : Matriz k × 2 que contiene las aristas a eliminar de T
Yopt : Matriz k × 2 que contiene las aristas a añadir de T
X : Matriz k × 2 que contiene en las filas las aristas a eliminar

Y : Matriz k × 2 que contiene en las filas las aristas a añadir

XY : Vector que representa el ciclo alterno de longitud 2k
E : Matriz n × 2 que contiene las aristas existentes en T

Eord : Array que contiene en las celdas una matriz k × 2, donde están concatenadas k aristas
e1,e2: Arista a evitar para que no forme ciclo esto en LK
M X : Matriz n × n que contiene los costes de las aristas que forman parte de la solución T.
El resto de elementos toman valor 0

M Y : Matriz n × n que contiene los costes de las aristas que NO forman parte de la solución
T. El resto de elementos toman valor INF
contador i : Valor que cuenta el número de t2i−1 que fueron probados
back i : Vector 1 × j que contiene las posibles opciones válidas para ti

sol ahead : Matriz n × 2 que contiene una solución válida x1 y1 x2 y2 xi yi∗

bool : Variable booleana que identifica si se encontró variable válida o no
txi : Vértice que determina X. Calculado para que yi rompa un xi+1
yi est ahead: Arista candidata a cerrar el ciclo

? ahead : Variables creadas para la comprobación de que yi rompa un xi+1 .

32
 A continuación se representa un ejemplo para tratar de ilustrar mejor este procedimiento.

4.4. Ejemplo ilustrativo para sTSP con Lin-Kernighan

Se escogió un grafo completo de 6 vértices representado por:

Figura 4.5: Grafo completo con 6 vértices

Por simplificar visualmente el grafo, no se representa el peso de cada una de las aristas, pero la
matriz de costes asociada es la siguiente:
 
∞ 1 8 6 4 3
1 ∞ 2 9 3 7
 
 8 2 ∞ 1 12 5 
C=  6 9 1 ∞ 4 13

 
 4 3 12 4 ∞ 1 
3 7 5 13 1 ∞
Se escoge de manera aleatoria una solución factible inicial:

Figura 4.6: Posible solución del problema (en rosa)

Como se puede ver en la Figura, la solución inicial es el ciclo T =[1 4 6 2 5 3] y el valor de

f (T ) = 6 + 13 + 7 + 3 + 12 + 8 = 49.

Ahora se comienza el proceso iterativo con X = [1] ,Y =[1] y XY = [1] :

1. NIVEL 1
G∗ = 0 .
Se escoge t1 = 1 y t2 = 4 =⇒ x1 = (1, 4).

33
 Ahora se elige t3 = 3, nodo no escogido, no adyacente a t2 en T y de peso mı́nimo. Luego
y1 = (4, 3).
g1 = |x1 | − |y1 | = 6 − 1 = 5 > 0
G1 = g1 = 5
X = [(1, 4)] , Y = [(4, 3)] y XY =[1 4 3]

2. NIVEL 2
Se debe probar si hay un x2 de modo que t4 no esté en XY , sea adyacente a t3 y que
(t4 , t1 ) ∈
/ T pueda formar un ciclo.
Se escoge t4 = 5 =⇒ x2 = (3, 5).
Ahora, se debe calcular el coste de cerrar el ciclo con (t4 , t1 ) en caso de que esto fuese necesario.
Sea y2∗ = (t4 , t1 ) = (5, 1) por lo que g2∗ = 12 − 4 = 8 .Se comprueba que G1 + g2∗ = 5 + 8 =
13 > 0 = G∗ y al cumplirse, se actualiza G∗ = G1 + g2∗ = 13 y k = 2 pudiéndose realizar, por
tanto, un intercambio 2 − opt.
Se escoge el vértice t5 = 6,nodo no escogido, no adyacente a t4 en T y de peso mı́nimo: se
define y2 = (5, 6). Este t5 se escoge sabiendo que existe el candidato a t6 = 2 que permite
cerrar ciclo con (t6 , t1 ) en el nivel 3.
g2 = 12 − 1 = 11
G2 = g1 + g2 = 5 + 11 = 16 > 13 = G∗ por lo que continuo el proceso.
X = [(1, 4), (3, 5)] , Y = [(4, 3), (5, 6)] y XY =[1 4 3 5 6]
3. NIVEL 3
Se define t6 = 2(determinado al escoger t5 ) =⇒ x3 = (6, 2). Continuando el proceso como en
el paso anterior, y3∗ = (2, 1) y g3∗ = 7 − 1 = 6 luego G2 + g3∗ = 16 + 6 = 22 > 13 = G∗ =⇒
G∗ = G2 + g3∗ = 22 y k = 3. El 3-intercambio es un 3 − opt.
Como tiene 6 vértices no se puede seguir y, por tanto, la elección debe ser y3∗ cerrando el ciclo
y acabando, por tanto, con un intercambio 3 − opt.

Se cambia T = T ′ con f (T ′ ) = f (T ) − G∗ = 49 − 22 = 27
X = [(1, 4), (3, 5), (6, 2)] , Y = [(4, 3), (5, 6), (2, 1)] y XY = [1 4 3 5 6 2]
Se obtiene el grafo de la Figura siguiente en donde las aristas representadas en amarillo son las
aristas añadidas:

Figura 4.7: Nueva solución T ′ =[1 3 4 6 5 2 1]

Sin entrar a explicarlo con tanto detalle, ahora se repetirı́a el proceso de igual forma desde el paso 1
con T = T ′ .
En el NIVEL 2 se tiene que X = [(1, 3), (2, 5)] , Y = [(3, 2)] y XY =[1 3 2 5], G∗ = 5 y k = 2. Al
escoger y2 , el valor de G2 = 5 = G∗ por lo que se escoge y2∗ en vez de y2 .
Los conjuntos de aristas quedarán entonces X = [(1, 3), (2, 5)] e Y = [(3, 2), (5, 1)].
Se genera, por tanto, un intercambio 2 − opt obteniendose una nueva solución T ′ con f (T ′ ) =
f (T ) − G∗ = 22 representada en el grafo siguiente:

34
 Figura 4.8: Nueva solución T ′ =[1 2 3 4 6 5 1 ]

De nuevo, tomando T = T ′ se repite el proceso. En este nuevo T comenzando por el vértice, 1

como se hizo hasta ahora, la elección de x1 está limitada, ya que si se escoge la arista x1 = (1, 2),
|x1 | = 1 y, por tanto, g1 ≤ 0 luego, x1 = (1, 5). Se sigue el proceso hasta que en el NIVEL 2 no
es posible encontrar un y2 que cumpla las condiciones. Como G∗ = 0 se realiza una técnica de
backtracking. Se escoge una nueva arista x2 (permitiendo que no cierre ciclo la arista (t4 , t1 )) hasta
conseguir después de realizar backtracking en repetidas ocasiones la solución final siguiente con
f (T ′ ) = 12:

Figura 4.9: Nueva solución T ′ =[1 2 3 4 5 6 1 ]

35
4.5. Pruebas en Matlab con implementación de Lin-Kernighan
En esta sección se hace un pequeño experimento para estudiar la eficiencia del método.

Primeramente, se hace un pequeño estudio para determinar cuál es el valor de k más utiliza-
do para tratar de probar si los 2−opt y 3−opt son los más comunes, como se expone en las referencias.

Se crean 100 matrices de costes simétricas aleatorias de tamaño 30 × 30 y se tomará T solución

factible de dos maneras distintas; creada por un método constructivo y de manera aleatoria. Te-
niendo en cuenta esto, se obtienen los resultados:

Con T aleatorio:

2 − opt: 6 veces
3 − opt: 50 veces
4 − opt: 36 veces
5 − opt: 5 veces
6 − opt: 3 veces
Con T de vecino más cercano:

2 − opt: 40 veces
3 − opt: 33 veces
4 − opt: 20 veces
5 − opt: 5 veces
6 − opt: 2 veces
A vistas de esto, se puede ver que la afirmación es cierta con el vecino más cercano, sin embargo,
escogiendo un T aleatorio son más comunes los intercambios 4 − opt que los 2 − opt.

Ahora, para determinar la eficiencia y comprobar si una solución es “buena” o no, se empleará
la fórmula de a continuación:
(coste BB − coste LK)
brecha ( %) = · 100
coste BB
de modo que cuanto mejor sea este valor, mejor será la solución obtenida por el método heurı́stico.
de modo que cuanto menor sea mejor será la solución.

Se hizo de nuevo una estudio con 10 matrices simétricas 30 × 30 determinadas de manera aletario
de modo que se aplicó el método branch and bound para poder comparar costes y tiempos de
ejecución con Lin Kernighan. Para ello se escogió la solución factible T de dos maneras distintas;
mediante la técnica del vecino más cercano y creándola de manera aleatoria. El algoritmo de Lin
Kernighan se ejecuta de manera itera mientras la brecha no supere 5 % y no supere el tiempo de
ejecucción de Branch and Bound. Los resultados obtenidos se muestran a continuación, donde la
primera tabla muestra los resultados y la segunda con una solución aleatoria del vecino más cercano:

Para ambos casos se puede ver que LK tiene un bajo tiempo de ejucción frente al de BB y la
brecha es inferior al 5 % y por tanto, bastante buena.
Las diferencias entre los modos de escoger T es notable pero aún ası́ se obtienen muy buenas
soluciones.

36
Número Matriz Coste inicial CosteLK TiempoLK CosteBB TiempoBB Brecha
1 1087 1045 1.1891561 996 6.1940975 4.9196
2 1187 1163 0.2038167 1111 22.7431762 4.68046
3 1133 1046 0.2945601 1014 8.0927479 3.1558
4 1229 1140 0.8224327 1087 5.6769973 4.8758
5 1135 1102 1.3933069 1052 14.3041019 4.7528
6 1165 1080 1.3105758 1029 6.2817252 4.95626
7 1012 1012 0.0002579 976 8.0437434 3.6885
8 999 965 0.0851515 926 9.2128799 4.2116
9 1116 1048 0.9280645 1011 10.0756143 3.6597
10 1218 1218 0.0002428 1162 28.8920664 4.819

Número Matriz Coste inicial CosteLK TiempoLK CosteBB TiempoBB Brecha

1 1087 1045 1.189156 996 6.194097 4.919679
2 1187 1163 0.203817 1111 22.743176 4.680468
3 1133 1046 0.294560 1014 8.092748 3.155819
4 1229 1140 0.822433 1087 5.676997 4.875805
5 1135 1102 1.393307 1052 14.304102 4.752852
6 1165 1080 1.310576 1029 6.281725 4.956268
7 1012 1012 0.000258 976 8.043743 3.688525
8 999 965 0.085152 926 9.212880 4.211663
9 1116 1048 0.928064 1011 10.075614 3.659743
10 1218 1218 0.000243 1162 28.892066 4.819277

37
Capı́tulo 5

Aplicación y conclusiones

En esta sección se detallan dos problemas TSP a gran escala que se resolverán mediante la
implementación en Matlab de los algoritmos presentados en los capı́tulos anteriores.
Para comprobar que los métodos implementados funcionaban adecuadamente y, por tanto, que daban
resultados correctos además de estos dos ejemplos que se muestran a continuación, se utilizaron
las instancias del TSP simétrico de la librerı́a de TSP [22] que almacena los costes de las mejores
soluciones hasta el momento y las matrices de costes asociadas.
Se mostrarán los resultados del análisis del siguiente modo:
Se obtendrá la solución óptima de los problemas junto a sus costes y tiempos de ejecución mediante
el solver ya implementado en Matlab intlinprog que se trata de un resolvedor de problemas de
Programación lineal de enteros mixtos (MILP). Se introducirá la función objetivo y las restricciones
del problema planteado en la sección 2.3.
Primeramente, se obtendrá una solución sin tener en cuenta la restricción (4) para clarificar que
esta condición es necesaria. A posteri se añadirá para comprobar su suficiencia.

El objetivo de la utilización de esta función es tener una referencia con la que poder comparar
como de buenas son las soluciones obtenidas por los otros métodos.

A continuación se expondrán los resultados de la implementación propia de los métodos.

La implementación del método de Branch and Bound se hizo con backtracking por lo que se
sabe de antemano que aunque el tiempo de ejecución sea mayor se obtendrá la solución óptima.

El método de Lin Kernighan aunque sea heurı́stico, se comprobó su eficiencia y gran aproximación
a la solución óptima en 4.5.

5.1. Ejemplo práctico TSP Euclı́deo

El primero de los problemas planteados considera 38 ciudades situadas en el mapa de Galicia
distribuidas según el mapa siguiente:

38
 Figura 5.1: Mapa de Galicia con las ciudades escogidas

Para este caso primer ejemplo, no se tendrá en cuenta la red de carreteras, sino que se va a
suponer que todas las ubicaciones están conectadas entre sı́ por una arista. Se trabajará, por tanto,
desde la perspectiva del TSP Euclı́deo.

Definición 5.1.1. Sea G = (V, A) un grafo completo con |V | = n y C = [cij ] una matriz simétrica,
se llama TSP métrico si se cumple la desigualdad triangular cik ≤ cij + cjk ∀i, j, k ≤ n.
Definición 5.1.2. Se llama TSP Euclı́deo al problema que considera n puntos en R2 con
distancia euclı́dea, es decir, cij = d(x, y) = ∥x − y∥2 donde (x, y) son las coordenadas geográficas
(latitud,
P longitud) y se debe encontrar el ciclo hamiltoniano con peso mı́nimo, es decir,
(x,y)∈T ∥x − y∥2 sea lo mı́nimo posible. El TSP Euclı́deo es simétrico y métrico.

Supongamos entonces que nos interesa recorrer un camino que una todas estas ciudades comen-
zando y acabando en el mismo lugar, de modo que se minimice la distancia recorrida.

Para la construcción de la matriz de costes, es necesario el peso de las aristas y, por tanto,
distancias entre las ciudades se ha calculado a partir de las coordenadas geográficas (x, y) =(latitud,
longitud) mediante la Fórmula de Haversine.[1]

La fórmula de Harversine es la siguiente:

a = sin2 ( φ2 −φ1
) + cos(φ1 ) · cos(φ2 ) · sin2 ( λ2 −λ 1
)
2 √ √ 2
c = 2 · atan2( a, 1 − a)
d=R·c (5.1)
Siendo la terminologı́a empleada:

R: radio de la Tierra (6.371 km)

φ1 y φ2 : latitudes en notación decimal
λ1 y λ2 : longitudes en notación decimal

Teniendo en cuenta todo esto y siguiente el esquema descrito en el comienzo de esta sección:

39
5.1.1. Solver intlinprog
Para este caso se determinan las soluciones al problema desde las dos perspectivas citadas.
No teniendo en cuenta la restricción de subciclos se obtendrı́a lo siguiente:

Figura 5.2: Subciclos que se pueden formar sobre el

grafo
Figura 5.3: Subciclos sobre el mapa

Sin embargo, una vez añadida la restricción se obtiene el ciclo hamiltoniano representado:

Figura 5.4: Solución óptima en grafo

Figura 5.5: Solución sobre mapa

5.1.2. Branch and Bound

Teniendo en cuenta lo discutido en la Sección 3 ası́ como su implementación, se hará una primera
búsqueda sin backtracking y a continuación se volverá sobre los nodos vivos para ramificar desde ellos.

La primera solución sin aplicar backtracking con coste 1.1827e+03 y tiempo de ejecución 0,0628
segundos es la siguiente:

40
 33 25 17 37 2 4 8 32 11 19 15 14 12 34 13 23 24 16 21
20 7 18 1 26 36 38 22 28 31 5 3 27 27 23 24 30 6 35

Aplicando backtracking mejora considerablemente el coste obteniendo ası́ la solución óptima:

13 34 27 19 31 15 12 3 14 5 4 36 8 28 11 22 32 38 20
33 7 25 26 1 2 37 18 17 10 9 29 35 6 21 30 16 24 23

5.1.3. Vecino más cercano

El método heurı́stico voraz da un coste de solución igual a 1.1949e+03 y la solución obtenida es:

1 2 37 18 17 7 25 26 36 4 15 3 5 14 12 31 19 28 11
8 38 32 22 23 13 34 27 24 16 30 21 6 35 29 9 10 33 20

5.1.4. Lin-Kernighan
En este caso, como se detalló en la sección de este método, se empleará el método constructivo
del vecino más cercano para tratar de mejorar la solución obtenida mediante una serie de k − opt.
Se obtiene como solución el ciclo que se muestra a continuación cuyo coste es 1,0760e + 03 difiriendo
tan solo un 1,3374 % de la solución y con un tiempo de ejecución de 11,6952.

32 20 33 7 25 26 1 2 37 18 17 10 9 29 35 6 21 30
16 24 13 23 22 34 27 19 31 15 12 3 14 5 4 36 8 28

En la siguiente tabla se compactan todos los resultados de coste y tiempo de ejecución para
comparar todas las técnicas a la vez:

L1: intlinprog con tres restricciones

L2: intlinprog del problema completo
VC: Vecino más cercano
BB: Branch and Bound

LK: Lin Kernighan

L1 L2 VC BB LK
Coste 1.049e+03 1.0618e+03 1.1948e+03 1.0618e+03 1.0760e+03
Tiempo de ejecución(seg) 0.8095 0.5628 0.0083 34.9388 11.6952

Tabla 5.1: Resultados de Coste y Tiempo de ejecución

41
5.2. Ejemplo práctico TSP
Este segundo ejemplo tratará de una manera más “realista” el problema. El grafo estará consti-
tuido por 18 nodos, 18 ciudades situadas en el mapa de Galicia conectadas entre sı́ por carretera.

Vértice Ubicación
1 Ferrol
2 A Coruña
3 Castrelo de Miño
4 Santiago de Compostela
5 Lugo
6 Ribadeo
7 Viveiro
8 Finisterre
9 Puebla de Trives
10 Tui
11 Villalba
12 Monforte de Lemos
13 Betanzos
14 Lalı́n
15 Guitiriz
16 Arzúa
17 Curtis
18 Cerceda

Supongamos que de nuevo nos interesa recorrer un camino que una todas estas ciudades comenzando
y acabando en el mismo lugar, de modo que se minimice los kilómetros recorridos.

Se tratará, por tanto, de un grafo ponderado donde las aristas serán las propias carreteras y los
pesos de cada una de ellas se corresponden con las distancias en km tomando de referencia la red
de carreteras del IGN y con la ayuda de Google Maps [17] escogiendo de ser posible la de menor
kilometraje, sin peajes y por autopista.
La ubicación de las ciudades en el mapa es la siguiente:

Figura 5.6: Mapa de Galicia con las ciudades escogidas creada con Matlab

42
 Aunque las carreteras no unan ciudades en lı́nea recta, se simplifica el problema para represen-
tarlo en un grafo de esta manera. Las aristas de este problema en caso de ser un grafo completo
serı́an 18∗17
2 = 153 sin embargo, la figuras siguientes muestran las aristas a tener en cuenta:

Figura 5.7: Mapa con 18 vértices y aristas Figura 5.8: Grafo con 18 vértices y aristas
existentes existentes

Para la implementación, las aristas que no existen tomarán un valor muy elevado (máximo valor
·1000) para que no sean escogidas. Es por esto que técnicas con la del vecino más cercano no tiene
sentido aplicarla en este caso.

5.2.1. Solver intlinprog

Primeramente, como se dijo anteriormente, se aplica intlinprog y se obtienen las soluciones
representadas a continuación.

En caso de no tener en cuenta la restricción de ruptura de subciclos se obtendrı́a la misma que

añadiéndola, ya que no se forman ciclos disjuntos de menor longitud en la solución. El coste de la
solución óptima representada es 1015.

Figura 5.10: Subciclos que se pueden formar sobre

el mapa
Figura 5.9: Subciclos que se pueden for-
mar en un grafo

43
5.2.2. Branch and Bound
Al igual que para el ejemplo anterior, primero se aplicó el algoritmo sin retroceder en la búsqueda
y a continuación se estudiaron los nodos vivos sin utilizar.
Tras haber aplicado backtracking el coste de la solución que se muestra a continuación, disminuye
considerablemente con un tiempo de 0,7573 segundos a diferencia de los 0,0612 segundos del paso
anterior.
6 7 1 13 2 18 8 4 10
3 9 12 14 16 17 15 5 11

5.2.3. Lin-Kernighan
Para este caso, como el coste de una solución aleatoria o tras aplicar la técnica del vecino más
cercano es muy elevado, se tratará de mejorar el “branching a derecha” de Branch and Bound, con
este método. Con esto, se obtiene la solución siguiente:
6 11 15 5 14 12 9 3 10
4 8 18 2 17 16 13 1 7
Esta solución es óptima, pues su coste coincide con el determinado al principio de esta sección. El
tiempo de ejecución es demasiado elevado para ser esta técnica, pero se debe a que aproximo del
todo la solución hasta hacerla óptima.
En la siguiente tabla se mostrarán una recopilación con los costes de las soluciones óptimas
obtenidas, ası́ como los tiempos de ejecución, donde la nomenclatura es la siguiente:

L1: intlinprog con tres restricciones

L2: intlinprog del problema completo

BB: Branch and Bound

LK: Lin Kernighan

L1 L2 BB LK
Coste 1015 - 1015 1015
Tiempo de ejecución(seg) 0.1203 - 0.1416 1.1822

Tabla 5.3: Resultados de Coste y Tiempo de ejecución

En conclusión para ambos ejemplos, como se expuso en los capı́tulo correspondiente, el método
heurı́stico Lin Kernighan da muy buenos resultados en comparación con la solución óptima a pesar
de ser muy modificado por su “baja eficiencia”.

44
Bibliografı́a

[1] ¿Cómo calcular la distancia entre dos puntos geográficos?(Fórmula de Haversine). url:
https : / / www . genbeta . com / desarrollo / como - calcular - la - distancia - entre - dos -
puntos-geograficos-en-c-formula-de-haversine.
[2] José Vicente Álvarez. TEMA 5:COMPLEJIDAD ALGORÍTMICA. url: https://www2.
infor.uva.es/~jvalvarez/docencia/tema5.pdf.
[3] Clay Mathematics Institute. P vs. NP. url: https://www.claymath.org/millennium/p-
vs-np/.
[4] Ismael Crehuet Lucas. ((El problema del viajante con grafos)). Universidad de Valladolid. url:
https://uvadoc.uva.es/bitstream/handle/10324/57985/TFG-G5977.pdf?sequence=1.
[5] Definición de Grafos. url: https://ccia.ugr.es/ ~jfv/ed1/tedi/cdrom/docs/grafos.
htm.
[6] Reinhard Diestel. Graph Theory. Graduate Texts in Mathematics 173. Springer-Verlag, 2006.
[7] Mario Fioravanti. Apuntes de Matemática Discreta. 2021-2022. url: https://personales.
unican.es/fioravam/MatDiscreta/Mat-Discreta_3ra-parte_nov2021.pdf.
[8] R.Fulkerson y S.Johson G.Dantzing. Solution of a large scale traveling salesman problem. 12
April 1954. url: https://www.rand.org/content/dam/rand/pubs/papers/2014/P510.
pdf.
[9] Keld Helsgaun. An Effective Implementation of the Lin-Kernighan Traveling Salesman Heu-
ristic. Inf. téc. Department of Computer Science, Roskilde University, Denmark, 2000. url:
http://www.akira.ruc.dk/~keld/research/LKH/.
[10] Keld Helsgaun. ((An effective implementation of the Lin-Kernighan traveling salesman heuris-
tic)). En: European Journal of Operational Research (2000).
[11] Keld Helsgaun. ((General k-opt submoves for the Lin–Kernighan TSP heuristic)). En: Mathe-
matical Programming Computation 1.2-3 (2009), págs. 119-163.
[12] Miguel Infantes Durán. ((Trabajo de Fin de Grado: El Problema del Viajante (TSP))). Tesis
de mtrı́a. Universidad de Sevilla (US).
[13] J. Lee. Imagen puentes de Konigsberg. 2020. url: https://www.ksakosmos.com/post/%EA%
B7%B8%EB%9E%98%ED%94%84%EC%9D%B4%EB%A1%A0- %EC%84%B8%EC%83%81%EC%9D%84-
%EB%B0%94%EA%BE%B8%EB%8B%A4.
[14] S. Lin y B. W. Kernighan. ((An Effective Heuristic Algorithm for the Traveling Salesman
Problem)). En: Bell Telephone Laboratories, Incorporated, Murray Hill, N.J. (1971). 15 October
1971.
[15] S. Lin y B. W. Kernighan. ((An effective heuristic algorithm for the traveling-salesman
problem)). En: Operations Research 21.2 (1973), págs. 498-516. url: http://www.jstor.org/
stable/169020.
[16] John D. C. Little et al. ((An Algorithm for the Traveling Salesman Problem)). En: Operations
Research (6 March 1963).
[17] Google Maps. Google Maps. url: https://www.google.es/maps.
[18] Cecilia Pola Méndez. Apuntes de Optimización I. 2022-2023.

45
[19] Milestones in the History of the TSP. 2005. url: https://www.math.uwaterloo.ca/tsp/
history/milestone.html.
[20] Universidad del Paı́s Vasco (EHU). Investigación Operativa. Programación Lineal. Tema 5.
Programación entera. url: https://ocw.ehu.eus/file.php/19/5._entera.pdf.
[21] Problema de los Puentes de Königsberg. url: https://es.wikipedia.org/wiki/Problema_
de_los_puentes_de_K%C3%B6nigsberg.
[22] Gerhard Reinelt. TSPLIB 95. 1995.
[23] Alexander Schrijver. A Course in Combinatorial Optimization. 1982. url: https://homepages.
cwi.nl/~lex/files/dict.pdf.
[24] Tema 8:Complejidad Temporal. url: https://www.uhu.es/francisco.moreno/gii_mac/
docs / Tema _ 8 . pdf# : ~ : text = Se % 20define % 20la % 20complejidad % 20temporal % 20de %
20Mcomo%20una,M%20para%20ejecutar%20una%20entrada%20de%20longitud%20n.
[25] Dafne Lucı́a Arias Vilaboa. ((Modelos y algoritmos para el problema del viajante. Una
aplicación en planificación socio sanitaria)). Universidade de Santiago de Compostela (USC),
2021.

46
 Apéndice A

Programas en Matlab

En este apéndice, se incluyen los programas necesarios para el cálculo de soluciones (BB y LK).
El resto se adjuntan en un archivo .zip a la memoria.

A. Vecino más cercano

1 function [ coste , tour ] = v e c i n o m a s p r o x i m o ( matrizcostes )
2
3 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
4 % Funcion que implementa el algoritmo de Vecino mas cercano
5 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
6 % INPUTS :
7 % matrizcostes : Matriz simetrica nxn con los costes entre ciudades
8 % OUTPUTS :
9 % tour : Vector fila con la mejor solucion encontrada
10 % coste : Valor del coste de la mejor solucion
11 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
12
13 %Inicializo la salida y variables necesarias
14
15 n = size ( matrizcostes , 1) ;
16 mejor_coste = Inf ;
17 mejor_tour = [];
18
19 % Intentar iniciar desde cada ciudad
20 for inicio = 1: n
21 % In icializa cion
22 tour = zeros (1 , n ) ;
23 ciud _visitad as = zeros (1 , n ) ;
24 ciudad = inicio ;
25 tour (1) = ciudad ;
26 ciud _visitad as ( ciudad ) = 1;
27 for i = 2: n
28 % Inicializo variables de test de parada
29 min_coste = Inf ;
30 c i u d a d _ m a s _ c e r c a n a = 0;
31 for j = 1: n
32 % Se debe encontrar la ciudad mas cercana ( con menor coste )
33 que no fue visitada
34 if ci ud_visit adas ( j ) ==0 && matrizcostes ( ciudad , j ) < min_coste
35 min_coste = matrizcostes ( ciudad , j ) ;
36 ciudad_mas_cercana = j;
37 end
38 end
39
40 % Si no se encuentra una ciudad mas cercana , romper el bucle
41 if c i u d a d _ m a s _ c e r c a n a == 0
42 break ;
43 end

47
44
45 tour ( i ) = c i u d a d _ m a s _ c e r c a n a ;
46 ciud_ visitad as ( c i u d a d _ m a s _ c e r c a n a ) = 1;
47 ciudad = c i u d a d _ m a s _ c e r c a n a ;
48 end
49
50 % Verificar si se ha completado un tour
51 if all ( ciud _visita das )
52 coste = calcula_coste ( tour , matrizcostes ) ;
53 if coste < mejor_coste
54 mejor_coste = coste ;
55 mejor_tour = tour ;
56 end
57 end
58 end
59 % Si no se encuentra un tour valido , retornar vacio
60 if isempty ( mejor_tour )
61 coste = Inf ;
62 tour = [];
63 else
64 coste = mejor_coste ;
65 tour = mejor_tour ;
66 end
67 end

B. Branch and Bound

B.1. BB TSP
1 function [ mejorsol , costesol , tiempo_branch , tiempo ] = BB_TSP ( matrizcostes )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que aplica el m t o d o de Branch and Bound hasta conseguir una
4 % soluci n ptima y su coste asociado
5 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
6 % INPUTS :
7 % matrizcostes : Matriz s i m t r i c a de t a m a o nxn
8 % OUTPUTS :
9 % mejorsol : Vector fila con las aristas que e s t n incluidas en la
10 % soluci n
11 % costesol : Valor con el coste de la s o l u c i n
12 % tiempo_branch : Valor en segundos del ‘‘ branching a derecha ’’
13 % tiempo : Valor en segundos co rr es p on di en t e al tiempo de e j e c u c i n
14 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
15 % Inicializamos la salidas y las variables necesarias
16 tic
17 n = size ( matrizcostes , 1) ; % N m e r o de v r t i c e s que s e r n las ciudades
18 % NODO RAIZ
19 % Crear matriz reducida en filas y columnas
20 [ matrizred , costered ] = reducematriz ( matrizcostes ) ;
21 % Coste inferior del nodo r a z y por tanto del problema
22 LB = costered ;
23 % RAMIFICACION A DERECHA
24 [ LB , ~ , Eord , LB_I_lista , matriz_reducida_lista , E_lista , Eord_lista ] = BB_base (
matrizred , LB ) ;
25 Z = LB ; % Cota superior
26 tiempo_branch = toc ;
27 %SE INICIA BACKTRACKING
28 %Emplear la t c n i c a de backtraking si fuese necesario
29 i =1;
30 while i < length ( m a t r i z _ r e d u c i d a _ l i s t a ) +1
31 %Eliminar del array los que no cumplan la condicion
32 if LB_I_lista ( i ) >= Z
33 m a t r i z _ r e d u c i d a _ l i s t a ( i ) =[];
34 LB_I_lista ( i ) =[];
35 E_lista ( i ) =[];
36 Eord_lista ( i ) =[];
37 else
38 i = i +1;
39 end

48
40 end
41 % Se ramifica a derecha en los nodos que es posible
42 while ~ isempty ( m a t r i z _ r e d u c i d a _ l i s t a )
43 %Encontrar el minimo de los no vacios en LB_I_lista
44 [ LB_izq , idx_min ] = min ( LB_I_lista ) ;
45 matriz = m a t r i z _ r e d u c i d a _ l i s t a { idx_min };
46 % Reduzco la matriz para poder ramificar
47 [ matrizred , ~] = reducematriz ( matriz ) ;
48 LB = LB_izq ;
49 % Ramificar sobre esa matriz reducida a derecha
50 [ LB_temp , E_temp , Eord_temp , LB_I_lista_temp , matriz_reducida_lista_temp ,
E_lista_temp , E or d _l is ta _ te mp ] = BB_base ( matrizred , LB ,Z , E_lista { idx_min } ,
Eord_lista { idx_min }) ;
51 %Se elimina ese nodo vivo de las listas ya que fue ramificado
52 m a t r i z _ r e d u c i d a _ l i s t a ( idx_min ) =[];
53 LB_I_lista ( idx_min ) =[];
54 E_lista ( idx_min ) =[];
55 Eord_lista ( idx_min ) =[];
56 %Se concatenan las listas
57 m a t r i z _ r e d u c i d a _ l i s t a = [ matriz_reducida_lista , m a t r i z _ r e d u c i d a _ l i s t a _ t e m p ];
58 LB_I_lista = [ LB_I_lista , L B_ I_ li s ta _t em p ];
59 E_lista =[ E_lista , E_lista_temp ];
60 Eord_lista =[ Eord_lista , E or d_ li s ta _t em p ];
61 % Se comprueba si la ramificacion a derecha de ese nodo vivo
62 % dio lugar a una s o l u c i n completa y con menor coste ,
63 % ac tualizan dola si es el caso
64 if size ( E_temp ,1) == n && LB_temp < Z
65 Z = LB_temp ;
66 Eord = Eord_temp ;
67 i =1;
68 % Eliminar de la lista nodos vivos tras la actualizacion
69 % del coste
70 while i < length ( m a t r i z _ r e d u c i d a _ l i s t a ) +1
71 %Eliminar del array los que no cumplan la condicion
72 if LB_I_lista ( i ) >= Z
73 m a t r i z _ r e d u c i d a _ l i s t a ( i ) =[];
74 LB_I_lista ( i ) =[];
75 E_lista ( i ) =[];
76 Eord_lista ( i ) =[];
77 else
78 i = i +1;
79 end
80 end
81 end
82 end
83 costesol = Z ;
84 mejorsol = Eord {1};
85 tiempo = toc ;
86 end

B.2. BB base
1 function [ LB , E , Eord , LB_I_lista , matriz_reducida_lista , E_lista , Eord_lista ] =
BB_base ( matrizred , LB ,Z ,E , Eord )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % Funcion que ramifica sobre los nodos a derecha . En primer lugar , se
4 % empleara para obtener una solucion factible pero no necesa riamente
5 % optima .
6 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
7 % INPUTS :
8 % matrizred : Matriz simetrica reducida de t a m a o nxn
9 % LB : Cota inferior de v ( P )
10 % Z : Cota superior de v ( P )
11 % E : Matriz nx2 que contiene las aristas incluidas en la solucion
12 % Eord : Array que contiene en cada celda una matriz kx2 , donde estan
13 % concatenadas k aristas
14 % NOTA : Si nargin ==2 solo son necesarios los dos primeros inputs ya que se
15 % encuentra en el " branching a derecha " del algoritmo
16 % Si nargin >2 se tiene una solucion factible y por tanto una cota
17 % superior Z ademss de una solucion parcial con

49
18 % aristas incluidas E y dichas aristas ordenadas Eord
19 % ( a excepcion del nivel i =1 que estaran vacios )
20 % OUTPUTS :
21 % LB : Cota inferior de la solucion factible obtenida
22 % E : Matriz nx2 que contiene las aristas incluidas en la solucion
23 % Eord : Array con una celda que contiene las aristas ordenadas de
24 % la solucion
25 % Listas para almacenar informacion de nodos " a izquierda " de este
26 % " branching a la derecha "
27 % LB_I_lista : Vector fila que almacena las cotas inferiores
28 % c o r r e s p o n di e n t e s a las soluciones parciales de los
29 % nodos de la izquierda
30 % m a t r i z _ r e d u c i d a _ l i s t a : Array con las matrices ya reducidas de
31 % los nodos de la izquierda
32 % E_lista : Array con las aristas que se van incluyendo para crear
33 % la solucion parcial por tanto , en el orden de inclusion
34 % Eord_lista : Array con las aristas de E ordenadas en cadenas
35 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
36 %Inicializo la salida para almacenar la informacion de los nodos a
37 %izquieda
38 LB_I_lista = [];
39 m a t r i z _ r e d u c i d a _ l i s t a = {};
40 E_lista ={};
41 Eord_lista ={};
42
43 % Almacenar las aristas que estan incluidas en el ciclo
44 if nargin ==2 %Todavia no inclui ninguna arista
45 E = []; % aristas incluidas
46 Eord = {};
47 end
48 % Bucle iterativo hasta se reduzca la matriz a 2 x2 con elementos no inf
49 while sum ( matrizred (:) ~= Inf ) > 2
50 % Calcular la penalizacion de los 0 para ver las elecciones posibles
51 [ max_penalizacion , fila , columna ]= c a l c u l a _ p e n a l i z a c i o n ( matrizred ) ;
52 LB_I = LB + m a x _ p e na l i z a c i o n ; %cota inferior nodo izquierda
53 %si el nodo a izquierda creado es vivo por tener menor LB se a a d e
54 if nargin ==2 || LB_I <Z
55 %Se reduce la matriz del nodo de la izquierda
56 matriznueva = matrizred ;
57 matriznueva ( fila , columna ) = Inf ;
58 matriznueva = reducematriz ( matriznueva ) ;
59 % Se a a d e a las listas la informacion
60 LB_I_lista = [ LB_I_lista , LB_I ];
61 m a t r i z _ r e d u c i d a _ l i s t a { end + 1} = matriznueva ;
62 E_lista { end +1}= E ;
63 Eord_lista { end +1}= Eord ;
64 end
65 % A a d i r la arista seleccionada a la lista de aristas incluidas
66 E = [ E ; fila , columna ];
67 % Para RAMIFICAR POR LA DERECHA
68 % Establecer infinito para la fila y columna seleccionadas para
69 % que no se puedan escoger elementos de estas
70 matrizred ( fila , :) = Inf ;
71 matrizred (: , columna ) = Inf ;
72 % A a d i d o manual para la no formacion de ciclos de menor longitud
73 matrizred ( columna , fila ) = Inf ; %evitar ciclo con la arista nueva a a d i d a
74 [ Eord , e1 , e2 ] = v e r i f i c a r _ e x c i c l o ( Eord ,[ fila , columna ]) ;
75 matrizred ( e1 , e2 ) = Inf ;
76 % Actualizacion de matrizred
77 [ matrizred , costered ] = reducematriz ( matrizred ) ;
78 % Coste
79 LB = LB + costered ;
80 if nargin > 2 && LB >= Z
81 return
82 end
83 end
84 % Finalmente , cuando queden dos valores distintos de infinito , se
85 % a a d e n ambos
86 [ filas , columnas ] = find ( matrizred ~= Inf ) ;
87 for i = 1: length ( filas )

50
88 E = [ E ; filas ( i ) , columnas ( i ) ];
89 Eord = v e r i f i c a r _ e x c i c l o ( Eord ,[ filas ( i ) , columnas ( i ) ]) ;
90 LB = LB + matrizred ( filas ( i ) , columnas ( i ) ) ;
91 end
92 end

B.3. reducematriz
1 function [ matrizred , costered ] = reducematriz ( matrizcostes )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que obtiene una matriz reducida en filas y columnas
4 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
5 %INPUTS :
6 % matrizcostes : matriz s i m t r i c a de t a m a o nxn
7 %OUTPUTS :
8 % matrizred : matriz s i m t r i c a reducida de t a m a o nxn
9 % costered : valor del coste de la r e d u c c i n
10 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
11 %Inicializamos los elementos de salida
12 n = size ( matrizcostes ,1) ;
13 matrizred = matrizcostes ;
14 %FILAS
15 u = min ( matrizcostes ,[] ,2) ;
16 for i = 1: n
17 % Si toda la fila es infinito , establecer ui en 0
18 if all ( isinf ( matrizcostes (i , :) ) )
19 u ( i ) = 0;
20 else
21 matrizred (i , :) = matrizred (i , :) - u ( i ) ;
22 end
23 end
24 %COLUMNAS
25 v = min ( matrizred , [] , 1) ;
26 for j = 1: n
27 % Si toda la columna es infinito , establecer vj en 0
28 if all ( isinf ( matrizred (: , j ) ) )
29 v ( j ) = 0;
30 else
31 matrizred (: , j ) = matrizred (: , j ) - v ( j ) ;
32 end
33 end
34 costered = sum ( u ) + sum ( v ) ;
35 end

B.4. calcula penalización

1 function [ max_penalizacion , fila , columna ]= c a l c u l a _ p e n a l i z a c i o n ( matriz )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que calcula las pe nalizac iones de la matriz reducida para saber
4 % que cero escoger y por tanto por que arista ramificar .
5 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
6 %INPUTS :
7 % matriz : matriz s i m t r i c a reducida de t a m a o nxn
8 %OUTPUTS :
9 % m a x _ p e n a l i z a c i o n : valor de la p e n a l i z a c i n m x i m a y por tanto la
10 % elegida
11 % fila : ndice de la fila de la matriz que contiene el 0 con mayor
12 % penalizacion
13 % columna : ndice de la columna de la matriz que contiene el 0 con
14 % mayor penalizacion
15 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
16 n = size ( matriz ,1) ;
17 penalizacion = zeros ( n ) ; %Matriz en donde los elementos ~= 0 de la matriz de
costes t e n d r n valor 0
18 % y los elementos = 0 la suma del elemento minimo de
la fila m s el elemento minimo de la columna
19 minimo_fila = zeros ( n ) ; % Matriz c o n s t r u d a de igual forma que la
p e n a l i z a c i n pero en los elementos = 0 s e r el m n i m o
20 % la fila que corresponda

51
21 mini mo_colum na = zeros ( n ) ; % Matriz c o n s t r u d a de igual forma que la
p e n a l i z a c i n pero en los elementos = 0 s e r el minimo
22 % la columna que corresponda
23
24 for i = 1: n %indice fila
25 for j =1: n %indice columna
26 if matriz (i , j ) ==0
27 minimo_fila (i , j ) = min ( matriz (i , setdiff (1: n , j ) ) ) ;
28 mini mo_colum na (i , j ) = min ( matriz ( setdiff (1: n , i ) ,j ) ) ;
29 penalizacion (i , j ) = minimo_fila (i , j ) + min imo_colu mna (i , j ) ;
30 end
31 end
32 end
33 [ max_penalizacion , ~] = max ( penalizacion (:) ) ;
34 % Indices del 0 con la p e n a l i z a c i n m s alta
35 [ fila , columna ] = find ( penalizacion == m a x _ p e n al i z a c i o n ) ;
36
37 % Si el vector fila tiene m s de un elemento , buscar el m x i m o de la
p e n a l i z a c i n de la fila
38 if length ( fila ) > 1
39 % Inicializo para test de parada
40 m a x _ m i n i m o _ f i l a _ v a l = - inf ;
41 mejor_ind = 0;
42
43 for k = 1: length ( fila )
44 if minimo_fila ( fila ( k ) , columna ( k ) ) > m a x _ m i n i m o _ f i l a _ v a l
45 m a x _ m i n i m o _ f i l a _ v a l = minimo_fila ( fila ( k ) , columna ( k ) ) ;
46 mejor_ind = k ;
47 end
48 end
49 % Seleccionar la fila y columna c o r r e s p o n d i e n t e s
50 fila = fila ( mejor_ind ) ;
51 columna = columna ( mejor_ind ) ;
52 else
53 % Seleccionar el nico valor encontrado
54 fila = fila (1) ;
55 columna = columna (1) ;
56 end
57 end

B.5. verificar exciclo

1 function [ Eord , e1 , e2 ] = v e r i f i c a r _ e x c i c l o ( Eord , arista_nueva )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que obtiene un array cuyas celdas contienen las cadenas disjuntas
4 % que forman las aristas a a d i d a s hasta ahora a la s o l u c i n . El objetivo
5 % es " colocar " la arista a a d i d a en la cadena c or r es po nd i en te ( uniendo dos
6 % cadenas si se diese el caso ) a d e m s de reconocer los extremos de los
7 % nuevos ciclos de longitud menor a n que se deben evitar .
8 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
9 % INPUTS :
10 % Eord : Array que contiene las cadenas tras ir ‘ pegando ’ las aristas
11 % a a d i d a s hasta el momento
12 % arista_nueva : Vector 2 x1 que se corresponde con la arista que es
13 % a a d i d a a la s o l u c i n en cada paso
14 % OUTPUTS :
15 % Eord : Array que contiene una arista m s que el de entrada en donde
16 % cada celda es una matriz kx2 , donde e s t n concatenadas k aristas
17 % e1 : Valor del v r t i c e de un extremo del recorrido formado para
18 % evitar ciclo
19 % e2 : Valor del v r t i c e del otro extremo del recorrido formado para
20 % evitar ciclo
21 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
22 %Se designa cada v r t i c e de la arista que se acaba de a a d i r
23 i0 = arista_nueva (1) ;
24 j0 = arista_nueva (2) ;
25 %Para la primera i t e r a c i n , se pone el valor de la arista entrante
26 if isempty ( Eord )
27 Eord {1}=[ i0 , j0 ];
28 e1 = j0 ;

52
29 e2 = i0 ;
30 return ;
31 end
32
33 %Si es una arista que cierra ciclo hamiltoniano
34 if length ( Eord ) ==1
35 if j0 == Eord {1}(1 ,1) && i0 == Eord {1}( end ,2) %MONI c a m b i a " end "
36 Eord {1}=[ Eord {1};[ i0 , j0 ]];
37 e1 = [];
38 e2 = [];
39 return
40 end
41 end
42
43 % C o n c a t e n a c i n a izquierda de la nueva arista introducida
44 for s =1: length ( Eord )
45 if j0 == Eord { s }(1 ,1)
46 Eord { s }=[ i0 , j0 ; Eord { s }];
47
48 % C o m p r o b a c i n de que no une dos cadenas
49 % C o n c a t e n a c i n a derecha con otro recorrido ya existente
50 for r =1: length ( Eord )
51 if i0 == Eord { r }( end ,2)
52 Eaux = [ Eord { r }; Eord { s }];
53 Eord { s } = Eaux ;
54 Eord ( r ) = [];
55 e1 = Eaux ( end ,2) ;
56 e2 = Eaux (1 ,1) ;
57 return ;
58 end
59 end
60 e1 = Eord { s }( end ,2) ;
61 e2 = i0 ;
62 return ;
63 end
64 end
65
66 % C o n c a t e n a c i n a derecha de la nueva arista introducida
67 for k =1: length ( Eord )
68 if i0 == Eord { k }( end ,2)
69 Eord { k }=[ Eord { k }; i0 , j0 ];
70 e1 = j0 ;
71 e2 = Eord { k }(1 ,1) ;
72 return ;
73 end
74 end
75 %Si no c o n c a t e n en n i n g n lado
76 Eord { length ( Eord ) +1}=[ i0 , j0 ];
77 e1 = j0 ;
78 e2 = i0 ;
79 end

C. Lin-Kernighan
C.1. LK TSP
1 function [ T_final , coste_inicial , coste_final ] = LK_TSP ( matrizcostes )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %5 % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que aplica de manera iterativa LK_base hasta conseguir la mejor
4 % soluci n
5 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
6 % INPUTS :
7 % matrizcostes : Matriz s i m t r i c a nxn que tiene los costes de ir de
8 % una ciudad a otra
9 % OUTPUTS :
10 % T_final : Vector fila con la mejor s o l u c i n T ’ encontrada con
11 % este m t o d o
12 % coste_inicial : Valor del coste de la s o l u c i n T
13 % coste_final : Valor del coste de la s o l u c i n final T_final

53
14 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
15
16 %Inicializo la salida
17 [ coste_inicial , T ] = v e c i n o m a s p r o x i m o ( matrizcostes ) ;
18
19 [ T_nuevo , G_est ,k , Xopt , Yopt ] = LK_base (T , matrizcostes ) ;
20 if G_est == 0
21 T_final = T ;
22 coste_final = calcula_coste ( T_final , matrizcostes ) ;
23 else
24 while G_est > 0
25 %toc
26 T = T_nuevo ;
27 [ T_nuevo , G_est ,k , Xopt , Yopt ] = LK_base (T , matrizcostes ) ;
28 end
29 T_final = T_nuevo ;
30 coste_final = calcula_coste ( T_final , matrizcostes ) ;
31 end
32 end

C.2. LK base
1 function [ T_nuevo , G_est ,k , Xopt , Yopt ] = LK_base (T , matrizcostes )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %5 % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que realiza un movimiento k - opt sobre T s o l u c i n factible
4 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
5 % INPUTS :
6 % T : Vector fila que es s o l u c i n factible del problema , contiene los
7 % v r t i c e s en el orden en que se visitan
8 % matrizcostes : Matriz s i m t r i c a nxn que tiene los costes de ir de
9 % una ciudad a otra
10 % OUTPUTS :
11 % T_nuevo : Vector fila con la s o l u c i n T mejorada tras haber
12 % intercambiado k aristas
13 % G_est : Valor de la mejora ganancia encontrada
14 % k : Valor del n m e r o de aristas interca mbiadas
15 % Xopt : Matriz kx2 que contiene las aristas a eliminar de T
16 % Yopt : Matriz kx2 que contiene las aristas a a a d i r de T
17 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
18 % Inicialico las variables necesarias y de salida
19 n = length ( T ) ;
20 E = soluc_aristas ( T ) ; % Matriz nx2 que contiene las aristas del tour T
21 G_est = 0; % Mejor ganancia hasta el momento
22 k =0;
23 Xopt =[];
24 Yopt =[];
25 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
26 % Creamos variables M_X y M_Y
27 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
28 % Crear una matriz en donde solo aparezcan los costes de los elementos que
29 % e s t n en T y el resto 0 ( cuando tenga que escoger a partir de x4 el m x
30 % que no tenga en cuenta los inf )
31 M_X = zeros ( n ) ;
32 for i = 1: n -1
33 M_X ( T ( i ) , T ( i +1) ) = matrizcostes ( T ( i ) , T ( i +1) ) ;
34 M_X ( T ( i +1) , T ( i ) ) = matrizcostes ( T ( i +1) , T ( i ) ) ;
35 end
36 M_X ( T ( n ) , T (1) ) = matrizcostes ( T ( n ) , T (1) ) ;
37 M_X ( T (1) , T ( n ) ) = matrizcostes ( T (1) , T ( n ) ) ;
38
39 % Crear una matriz en donde solo aparezcan los costes de los elementos que
40 % NO e s t n en T ( y sus inversas ) y el resto inf
41 M_Y = matrizcostes ;
42 for i = 1: n -1
43 M_Y ( T ( i ) , T ( i +1) ) = inf ; M_Y ( T ( i +1) , T ( i ) ) = inf ;
44 end
45 M_Y ( T ( n ) , T (1) ) = inf ; M_Y ( T (1) , T ( n ) ) = inf ;
46
47 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
48 % Empiezo el proceso de LK

54
49 back_1 = 1: n ; %Creo una lista con los v r t i c e s existentes para t1
50 while ~ isempty ( back_1 )
51 %% % % % % % % % % % % % % % % % % % % % % % % NIVEL 1 % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
52 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % t1
53 % Escoger t1 ale atoriam ente
54 t1 = back_1 ( randi ( length ( back_1 ) ) ) ; % Escoger t1 aleatorio
55 back_1 = setdiff ( back_1 , t1 ) ; % Eliminar t1 de la lista de opciones a
explorar
56 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % t2
57 % Obtener los valores adyacentes a t1 en V , dos opciones para t2
58 back_2 = find ( M_X ( t1 ,:) ~= 0) ;
59 while ~ isempty ( back_2 )
60 % Obtener los valores adyacentes a t1 en V , dos opciones para t2
61 %Escoger un t2 de manera aleatoria
62 t2 = back_2 ( randi ( length ( back_2 ) ) ) ;
63 back_2 = setdiff ( back_2 , t2 ) ; %guardo la otra o p c i n para el
backtraking
64 X = [ t1 , t2 ];
65 contador_3 = 0;
66 %% % % % % % % % % % % % % % % % % % % % % % % % % % % t3
67 % Buscar y1 que cumpla las condiciones :
68 % 1. No este en T y no sea de los v r t i c e s ya usados en X
69 % 2. Tenga peso m n i m o
70 aux =[1: n ; M_Y ( t2 ,:) ];
71 aux = aux (: , setdiff (1: n , X ) ) ;
72 infs = isinf ( aux (2 ,:) ) ;
73 aux (: , infs ) =[];
74 [~ , ord_t3 ] = sort ( aux (2 ,:) ) ;
75 back_3 = aux (1 , ord_t3 ) ;
76 while ~ isempty ( back_3 ) && contador_3 < 5
77 t3 = back_3 (1) ;
78 contador_3 = contador_3 + 1;
79 %Eliminar la o p c i n ya utilizada
80 back_3 (1) =[];
81 ganancia = M_X ( t1 , t2 ) - M_Y ( t2 , t3 ) ; % G1 = g1
82 if ganancia <1
83 break %vuelve a buscar otro t2 porque no genera ganancia
84 end
85 Y = [ t2 , t3 ];
86 XY = [ t1 , t2 , t3 ]; %ciclo alternado
87 %% % % % % % % % % % % % % % % % % % NIVEL 2 % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
88 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % t4
89 %Opciones v l i d a s para t4
90 opt_t4 = find ( M_X ( XY ( end ) ,:) ~= 0) ;
91 % Filtrar las condiciones que tiene que cumplir t4 para
92 % Escojo t4 de modo que cumpla las siguientes condiciones :
93 % 1. M_X ( t4 , t1 ) == 0
94 % 2. t4 no es un v r t i c e que ya e s t en XY
95 valid_opts = opt_t4 ( M_X ( opt_t4 , XY (1) ) == 0) ;
96 back_4 = valid_opts (~ ismember ( valid_opts , XY ) ) ;
97 while ~ isempty ( back_4 )
98 % Declaro las variables para entrar en caso de que
99 % fuesen modificadas
100 X =[ t1 , t2 ];
101 Y =[ t2 , t3 ];
102 XY =[ t1 , t2 , t3 ];
103 [ T_2 ,X , XY , back_4 , bool ]= encontrar_t4 (X ,Y , XY ,E , M_X , M_Y , back_4 ) ;
104 t4 = XY ( end ) ;
105 % Solo si bool ==1 ( en caso de backtracking )
106 if bool ==1
107 ganancia_est = M_X ( t3 , t4 ) - M_Y ( t4 , t1 ) ; % x1y1x2y2 *
gana ncia_2_e st
108 ganancia_temp = ganancia + ganancia_est ;
109 if ganancia_temp > G_est % T_2 nueva candidata
110 G_est = ganancia_temp ;
111 T_prov = T_2 ;
112 Xopt =[ t1 , t2 ; t3 , t4 ];
113 Yopt =[ t2 , t3 ; t4 , t1 ];
114 k =2;
115 end

55
116 end
117 contador_5 = 0;
118 %% % % % % % % % % % % t5 y t6
119 aux =[1: n ; M_Y ( t4 ,:) ];
120 aux = aux (: , setdiff (1: n , X ) ) ;
121 infs = isinf ( aux (2 ,:) ) ;
122 aux (: , infs ) =[];
123 [~ , ord_t5 ] = sort ( aux (2 ,:) ) ;
124 back_5 = aux (1 , ord_t5 ) ;
125 while ~ isempty ( back_5 ) && contador_5 < 5
126 % Declaro las variables para entrar en caso de que
127 % fuesen modificadas
128 X =[ t1 , t2 ; t3 , t4 ];
129 Y =[ t2 , t3 ];
130 XY =[ t1 , t2 , t3 , t4 ];
131 [ t6 , sol_ahead ,Y , XY , back_5 , contador_5 ]= encontrar_t5 (X ,Y , XY ,
E , M_X , M_Y , contador_5 , back_5 ) ;
132 if isempty ( t6 ) %Si no e n c o n t r t5 que crease t6
133 break ; %busco otro t4 porque p r o b entre todas
134 %opciones de t5 para ese t4 ( puede no
135 %cerrar ciclo
136 end
137 t5 = XY ( end ) ;
138 ganancia = ganancia + ( M_X ( t3 , t4 ) - M_Y ( t4 , t5 ) ) ; % G2 = g1 + g2
139
140 % Ganancia de cerrar con t6
141 ganancia_est = M_X ( t5 , t6 ) - M_Y ( t6 , t1 ) ; %x1y1x2y2x3y3 *
142 ganancia_temp = ganancia + ganancia_est ;
143 if ganancia_temp > G_est % sol_ahead nueva candidata
144 G_est = ganancia_temp ;
145 T_prov = sol_ahead ;
146 k =3;
147 Xopt =[ t1 , t2 ; t3 , t4 ; t5 , t6 ];
148 Yopt =[ t2 , t3 ; t4 , t5 ; t6 , t1 ];
149 end
150 if ganancia <= G_est
151 if G_est >0 % Salgo con intercambio 2 - opt o 3 - opt
152 T_nuevo = T_prov ;
153 return
154 else %pruebo con otro t5
155 continue
156 end
157 end
158 %% % % % % % % % % % % % % % % % % % NIVEL 3 y SUCESIVOS % % % % % % % % % % %
159 seguir = 1; % Mientras se pueda seguir iterando
160 txi = t6 ;
161 while seguir
162 [X ,Y , XY , T_i , txi ] = a um en t ar _c am i no (T , M_X , M_Y ,X ,Y , XY , txi
);
163 s = size (X ,1) ;
164 if ~ isempty ( txi )
165 ganancia = ganancia +( M_X ( X ( end ,1) ,X ( end ,2) ) - M_Y ( Y ( end
,1) ,Y ( end ,2) ) ) ; % G_i = G_ {i -1}+ gi
166
167 ganancia_est = M_X ( Y ( end ,2) , txi ) - M_Y ( txi , t1 ) ; %
x1y1x2y2 ... xiyi *
168 ganancia_temp = ganancia + ganancia_est ;
169 if ganancia_temp > G_est % sol_ahead nueva candidata
170 G_est = ganancia_temp ;
171 T_prov = T_i ;
172 Xopt =[ X ; Y ( end ,2) , txi ];
173 Yopt =[ Y ; txi , t1 ];
174 k = s +1;
175 end
176 if ganancia <= G_est
177 if G_est >0
178 T_nuevo = T_prov ;
179 k;
180 % end
181 seguir =0;

56
182 return
183 else %hago backtracking a t5
184 break
185 end
186 end
187 else % Me q u e d a r a con el ltimo intercambio porque
no se pudo hacer m s
188 if G_est >0
189 T_nuevo = T_prov ;
190 return
191 else %bactraking
192 end
193 seguir =0;
194 end % txi
195 end %seguir
196 end %back_5
197 end %back_4
198 end %back_3
199 end %back_2
200 end %back_1
201 if isempty ( back_1 )
202 Xopt = [];
203 Yopt = [];
204 T_nuevo = T ;
205 k =0;
206 end
207 end

C.3. encontrar t4
1 function [ T_nuevo ,X , XY , back_4 , bool_prov ] = encontrar_t4 (X ,Y , XY ,E , M_X , M_Y , back_4 )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %5 % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que determina el v r t i c e t4 de x2 de modo que pueda cerrar ciclo
4 % en caso de ser necesario
5 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
6 % INPUTS :
7 % X : Matriz kx2 que contiene en las filas las aristas a eliminar
8 % Y : Matriz kx2 que contiene en las filas las aristas a a a d i r
9 % XY : Vector que representa el ciclo alterno de longitud 2 k
10 % E : Matriz nx2 que contiene las aristas existentes en T
11 % M_X : Matriz nxn que contiene los costes de las aristas que forman
12 % parte de la solucion T . El resto de elementos toman valor 0
13 % M_Y : Matriz nxn que contiene los costes de las aristas que NO
14 % forman parte de la solucion T . El resto de elementos toman valor
INF
15 % back_4 : Vector 1 xi que contiene las posibles opciones validas para
16 % t4
17 % OUTPUTS :
18 % T_nuevo : Matriz nx2 que contiene una s o l u c i n valida x1y1x2y2 *
19 % X : Matriz k +1 x2 que contiene en las filas las aristas a eliminar
20 % XY : Vector que representa el ciclo alterno de longitud 2( k +1)
21 % back_4 : Vector 1 x (i -1) que contiene las posibles opciones validas para
22 % t4
23 % bool_prov : Variable booleana que identifica si se encontro variable
24 % valida o no
25 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
26 bool_prov =0;
27 i =1;
28 while bool_prov ==0 && i < length ( back_4 ) +1
29 t4 = back_4 ( i ) ;
30 X_temp =[ X ; XY ( end ) , t4 ];
31 XY_temp =[ XY , t4 ];
32 % Se calcula posible Y en caso de ser necesario cerrar ciclo
33 y2_est = [ t4 , XY (1) ];
34 Y_temp = [ Y ; y2_est ];
35 [ bool_prov , sol ] = es_tour (E , X_temp , Y_temp ) ;
36 i = i +1;
37 end
38 if bool_prov ==1
39 T_nuevo = sol (: ,1) ’;

57
40 else
41 %No encuentro t4 que cierre
42 T_nuevo = [];
43 end
44 X = X_temp ;
45 XY = XY_temp ;
46 back_4 = setdiff ( back_4 , t4 ) ; % Elimino el que ya probe
47 end

C.4. encontrar t5
1 function [ t6 , sol_ahead ,Y , XY , back_5 , contador_5 ]= encontrar_t5 (X ,Y , XY ,E , M_X , M_Y , cont
, back_5 )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %5 % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que determina el v r t i c e t5 de y2 de modo que rompa un x3
4 % verificando la existencia de y3 *
5 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
6 % INPUTS :
7 % X : Matriz kx2 que contiene en las filas las aristas a eliminar
8 % Y : Matriz kx2 que contiene en las filas las aristas a a a d i r
9 % XY : Vector que representa el ciclo alterno de longitud 2 k
10 % E : Matriz nx2 que contiene las aristas existentes en T
11 % M_X : Matriz nxn que contiene los costes de las aristas que forman
12 % parte de la s o l u c i n T . El resto de elementos toman valor 0
13 % M_Y : Matriz nxn que contiene los costes de las aristas que NO
14 % forman parte de la s o l u c i n T . El resto de elementos toman valor
INF
15 % cont : Valor que cuenta el n m e r o de t5 que fueron probados
16 % back_5 : Vector 1 xi que contiene las posibles opciones validas para
17 % t5
18 % OUTPUTS :
19 % t6 : V r t i c e que determina x4
20 % sol_ahead : Matriz nx2 que contiene una s o l u c i n v l i d a x1y1x2y2x3y3 *
21 % Y : Matriz k +1 x2 que contiene en las filas las aristas a a a d i r
22 % XY : Vector que representa el ciclo alterno de longitud 2( k +1)
23 % contador_5 : Valor que cuenta el numero de t5 que fueron
24 % probados ( implementado en una unidad si se encontro t5 valido )
25 % back_5 : Vector 1 x (i -1) que contiene las posibles opciones v l i d a s para
26 % t5
27 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
28 n = size ( M_X ,1) ;
29 % Se calcula t5 ( y2 ) de modo que rompa la arista x3 para entrar en el bucle con
un posible t6
30 %Buscar y2 que rompa un x3
31 t5_encontrado =0;
32 contador_5 = cont ;
33 while ~ isempty ( back_5 ) && t5_encontrado ==0
34 t5 = back_5 (1) ;
35 back_5 (1) =[];
36 Y_pos = [ Y ; XY ( end ) , t5 ];
37 XY_pos = [ XY , t5 ];
38
39 opt_t6 = find ( M_X ( t5 ,:) ~= 0) ;
40 % Escojo t6 de modo que cumpla las siguientes condiciones :
41 % 1. M_X ( t6 , t1 ) == 0
42 % 2. t6 no es un v r t i c e que ya e s t en XY
43 opt_val = opt_t6 ( M_X ( opt_t6 , XY (1) ) == 0) ;
44 back_6 = opt_val (~ ismember ( opt_val , XY_pos ) ) ;
45 if isempty ( back_6 ) % No se encuentran opciones v l i d a s para t6 -> otro t5
46 continue ;
47 else
48 % Seleccionar un t6 aleatorio de las opciones v l i d a s
49 for j = 1: length ( back_6 )
50 t6 = back_6 ( j ) ;
51 % Se actualiza X y XY con la s e l e c c i n
52 X_ahead = [ X ; XY_pos ( end ) , t6 ];
53 XY_ahead = [ XY_pos , t6 ];
54 % Calcular y3_est y verificar el ciclo con es_tour
55 y3_est_ahead = [ t6 , XY (1) ];
56 Y_temp_ahead = [ Y_pos ; y3_est_ahead ];

58
57 [ bool , sol2 ] = es_tour (E , X_ahead , Y_temp_ahead ) ;
58 if bool ==1 %Encontro cierre salgo del for
59 t5_econtrado = 1;
60 break ;
61 end
62 end
63 %Sale del for anterior si ha encontrado un bool =1 o si todos eran bool
=0
64 if bool ==1
65 %Esto significa que e n c o n t r un t6 v l i d o para un t5
66 % Se e n c o n t r un t5 v l i d o
67 contador_5 = cont + 1;
68 % Actualizar Y con el par [ t4 , t5 ] de menor peso
69 Y = Y_pos ;
70 XY = XY_pos ;
71 sol_ahead = sol2 (: ,1) ’;
72 return ;
73 else %para ese t5 no hay t6 que cierre
74 continue ;
75 end
76 end
77 end
78 if contador_5 == cont
79 t6 =[];
80 sol_ahead =[];
81 end
82 end

C.5. aumentar camino

1 function [ Xnew , Ynew , XYnew , T_nuevo , txi ] = au m en ta r_ c am in o (T , M_X , M_Y ,X ,Y , XY , txi )
2 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3 % F u n c i n que devuelve los conjuntos X , Y de aristas para realizar un k - opt

4 %intercambio donde k es el n m e r o de aristas in tercamb iadas
5 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

6 % INPUTS :
7 % T : Vector fila que contiene la s o l u c i n inicial
8 % M_X : Matriz nxn que contiene los costes de las aristas que forman
9 % parte de la s o l u c i n T . El resto de elementos toman valor 0
10 % M_Y : Matriz nxn que contiene los costes de las aristas que NO
11 % forman parte de la s o l u c i n T . El resto de elementos toman
valor INF
12 % X : Matriz kx2 que contiene en las filas las aristas a eliminar
13 % Y : Matriz kx2 que contiene en las filas las aristas a a a d i r
14 % XY : Vector que representa el ciclo alterno de longitud 2 k
15 % txi : V r t i c e que determina X , calculado en el paso anterior
16 % OUTPUTS :
17 % Xnew : Matriz ( k +1) x2 que contiene en las filas las aristas a
eliminar
18 % Ynew : Matriz ( k +1) x2 que contiene en las filas las aristas a
a adir
19 % XYnew : Vector que representa el ciclo alterno resultante de
longitud 2( k +1)
20 % T_nuevo : Array en el que cada celda tiene la s o l u c i n encontrada
21 % con un intercambio k
22 % txi : V r t i c e que determina la siguente arista a introducir
23 % en X
24 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

25 seguir =0;
26 n = length ( T ) ;
27 E = soluc_aristas ( T ) ;
28 %Inicializo las variables de salida
29 Xnew = X ;
30 Ynew = Y ;

59
31 XYnew = XY ;
32
33 %Declaro xi con lo calculado en el paso anterior
34 X = [ X ; XY ( end ) , txi ];
35 XY = [ XY txi ];
36 %% % % % % % % % % % % % % % % % % % % % % CAMBIADO
37 % Buscar un yi que rompa un x_ { i +1}
38 %ind_rest = setdiff (1: n , XY ) ;
39 % Ordenar de menor a mayor ind_rest por los pesos M_Y ( XY ( end ) , tyi )
40 %[~ , ind_ord ] = sort ( M_Y ( XY ( end ) , ind_rest ) ) ;
41 %i n d _ r e s t _ o r d _ t e m p = ind_rest ( ind_ord ) ;
42 %ind_rest_ord = i n d _ r e s t _ o r d _ t e m p (~ isinf ( M_Y ( XY ( end ) , i n d _ r e s t _ o r d _ t e m p ) ) ) ; %Le
quito los infinitos porque esas aristas e s t n en E
43 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
44 aux =[1: n ; M_Y ( XY ( end ) ,:) ];
45 aux = aux (: , setdiff (1: n , XY ) ) ;
46 infs = isinf ( aux (2 ,:) ) ;
47 aux (: , infs ) =[];
48 [~ , ind_ord ] = sort ( aux (2 ,:) ) ;
49 ind_rest_ord = aux (1 , ind_ord ) ;
50
51 for i = 1: length ( ind_rest_ord )
52 tyi = ind_rest_ord ( i ) ;
53 Y_pos = [ Y ; XY ( end ) , tyi ];
54 XY_pos = [ XY , tyi ];
55 % Busco que cumpla la rotura de x_ { i +1}
56 opt_txi = find ( M_X ( tyi ,:) ~= 0) ;
57 % Escojo txi de modo que cumpla las siguientes condiciones :
58 % 1. M_X ( txi , t1 ) == 0
59 % 2. txi no es un v r t i c e que ya e s t en XY
60 valid_opts = opt_txi ( M_X ( opt_txi , XY (1) ) == 0) ;
61 opt_val = valid_opts (~ ismember ( valid_opts , XY_pos ) ) ;
62 if isempty ( opt_val )
63 continue
64 end
65 % Ordenar txi por peso y tomar hasta dos opciones
66 [~ , ind ] = sort ( M_X ( tyi , opt_val ) , ’ descend ’) ; %ordena de mayor a menor
67 back_i = opt_val ( ind ) ; % Reordenamos opt_val s e g n los pesos
68 for j = 1: length ( back_i ) % Iterar sobre las dos opciones y quedarse con la
mejor
69 txi = back_i ( j ) ;
70 X_ahead = [ X ; tyi , txi ];
71 XY_ahead = [ XY_pos , txi ];
72 yi_est_ahead = [ txi , XY (1) ];
73 Y_temp_ahead = [ Y_pos ; yi_est_ahead ];
74 % Verificar el ciclo con es_tour
75 [ bool , sol_ahead ] = es_tour (E , X_ahead , Y_temp_ahead ) ;
76 if bool ==1
77 break
78 end
79 end
80 if bool ==1
81 Xnew = X ;
82 Ynew = [ Y ; XY ( end ) , tyi ];
83 XYnew = [ XY , tyi ];
84 T_nuevo = sol_ahead (: ,1) ’;
85 seguir = 1;
86 return ; % Termina la f u n c i n si encuentra una s o l u c i n v l i d a
87 end
88 end
89 if seguir ==0
90 txi =[];
91 T_nuevo = [];
92 end
93 end

C.6. calcula coste

1 function coste = calcula_coste (T , matrizcostes )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

60
3 %Funcion que calcula el coste del ciclo completo
4 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
5 % INPUTS :
6 % ciclo : Vector fila que contiene los v r t i c e s en el orden en el que se
7 % recorren
8 % matrizcostes : Matriz nxn que contiene los costes de las aristas
9 % OUTPUTS :
10 % coste : Valor que representa el coste del ciclo
11 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
12 %Inicializo la salida
13 coste = 0;
14 n = length ( T ) ;
15 % Calcula el coste total del tour
16 for i = 1: n -1
17 coste = coste + matrizcostes ( T ( i ) , T ( i +1) ) ;
18 end
19 coste = coste + matrizcostes ( T ( n ) ,T (1) ) ;

C.7. es tour
1 function [ bool , sol ] = es_tour (E ,X , Y )
2 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
3 % F u n c i n que comprueba si se crea un tour al eliminar las aristas X y
4 % a a d i r las de Y
5 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
6 % INPUTS :
7 % E : Matriz 2* nx2 que contiene las aristas ( en las dos direcciones posibles )
8 % que formar parte del ciclo
9 % X : Matriz kx2 cuyas filas s e r n las aristas a eliminar del ciclo inicial
10 % Y : Matriz kx2 cuyas filas s e r n las aristas a introducir al ciclo
11 % inicial OUTPUTS : bool : Booleano que representa si se econtro o no un
12 % tour correcto sol : Ciclo ordenado construido con E \ X U Y
13 %% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
14 % Inicializo la salida
15 n = size (E ,1) ;
16
17 % Com probaci ones previas
18 % Comprobar que las aristas de Y ( o sus inversas ) no forman parte del ciclo
19 for k =1: size (Y ,1)
20 idx_Y = ismember (E , Y (k ,:) , ’ rows ’) | ismember (E , flip ( Y (k ,:) ) , ’ rows ’) ;
21 E ( idx_Y , :) = [];
22 end
23 if size (E ,1) ~= n
24 disp ( ’ error : Alguna arista de Y esta en E ’)
25 return
26 end
27
28 % Comprobar que las aristas de X ( o sus inversas ) forman parte del ciclo
29 for k =1: size (X ,1)
30 idx_X = ismember (E , X (k ,:) , ’ rows ’) | ismember (E , flip ( X (k ,:) ) , ’ rows ’) ;
31 E ( idx_X , :) = [];
32 end
33
34 if size (E ,1) ~= n - size (X ,1)
35 disp ( ’ error : No todas las aristas de X estan en E ’)
36 return
37 end
38 % C o m p r o b a c i n de que el intercambio es posible
39 % Esto relmente solo se cumple si es para cerrar el ciclo ( la segunda
40 % condicion
41 if all ( Y (: ,1) == X (: ,2) ) && all ( Y (: ,2) ==[ X (2: end ,1) ; X (1 ,1) ])
42 else
43 disp ( ’ Las X e Y no hacen ciclo ’)
44 return
45 end
46 % C o n s t r u c c i n del ciclo ordenado con las aristas c o r r e s po n d i e n t e s
47 %interca mbiadas
48 Eord ={}; %Array fila donde cada celda contiene una matriz (k ,2) con una cadena
49 % tras haber extraido las aristas de X
50 for i =1: size (E ,1)

61
51 Eord = v e r i f i c a r _ e x c i c l o ( Eord , E (i ,:) ) ;
52 end
53
54 % " Pegado " de las aristas de Y a alguna de las cadenas
55 for j =1: size (Y ,1)
56 pegado = 0;
57 for i = 1: length ( Eord )
58 if Y (j ,2) == Eord { i }(1 , 1) && pegado == 0
59 Eord { i } = [ Y (j ,:) ; Eord { i }];
60 pegado = 1;
61 elseif Y (j ,2) == Eord { i }( end , 2) && pegado == 0
62 Eord { i } = [ Eord { i }; flip ( Y (j ,:) ) ];
63 pegado = 1;
64 elseif Y (j ,1) == Eord { i }( end , 2) && pegado == 0
65 Eord { i } = [ Eord { i }; Y (j ,:) ];
66 pegado = 1;
67 elseif Y (j ,1) == Eord { i }(1 , 1) && pegado == 0
68 Eord { i } = [ flip ( Y (j ,:) ) ; Eord { i }];
69 pegado = 1;
70 end
71 end
72 end
73
74
75 %Pegar las posibles cadenas que pueda haber en Eord entre si . Habiendo 4
76 % posibilidades de pegado
77 i mp os ib l e_ pe ga r =0;
78 while length ( Eord ) >1 && i mp os i bl e_ pe g ar ==0
79 for k =2: length ( Eord )
80 encontrado =0;
81 if Eord {1}( end ,2) == Eord { k }(1 ,1) % ltima con primera
82 encontrado =1;
83 elseif Eord {1}(1 ,1) == Eord { k }(1 ,1) % Primera con primera
84 Eord {1}= flip ( Eord {1}) ;
85 Eord {1}= flip ( Eord {1} ,2) ;
86 encontrado =1;
87 elseif Eord {1}(1 ,1) == Eord { k }( end ,2) % Primera con ltima
88 Eord {1}= flip ( Eord {1}) ; Eord {1}= flip ( Eord {1} ,2) ;
89 Eord { k }= flip ( Eord { k }) ; Eord { k }= flip ( Eord { k } ,2) ;
90 encontrado =1;
91 elseif Eord {1}( size ( Eord {1} ,1) ,2) == Eord { k }( size ( Eord { k } ,1) ,2) % ltima con
ltima
92 Eord { k }= flip ( Eord { k }) ;
93 Eord { k }= flip ( Eord { k } ,2) ;
94 encontrado =1;
95 end
96
97 if encontrado ==1
98 Eord {1}=[ Eord {1}; Eord { k }]; %concatena las dos cadenas
99 Eord ( k ) =[];
100 break ;
101 else
102 i mp os ib l e_ pe ga r =1;
103 end
104 end
105 end
106 % Definir la salida si fue o no posible pegar y por tanto generando una
107 % unica celda en Eord
108 if length ( Eord ) ==1
109 bool =1;
110 sol = Eord {1};
111 else
112 bool =0;
113 sol = Eord ;
114 end
115
116 end

C.8. soluc aristas

62
1 function E = soluc_aristas ( T )
2 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3 % F u n c i n que crea a partir de una s o l u c i n representada por los v r t i c e s

4 % en orden la s o l u c i n con la secuencia de aristas en la s o l u c i n
5 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

6 % INPUTS :
7 % T : Vector fila que contiene la s o l u c i n con los v r t i c e s en orden
8 % s e g n el recorrido
9 % OUTPUTS :
10 % E : Matriz 2 nx2 que contiene las aristas de la s o l u c i n y sus
11 % inversas
12 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

13 % Inicializamos la salida
14 n = length ( T ) ;
15 E = zeros (n , 2) ; %MONI
16 % Aristas del vector
17 for i = 1: n -1
18 E (i , 1) = T ( i ) ;
19 E (i , 2) = T ( i +1) ;
20 end
21 E (n ,1) = T ( n ) ; E (n ,2) = T (1) ;
22 end

63
Apéndice B

Tablas de datos

A. Tabla de coordenadas de las 38 ciudades

Vértice Ubicación Latitud Longitud

1 Ferrol 43,48268 -8,2238
2 A Coruña 43,36762 -8,42287
3 Nois 43,61552 -7,31274
4 Villalba 43,29748 -7,68077
5 Chavı́n 43,61076 -7,5798
6 Castrelo de Miño 42,29728 -8,05742
7 Santiago de Compostela 42,87686 -8,54417
8 Lugo 43,00973 -7,55675
9 Samil 42,20962 -8,77215
10 Sanxenxo 42,4033 -8,81135
11 Sarria 42,78083 -7,41407
12 Ribadeo 43,53367 -7,04034
13 Comarca de Valdeorras 42,37886 -6,96386
14 Viveiro 43,66426 -7,59453
15 Mondoñedo 43,42857 -7,36384
16 Verı́n 41,94025 -7,43495
17 Finisterre 42,9078 -9,26503
18 Lage 43,22121 -8,94043
19 Fonsagrada 43,12471 -7,06885
20 Lalı́n 42,66115 -8,11102
21 Lobios 41,90113 -8,08335
22 Monforte de Lemos 42,51902 -7,5158
23 Puebla de Trives 42,33934 -7,2538
24 A Gudiña 42,06057 -7,1393
25 Ordes 43,07637 -8,40806
26 Curtis 43,12533 -8,14599
27 Piornedo 42,85648 -6,87581
28 Baralla 42,8939 -7,25042
29 Tui 42,04915 -8,6466
30 Xinzo de Limia 42,06311 -7,72594
31 Meira 43,21342 -7,29428
32 Melide 42,58948 -7,75338
33 Cerdedo 42,53219 -8,38949
34 Visuña 42,60963 -7,06917
35 Mondariz 42,22577 -8,46878

64
 36 Baamonde 43,17667 -7,75658
37 Carballo 43,21282 -8,69152
38 Taboada 42,71558 -7,76208

B. Tabla con las distancias entre las aristas existentes de las

18 ciudades

Vértice 1 Ciudad 1 Vértice 2 Ciudad 2 Distancia en km

1 Ferrol 7 Viveiro 76
13 Betanzos 1 Ferrol 35
1 Ferrol 11 Villalba 58
2 A Coruña 8 Finisterre 101
13 Betanzos 2 A Coruña 22
4 Santiago de Compostela 8 Finisterre 78
4 Santiago de Compostela 10 Tui 102
4 Santiago de Compostela 13 Betanzos 63
4 Santiago de Compostela 14 Lalı́n 50
4 Santiago de Compostela 3 Castrelo de Miño 103
3 Castrelo de Miño 10 Tui 70
3 Castrelo de Miño 14 Lalı́n 55
3 Castrelo de Miño 12 Monforte de Lemos 68
3 Castrelo de Miño 9 Puebla de Trives 93
9 Puebla de Trives 12 Monforte de Lemos 49
12 Monforte de Lemos 14 Lalı́n 65
12 Monforte de Lemos 5 Lugo 62
5 Lugo 11 Villalba 38
5 Lugo 6 Ribadeo 90
5 Lugo 14 Lalı́n 73
5 Lugo 15 Guitiriz 42
15 Guitiriz 11 Villalba 30
15 Guitiriz 13 Betanzos 33
6 Ribadeo 7 Viveiro 59
6 Ribadeo 11 Villalba 72
7 Viveiro 11 Villalba 52
16 Arzua 5 Lugo 66
16 Arzua 4 Santiago de Compostela 37
16 Arzua 14 Lalı́n 39
17 Curtis 15 Guitiriz 29
17 Curtis 13 Betanzos 22
17 Curtis 4 Santiago de Compostela 49
17 Curtis 16 Arzua 27
18 Cerceda 13 Betanzos 36
18 Cerceda 2 A Coruña 26
18 Cerceda 4 Santiago de Compostela 50
18 Cerceda 8 Finisterre 93

65