El problema paramétrico del emparejamiento en grafos y problema de emparejamiento con dos objetivos

TRABAJO FIN DE MÁSTER El problema paramétrico del emparejamiento en grafos y problema de emparejamiento con dos objetivos Presentado por: Rafael González López Supervisado por: DR. JUSTO PUERTO ALBONDOZ FACULTAD DE MATEMÁTICAS Departamento de Estadística e Investigación Operativa Sevilla, Junio 2019 Índice general Abstract 5 Introducción 7 1. El problema del emparejamiento 9 1.1. El problema del b-emparejamiento . . . . . . . . . . . . . . . . . . . 9 1.1.1. El problema del emparejamiento biobjetivo . . . . . . . . . . 12 1.2. El algoritmo de Grötschel-Holland . . . . . . . . . . . . . . . . . . . 13 1.3. El problema de separación del matching . . . . . . . . . . . . . . . . 17 1.3.1. Cortes y árbol de Gomory-Hu . . . . . . . . . . . . . . . . . 17 1.3.2. El problema de separación . . . . . . . . . . . . . . . . . . . 18 1.4. Aplicaciones y problemas relacionados . . . . . . . . . . . . . . . . . 19 1.4.1. El problema de la asignación . . . . . . . . . . . . . . . . . . 19 1.4.2. Aproximación al problema del viajante . . . . . . . . . . . . 20 1.4.3. Cobertura de aristas . . . . . . . . . . . . . . . . . . . . . . . 21 1.4.4. Coloración de aristas . . . . . . . . . . . . . . . . . . . . . . 22 1.4.5. T-Enlaces y el problema del cartero chino . . . . . . . . . . . 23 1.4.6. Localización de objetos en el espacio . . . . . . . . . . . . . 24 1.4.7. Movimiento de objetos . . . . . . . . . . . . . . . . . . . . . 25 1.4.8. Orden óptimo del inventario . . . . . . . . . . . . . . . . . . 25 1.4.9. Envolvente convexa del b-matching .............. 26 2. Algoritmo SAP 27 2.1. Propiedades teóricas . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.1.1. El problema del matching y las familias de contracción . . . . 29 2.2. Fundamentos del algoritmo . . . . . . . . . . . . . . . . . . . . . . . 30 2.2.1. Transformaciones admisibles . . . . . . . . . . . . . . . . . . 31 2.3. Algoritmo del camino de aumento más corto . . . . . . . . . . . . . . 33 2.4. Método de reoptimización . . . . . . . . . . . . . . . . . . . . . . . 36 2.4.1. Problema biobjetivo . . . . . . . . . . . . . . . . . . . . . . 38 2.5. Manteniendo optimalidad . . . . . . . . . . . . . . . . . . . . . . . . 39 2.5.1. Problema biobjetivo . . . . . . . . . . . . . . . . . . . . . . 39 3 4 Índice general 3. Análisis computacional 41 3.1. Implementación en Python . . . . . . . . . . . . . . . . . . . . . . . 41 3.1.1. Método de Padberg-Rao modificado . . . . . . . . . . . . . . 43 3.1.2. Algoritmo de Grötschel-Holland . . . . . . . . . . . . . . . . 45 3.1.3. Método del simplex sobre el poliedro del matching . . . . . . 49 3.1.4. Resolución como problema de Programación Entera . . . . . 50 3.1.5. Grafos de Erdös-Rényi . . . . . . . . . . . . . . . . . . . . . 52 3.1.6. Reoptimización......................... 54 3.1.7. Dibujando matchings . . . . . . . . . . . . . . . . . . . . . . 56 3.1.8. SCIP y soluciones no soportadas . . . . . . . . . . . . . . . . 58 3.2. Resultados computacionales . . . . . . . . . . . . . . . . . . . . . . 60 3.3. Reoptimización............................. 63 4. Conclusiones 67 Abstract The minimum (maximum) matching problem is a fundamental problem in combinatorial optimization. The special structure of matching found by Edmonds allows to use different approximations. The idea of this work is to present a general insight of the matching problem, different sorts of algorithms for solving its original and the bi-objective formulations and a computational approach. In section 1 we describe the problem, its applications and the Grötschel-Holland algorithm (for solving using linear programming). In section 2 we introduce the theoretical background necessary for understanding the SAP algorithm and the sensitive analysis. In section 3 we implement some methods in Python 3.7 with Gurobi to study their computational properties and sensitive analysis. In section 4 we present our conclusions. 5 Introducción El problema del emparejamiento es una de las cuestiones clásicas de la Programación Matemática. En este trabajo desarrollamos un extenso análisis sobre el trasfondo teórico de ciertos métodos para resolverlo, tanto en su formulación original, como en la formulación biobjetivo. Tras una presentación inicial del problema, presentamos una primera aproximación a la resolución del problema mediante la Programación Lineal conocida como el algoritmo de Grötschel-Holland. Seguidamente, veremos numerosas aplicaciones que motivan el interés del problema. En el segundo capítulo nos centramos en los conceptos y la estructura intrínseca del problema del emparejamiento, la cuál nos permite finalmente detallar el algoritmo del SAP. Además, presentamos el método de reoptimización propuesto por Derigs para este algoritmo. El tercer epígrame lo dedicamos a realizar un análisis computacional. Implementamos en Python 3.7 distintos métodos de resolución del problema del emparejamiento con el fin de reportar una comparativa entre ellos. Finalmente, llevamos a cabo un análisis de sensibilidad del problema original y la versión biobjetivo. Finalmente, en el cuarto capítulo presentamos las conclusiones. 7 Capítulo 1 El problema del emparejamiento El problema del emparejamiento es uno de los problemas clásicos de la optimización combinatoria. Su solución debido a Edmonds representa uno de los momentos célebres en la historia de la optimización por ser el primer problema de Programación Entera, descontando aquellos que se resuelven usando la relajación lineal, en ser resuelto en tiempo polinomial. El caso bipartito puede verse como un problema de asginación y fue resuelto en tiempo polinómico por Harowld W. Kuhn con su famoso algoritmo húngaro en 1955 [2], predecesor de los algoritmos dual-primal. Sin embargo, hubo que esperar hasta 1961 para que Jack Edmonds descubriera la solución del caso general. En su trabajo original (publicado en 1965 [1]) no solamente caracterizó linealmente el poliedro del matching -abriendo el camino para métodos basados en el simplexsino que utilizando la estructura especial del problema construyó un algoritmo iterativo denominado algoritmo de blossom que resuelve el problema polinomialmente. 1.1. El problema del b-emparejamiento Dado que nuestro objetivo es el estudio paramétrico del problema del emparejamiento y la versión multiobjetivo del mismo, comenzamos este trabajo exponiendo en qué consiste el problema del ejemparejamiento o problema del matching. Para ello, comenzamos definiendo algunos conceptos básicos en el marco de la teoría de grafos. Definición 1.1.1. Sea G= (V, E)un grafo y sea S⊂V, definimos δ(S)como el conjunto de aristas con un único extremo en S. En el caso de un conjunto unitario {i}, denotamos δ(i) := δ({i}). Usualmente se denomina grado del vértice ial cardinal de δ(i). Definimos además γ(S)como el conjunto de aristas que tienen ambos extremos en S. En adelante consideramos siempre G= (V, E)un grafo no dirigido. El problema del matching consiste en encontrar un subconjunto M⊂Econ la propiedad de que en 9 16 1.2. El algoritmo de Grötschel-Holland necesitamos determinar planos de corte para obtener una nueva solución x∗, es decir, vamos al Paso 7. Paso 7. Si la solución x∗no es entera, entonces no puede inducir un matching en nuestro grafo. Necesitamos encontrar planos de corte y, naturalmente, deseamos encontrarlos en el menor tiempo posible. Comenzamos utilizando búsqueda en profundidad, comprobamos si Gx∗tiene componentes conexas con un número impar de nodos. Si lo hay, sean (Vi, Ei) estas componentes de Gx∗, estas dan lugar a unas desigualdades X (u,v)∈Ei xuv ≤1 2(|Vi| − 1) que son violadas por x∗. En este caso no necesitamos buscar más planos de corte, sino directamente utilizamos estos y continuamos al Paso 8. Si nuestra solución no viola las desigualdades anteriores, consideramos el grafo G2 x∗obtenido al eliminar las aristas de Ex∗cuyos costes -recordemos, los kuv definidos anteriormentesean menores que un cierto ε. Aplicamos ahora el procedimiento inicial del Paso 7 al grafo G2 x∗. Incluso si encontramos componentes conexas de cardinalidad impar, estas no tendrían por qué inducir desigualdades que se violasen necesariamente, por lo que tenemos que comprobarlo. Si es así, hemos encontrado planos de corte y pasamos al Paso 8. Si después de utilizar esta segunda heurística tampoco hemos encontrado planos de corte, consideremos V1como el conjunto de vértices de Gx∗incidentes con alguna arista cuyo peso kuv = 1. En tal caso, eliminamos de Gx∗el conjun- to de vértices V1así como todas las aristas que les sean incidentes. Denotamos por G0 x∗= (V0, E0 x∗)al grafo resultante de esta operación. Padberg y Rao presentaron un método para encontrar un plano de corte adecuado [13] basado en determinar el árbol de Gomory-Hu de G0 x∗y encontrar un eje de este árbol de peso mínimo entre aquellos que al ser eliminados dividen el árbol en dos componentes de cardinalidad impar. Notemos que, en particular, G0 x∗tiene un número par de nodos. Una vez encontrado el plano de corte, continuamos con el Paso 8. Notemos que el algoritmo de Padberg y Rao que utilizan los autores para encontrar el plano de corte, a pesar de ser polinomial, es O(n4). Para grafos de gran tamaño, este procedimiento también puede resultar prohibitivo, razón por la cuál se intentan encontrar planos de corte con heurísticas más sencillas. Entre ellas se recomienda el algoritmo de Gusfield. Paso 8. Añadimos los planos de corte encontrados en el paso anterior a nuestro problema y volvemos al Paso 4. Capítulo 1. El problema del emparejamiento 17 Paso 9. Hemos obtenido una solución entera que induce un matching perfecto e (V, E0), pero naturalmente no tiene por qué ser un matching óptimo en Kn. Para comprobarlo, calculamos el coste reducido a cada eje (u, v)∈E\E0. Sea V AR el conjunto de variables que podrían mejorar la solución óptima. Si todos los costes tienen el signo adecuado, entonces V AR =∅y nuestra solución actual es de hecho óptima en Kn, por lo que hemos terminado. Si V AR 6=∅, entonces establecemos E0=E0∪V AR y añadimos las correspondientes variables. Volvemos al Paso 4. Los detalles nuestra implementación pueden verse comentados en el código presente en la tercera parte de este trabajo. 1.3. El problema de separación del matching Para aprovechar el potencial del software expresamente diseñado para resolver problemas de programación lineal, nosotros implementaremos el algoritmo de Padberg- Rao, a pesar de que para el análisis de sensibilidad vendrá dado por el algoritmo del SAP. Esto está justificado, pues consideramos más interesante intentar utilizar la estructura combinatoria propia del problema para dicho análisis, frente a la aproximación clásica cuyos resultados son conocidos y difícilmente interpretables. Para entender mejor cómo funciona el algoritmo anterior y en qué se basa el procedimiento de Padberg-Rao, así como qué relación guarda con el problema del bemparejamiento, discutimos una serie de conceptos y problemas en esta sección. 1.3.1. Cortes y árbol de Gomory-Hu Definición 1.3.1. Sea G= (V, E)un grafo. Decimos que C= (S, T)es un corte de Gsi V=S∪TyS∩T=∅. Definición 1.3.2. Sea G= (V, E)un grafo y sea Cun corte de G. Se define el conjunto de corte de Ccomo {(u, v)∈E|u∈S, v ∈T} Naturalmente, podemos plantearnos los problemas de encontrar el corte cuyo conjunto de corte tenga mínima o máxima cardinalidad o, en caso de que las aristas tengan peso, los conjuntos de cortes de mínimo o máximo peso. Definición 1.3.3. Sea G= (V, E)un grafo y sean s, t ∈V. Decimos que un corte C es un s−tcorte si sytestán en distintos conjuntos de la partición. Además, notemos que si tenemos el conjunto corte se corresponde biunivocamente con el corte asociado. Una propiedad de los conjuntos de corte de un grafo es que tienen una estructura algebraica. 18 1.3. El problema de separación del matching Proposición 1.3.1. Sea G= (V, E)un grafo y sea C(G)el conjunto de los conjuntos de corte de G. Entonces con la operación ∆diferencia simétrica y el cuerpo Z/Z2 forma un espacio vectorial, donde 0·H=∅y1·H=H∀H∈C(G). El concepto de corte es importante en Teoría de Grafos, pues no solamente son centrales en el algoritmo que nos ocupa, sino que además están estrechamente conectados con los denominados problemas de flujo. Gomory y Hu introdujeron en 1961 un algoritmo polinomial para encontrar un grafo de árbol que representa una base del espacio vectorial anterior. Esta representación se conoce como el árbol de Gomory-Hu, y es especial por verificar la siguiente propiedad. Proposición 1.3.2. Sea G= (V, E)un grafo cuyas aristas tienen coste c, sean s, t ∈V y sea Tel árbol de Gomory-Hu asociado al grafo y costes anteriores. Sea Pel único camino dentro de Tque tiene por extremos syt. Sea ela arista de menor peso dentro de P. Entonces la partición de los vértices de Tobtenida al eliminar eproporciona precisamente el corte de coste mínimo tanto en Tcomo en G. 1.3.2. El problema de separación Presentamos ahora el problema de la separación y qué relación tiene con el problema del matching y el procedimiento de Padberg-Rao en el algoritmo de Grötschel- Holland. Sea P⊂Rnun poliedro. El problema de la separación asociado al poliedro P consiste en, dado un vector y∈Rndecidir si y∈Po bien encontrar π∈Rn,π0∈R tales que πx ≤x0∀x∈Ppero πy > π0. Este problem es interesante, pues en [15] podemos encontrar una demostración del siguiente resultado. Teorema 1.3.1. Dado un poliedro P⊂Rny un vector c∈Rn, el problema de optimización m´ın{c0x|x∈P}puede resolverse en tiempo polinómico si y solo si el problema de la separación para el poliedro Ppuede ser resulto en tiempo polinómico. En [13] Padberg y Rao presentan un algoritmo que resuelve el problema de la separación del poliedro del b-matching en tiempo polinómico, pero para ello utilizan una caracterización del poliedro del matching que no hemos visto todavía, mediante las llamadas desigualdades de blossom. Esta caracterización fue probada por Edmonds en [1] y pasamos a continuación a presentarla. Teorema 1.3.2. Sea G= (V, E)un grafo. Entonces, la envolvente convexa del politopo del b-matching viene dada por X (u,v)∈δ(u) xuv ≤1,∀u∈V xuv ≥0 X (u,v)∈δ(W)\F xuv −X (u,v)∈F xuv ≥1− |F| Capítulo 1. El problema del emparejamiento 19 donde W⊂V,F⊂δ(W)y|F|+Pv∈Wbves impar. Además, si nos restringimos a los emparejamientos perfectos, entonces basta considerar únicamente la igualdad en el primer bloque de restricciones. En el caso del problema del emparejamiento que estamos tratando bv= 1 y simplemente |W|+|F|tiene que ser impar. En nuestra implementación del algoritmo de Grötschel-Holland nosotros para el análisis computacional no vamos a utilizar directamente el algoritmo de Padberg-Rao, sino un modificación del mismo que presentan Letchford, Reinelt y Theis en [14], junto con una demostración de la corrección del mismo. Inicialización: Un grafo Gy un vector x Paso 1: Computamos el árbol de Gomory-Hu de Gpara los pesos wuv = m´ın{xuv,1−xuv} Paso 2: Para cada uno de los |V| − 1cortes δ(W)realizamos: Encontramos el mejor conjunto Ftal que F⊂δ(W)con |V|+|W|impar y minimizando el valor X (u,v)∈δ(W)\F xuv −X (u,v)∈F xuv +|F| Si se verifica X (u,v)∈δ(W)\F xuv −X (u,v)∈F xuv +|F|<1 Entonces hemos encontrado la desigualdad deseada. Si tras terminar el bucle no hemos encontrado ninguno, entonces xestá en el poliedro del matching. 1.4. Aplicaciones y problemas relacionados En esta sección presentamos algunas aplicaciones, como la resolución del problema de la asignación o un algoritmo de aproximación del problema del viajante, y otros problemas interesantes, como el problema de la cobertura por aristas, el de la coloración de aristas o los T-enlaces, que guardan cierta relación con el problema del emparejamiento. 1.4.1. El problema de la asignación Para ilustrar las distintas aplicaciones que tiene el problema del matching, comenzamos con un problema clásico de la Investigación Operativa. El problema de asignación o assignment problem consiste en encontrar la forma de asignar recursos a un 20 1.4. Aplicaciones y problemas relacionados cierto conjunto de tareas determinadas con un cierto criterio para evaluar las asignaciones. Se supone que cada recurso se destina a una sola tarea, y que cada tarea es ejecutada por uno solo de los recursos. En muchos contextos deseamos hacer asignaciones, por ejemplo, entre objetos y personas. Por ejemplo, podemos desear asignar tareas, estancias, maquinas, procesos, asientos, horarios, etc. Algunas posibilidades: 1. Una empresa desea contratar nempleados para npuestos de trabajo. Tras realizar unos test de aptitud la empresa asigna a cada persona iy tarea jun coeficiente de competencia cij. El objetivo del problema es encontrar la asignación que maximiza la competencia total. 2. En las fuerzas armadas muchos soldados están cualificados para realizar cierto tipo de tareas. A las fuerzas armadas les gustaría asignar al personal minimizando el coste de movimiento. Las reglas especifican la necesidad de que a ciertos destinos solo pueden ir personal con cierto tipo de cualificación. Para cada posible asignación tenemos asociado un coste de movimiento asociado a enviar a una persona, su familia y sus pertenencias a una nueva residencia. En este caso deseamos minimizar este coste. 3. El dueño de un hostal quiere asignar parejas para ser compañeros de habitación. Dos clientes se pueden emparejar si comparten nacionalidad, religión, trasfondo cultural y hobbys. El problema consiste en encontrar el emparejamiento de máxima cardinalidad. 1.4.2. Aproximación al problema del viajante El problema del viajante o travelling salesman problem (TSP) en un grafo no dirigido consiste en encontrar el ciclo hamiltoniano de peso mínimo. Computacionalmente hablando, se sabe que es NP-duro. Este problema aparece para responder a la pregunta: ¿Dada una lista de ciudades y las distancias entre cada pareja de ciudades, cuál es la ruta más corta para visitar todas las ciudades y volver a la ciudad de origen? Este problema puede formularse como problema de programación entera, como puede encontrarse en [6]. En el algoritmo se utiliza un conocido teorema de teoría de grafos, denominado lema del apretón de manos. Aunque no puede utilizarse el problema del matching para resolver este problema, Nicos Christofides encontró un algoritmo [7] para obtener, bajo ciertas condiciones, una aproximación no mayor que 1,5el valor óptimo del TSP. Para este resultado es necesario imponer que las distancias estén en el marco de un espacio métrico, es decir, han de ser simétricas y verificar la desigualdad Capítulo 1. El problema del emparejamiento 21 triangular. El algoritmo consiste en: Paso 1: Obtener el árbol recubrido mínimo Tde G. Paso 2: Consideremos Oel conjunto de vértices con grado impar Por el Lema del apretón de manos, Otiene cardinalidad par. Hallar el MCPM sobre el grafo inducido por O. Paso 3: Combinar los ejes del árbol y del matching obteniendo un multigrafo euleriano. Encontrar un circuito euleriano. Paso 4: Obtener un circuito hamiltoniano descartando nodos ya visitados 1.4.3. Cobertura de aristas Sea G= (V, E)un grafo. Una cobertura por aristas del grafo Ges un conjunto Cde aristas de manera que todo vértice de Ges incidente con, al menos, un arista de C. Podemos definir además el número de cobertura de aristas de un grafo Gcomo el tamaño de una cobertura de aristas de tamaño mínimo. Se verifica el siguiente resultado. Proposición 1.4.1. Sea Mun matching de cardinalidad máxima en Gy sea Cuna cobertura de aristas de cardinalidad mínima en el grafo G= (V, E). Entonces |M|+ |C|=|V|. Demostración. La demostración de este resultado es sencilla. Sea Uel conjunto de nodos con grado 0relativo a M, entonces |U|=|V| − 2|M|. Dado que añadiendo las aristas de UaMobtendríamos una cobertura por aristas de G, necesariamente |C|≤|M|+|U|=|M|+|V| − 2|M|=|V|−|M| Por otro lado, consideremos el grafo inducido por C, y sea M0un matching de cardinalidad máxima en este subgrafo. Sea U0el conjunto de nodos con grado 0relativo a M0en este subgrafo, entonces |C|=|M0|+|U0|=|V|−|M0| y, por tanto, |M| ≥ |M0|=|V|−|C| de donde se deduce el resultado. Además de la relación anterior, tenemos un teorema análogo al Teorema 1.1.1. para la cobertura de aristas. 22 1.4. Aplicaciones y problemas relacionados Teorema 1.4.1. La envolvente convexa de las coberturas en un grafo G= (V, E)está dada por X (u,v)∈δ(u) xuv ≥1,∀u∈V xuv ≥0 X (u,v)∈γ(S)∪δ(S) xuv ≥1 2(|S|+ 1) ∀S∈B 1.4.4. Coloración de aristas Consideremos ahora el problema de la coloración de aristas: dado un grafo G= (V, E), colorear las aristas de G, con el mínimo número de colores posibles, de manera que aristas que sean incidentes en un mismo vértice tengan colores diferentes. El problema de la coloración por aristas está relacionado con el problema del matching dado que una coloración es factible si y solo si cada conjunto de ejes de un mismo color forma un matching. Por tanto, podemos formular el problema de la cobertura de aristas como el problema de cubrir las aristas de Econ el mínimo número posible de matching maximales. Consideremos χ(G)como el número de colores mínimo con el que se puede colorear G.χ(G)se conoce usualmente como número cromático. Sea además ∆(G) = m´axv∈V|δ(v)|. Veamos algunos resultados relacionados con estos conceptos. Proposición 1.4.2. Si Ges un grafo bipartito, entonces χ(G) = ∆(G). Más generalmente, Vinzing probó el siguiente teorema. Teorema 1.4.2. En general, para cualquier grafo G,χ(G)es ∆(G)o∆(G)+1. La demostración del mismo presenta un algoritmo -bastante rápidopara colorear Gcon ∆(G)+1colores. Esperaríamos que el teorema anterior nos permitiese decidir con facilidad si χ(G)es o no ∆(G), pero en realidad este problema es NP-completo. Es más, determinar si χ(G)es relativamente pequeño puede ser realmente costoso, como muestra el siguiente resultado. Teorema 1.4.3. El problema de decidir si χ(G)≤3es NP-completo. Ahora bien, sea Ala matriz cuyas filas corresponden a los matching maximales de G, entonces podemos formular χ(G) = m´ın y m X i=1 yi s.a. yA·i≥1∀i= 1, . . . , m y∈ {0,1}m Capítulo 1. El problema del emparejamiento 23 La relajación lineal de este problema, así como el dual de la misma, están dados por χLP (G) = m´ın y m X i=1 yi s.a. yA·i≥1∀i= 1, . . . , m yi≥0∀i= 1, . . . , m ∆LP (G) = m´ax x n X i=1 xi s.a. Ax·i≤1∀i= 1, . . . , n xi≥0∀i= 1, . . . , n donde |E|=n. El problema dual puede ser resuelto en tiempo polinomial, pues x∗ (0≤x∗ i≤1) es una solución factible del mismo si y solo si un matching de coste máximo en Gtomando como pesos x∗tiene valor no mayor que 1. 1.4.5. T-Enlaces y el problema del cartero chino Continuamos presentando la relación entre los emparejamientos y los T-enlaces, el problema de emparejamiento y el problema del cartero chino. Definición 1.4.1. Dado un grafo G= (V, E), consideremos T⊂Vcon cardinalidad par y E0⊂E. Decimos que E0es un T-enlace si en G0= (V, E0)si v∈Vtiene grado impar si y solo si v∈T. Proposición 1.4.3. Los T-enlaces de cardinalidad mínima son árboles. Proposición 1.4.4. El problema de encontrar un T-enlace de peso mínimo puede resolverse en tiempo polinomial. Demostración. Encontrar un T-enlace de coste mínimo. Para verlo, veamos que podemos reducir este problema al de encontrar un matching perfecto. Sea G= (V, E)y T, reemplazamos cada v∈Vpor un grafo completo Cvcon δ(v) + αvvértices, donde αv=     0Si v∈Ty|δ(v)|es impar 0Si v /∈Ty|δ(v)|es par 1En otro caso Entonces, para cada (u, v)∈E, unimos un nodo de Cva otro de Cide manera que dos de estos nuevos ejes no sean incidentes en un mismo vértice. Sea G0= (V0, E∪E0) este nuevo grafo, donde E0son las aristas de los grafos completos. Tenemos entonces 24 1.4. Aplicaciones y problemas relacionados Proposición 1.4.5. Si M⊂E∪E0es un matching perfecto en G0, entonces M∩E es un T-enlace en G. Recíprocamente, si E∗es un T-enlace en G, entonces existe K⊂E0tal que E∗∪Kes un matching perfecto en G0. Veamos ahora la importancia de encontrar estos conjuntos de cara a la resolución del problema conocido como problema del cartero chino. Este problema trata de encontrar un el circuito de menor coste que visite todas las aristas. Otro nombre para este problema es el problema de inspección de rutas. Si el grafo es euleriano, necesariamente la solución ha de ser cualquier circuito euleriano. En caso de que nuestro grafo no sea euleriano, podemos utilizar el siguiente algoritmo polinomial basado precisamente en el T-enlace. Paso 1: Si el grafo es euleriano, entonces tomamos un circuito euleriano como solución. En otro caso, vamos al Paso 2. Paso 2: Obtener el T-enlace de coste mínimo, donde consideramos Tcomo el conjunto de los vértices de grado impar. Como hemos visto, este puede obtenerse usando un algoritmo que resuelva el problema del emparejamiento. Paso 3: Consideramos el multigrafo obtenido a partir de duplicar las aristas de los vértices de T. Se demuestra que este nuevo grafo es euleriano. Obtenemos un circuito euleriano. Paso 4: El circuito euleriano de este grafo induce una solución del problema. 1.4.6. Localización de objetos en el espacio Para identificar un objeto en un espacio tridimensional podemos utilizar dos sensores infrarrojos colocados en lugares distintos. Cada sensor nos da un ángulo de visión del objeto y, por tanto, una línea en la cuál debe encontrarse. La intersección de ambas lineas, suponiendo que los dos sensores y el objeto no sean colineales, nos de manera única la posición del objeto. Consideremos la situación de determinar la localización de pobjetos utilizando dos sensores. El primer sensor nos determina L1, . . . , Lplíneas para los pobjetos mientras que el segundo nos determina L0 1, . . . , L0 p. Para identificar los objetos, utilizando que dos líneas identifican un objeto si intersecan, tenemos que emparejar las líneas del primer conjunto con las del segundo. Esta aproximación tiene dos inconvenientes. Por una lado, una línea de un sensor podría intersecar con más de una lineal del otro, por lo que el matching no sería único. Por otro lado, los errores de medición podrían provocar que dos líneas que identifican un objeto no intersecaran. Para lidiar con estos problemas podemos utilizar la siguiente aproximación. Capítulo 1. El problema del emparejamiento 25 Planteemos un problema de asignación de las lineas de un conjunto a las del otro. Definimos el coste cij de la arista (i, j)como el mínimo de las distancias euclídeas entre LiyL0 j. Estudios con simulaciones han descubierto que esta aproximación identifica correctamente los objetos en la gran mayoría de casos. 1.4.7. Movimiento de objetos En muchos contextos diferentes nos interesa intentar aproximar la velocidad y dirección de movimiento de una serie de pobjetos que se mueven en el espacio (por ejemplo: aeronaves, misiles). Para la localización de los objetos podríamos emplear el método visto anteriormente. Una posibilidad para estimar la dirección y velocidad del objeto podría ser calcular la localización de dos objetos en dos instantes distintos de tiempo y emparejar los primeros puntos con los segundos. Si este emparejamiento es correcto, es claro que podemos hacer una estimación de los parámetros. Sean (xi, yi, zi)las coordenadas de los objetos del primer conjunto y sean (x0 i, y0 i, z0 i) las del segundo. Las formas de emparejar ambos conjuntos son muy diversas. Una posibilidad es minimizar la suma de las diferencias al cuadrado entre los puntos emparejados. Esta idea se adecua bastante a la realidad, pues penaliza la posibilidad de emparejar objetos a gran distancia. Si la diferencia de tiempo en la que tomamos las localizaciones es suficientemente pequeña, el emparejamiento óptimo será probablemente correcto. Si denotamos por {1, . . . , p}al primer conjunto de puntos y por {10, . . . , p0} al segundo, entonces un arco (i, j0)tendría un coste asociado cij0= (xi−x0 j)2+ (yi−y0 j)2+ (zi−z0 j)2 1.4.8. Orden óptimo del inventario En otros contextos, necesitamos almacenar objetos que pierden o ganan valor con el paso del tiempo. Supongamos que tenemos una pila consistente en pobjetos del mismo tipo. Cuando hablamos de una pila, nos referimos al concepto de lista ordenada de manera que solo podemos sacar el último elemento en ser apilado cada vez. Cada objeto itiene una edad ai. Una función v(t)nos da la utilidad (o valor) esperada para un objeto de edad tcuando lo retiramos del almacén. Tenemos que cumplir con un determinado horario que especifica la hora en la que necesitamos cada elemento. El problema consiste pues en encontrar el orden de emisión de los artículos que maximiza la suma de las utilidades esperadas en los partículos. Un ejemplo de este tipo situaciones se da cuando almacenamos varias tinas de alcohol, ya que al ser éste un líquido volátil la evaporación del mismo deprecia el valor de cada tina con el tiempo. Algunas instancias de este problema son particularmente fáciles de resolver. Por ejemplo, cuando la función v(t)es convexa o cóncava. Cuando tenemos una función v(t)arbitraria, entonces podemos resolver el problema como problema de asignación. Sean t1, . . . , tp. Dado que cada objeto itiene una edad aien un tiempo cero, la utilidad 32 2.2. Fundamentos del algoritmo A partir de la definición de transformación admisible podemos probar un criterio de optimalidad para que un camino de aumento sea de longitud mínima. Teorema 2.2.3. Sea Tuna transformación admisible y P0∈P(M)con l0(P0) = 0 y d(P0)≤d(P)∀P∈P(M). Entonces P0es el camino de aumento con menor longitud con respecto a la función de coste cy matching M. Demostración. Se deduce inmediatamente a partir de la definición y las hipótesis a través de la siguiente cadena de igualdades / desigualdades. l(P0) = l0(P0) + d(P0) =d(P0)≤d(P)≤l0(P) + d(P) = l(P)∀P∈P(M) Definición 2.2.4. Sea B={S⊂V| |S|es impar, |S| ≥ 3}. A cada u∈Vy S∈Bpodemos asociar unos pesos yu∈ReyS∈R≥0. Definimos la transformación T:cuv →c0 uv tal que c0 uv =cuv −yu−yv−X S∈B,(u,v)∈δ(S) yS Nota 2.2.1. Si consideramos yuna solución factible del problema dual del matching, la transformación anterior asocia a los costes de cada arista los costes reducidos asociados a y. Esta transformación pues puede establecerse a través de las condiciones de holgura complementaria que si Sno pertenece a una familia de contracción, yS= 0. Definición 2.2.5. En las condiciones de la definición anterior, el grafo de admisibilidad G0= (V, E0)es aquel en el que E0={(u, v)∈E|c0 uv = 0}. Definición 2.2.6. Una transformación admisible Tdel tipo definida en la Definición 2.2.4 se dice adecuada con respecto a Msi existe una familia de contracción Aen G0 tal que (u, v)∈M⇒c0 uv = 0 yS>0⇒S∈Ay|M∩γ(S)|=1 2(|S| − 1) Con estas definiciones en mente pasamos a enunciar y probar el criterio de optimalidad que utilizaremos en el nuestro algoritmo. Teorema 2.2.4. Sea Mun matching y sea Tuna transformación adecuada. Sea A= {S∈B|yS>0}. Si existe P∈P(M)tal que c0(P) = 0 yPAes un camino de aumento respecto MA, entonces Pes el camino alternante de menor costo en P(M). Capítulo 2. Algoritmo SAP 33 Demostración. Sea s∈V, consideremos todos los caminos de aumento para Mque tienen como nodo inicial s. Sea P0un camino de este tipo. Sea S∈A, entonces s∈S⇒|(P0∩M)∩δ(S)|= 0 y |(P0\M)∩δ(S)|= 1 s /∈S⇒|(P0\M)∩δ(S)|−|(P0∩M)∩δ(S)| ≥ 0 En particular, l(P0)≥ys+Ps∈SyS. Dado que que PAes un camino de aumento con respecto a MA, se verifica s∈S⇒|(P\M)∩δ(S)| s /∈S⇒|(P\M)∩δ(S)|=|(P∩M)∩δ(S)|= 1 y, por tanto, l(P) = ys+Ps∈SyS. Aplicando el Teorema 2.2.3 tenemos el resultado. 2.3. Algoritmo del camino de aumento más corto Durante el procedimiento un emparejamiento M, una familia de contracción Ay unos peso yu(u∈V) e yS(Sconjunto impar de vértices) definen una transformación admisible T:cuv →c0 uv que verifica además que c0 uv = 0 si (u, v)∈M. Todos los pasos del algoritmo están constituidos en el grafo G×A. En [8] se utilizan distintas notaciones para los costes (πvpara los nodos e ySpara los conjuntos impares). Sin embargo, los costes no nulos ySasociados con un pseudonodo vSrepresenta ya de por sí al conjunto S∈A. Que estos costes puedan entenderse como asociados a un pseudonodo justifica que nosotros vayamos a utilizar la notación ypara ambos tipos. Antes de dar el algoritmo es preciso notar que utilizaremos un etiquetado de manera que todo nodo de G×Atendrá tres etiquetas temporales d+ v, d− vyp(v)con la información: d+ ves el límite superior de la longitud del camino alternante más corto que conecta vcon un cierto nodo sy tiene un número par de ejes. d− ves análogo para un número impar de ejes. p(v)es el predecesor de ven el camino alternante que define d− v. Estas distancias y medidas están tomadas respecto de los costes c0 uv. Pasemos a describir el algoritmo. Inicialización: yu:= yS:= 0 ∀u∈V, S ∈B M:= A:= ∅ c0 uv =cuv ∀(u, v)∈E Paso 1: Determinar si existe un nodo sno emparejado en G×Acon respecto a MA. Si no existe ninguno, ir al Paso 8. En caso contrario, ir al Paso 2. 34 2.3. Algoritmo del camino de aumento más corto Paso 2: [Inicialización de una fase.] d+ v:= ∞ ∀j∈G×A, v 6=s d+ s:= 0 d− v:= c0 sv, p(v) := s∀j∈G×A, v 6=s d− s:= ∞ Etiquetamos scon una Sy avanzamos al Paso 3. Paso 3: Determinamos: δ1= m´ın{d− v|vno está etiqetado ni emparejado}, δ2= m´ın{d− v|vno está etiquedato pero sí emparejado}, δ3= m´ın{1 2(δ+ v+δ− v)|vestá S-etiquetado} δ4= m´ın{yB+d− B|vBes un pseudonodo T-etiquetado} δ= m´ın{δ1, δ2, δ3, δ4} Si δ=δ1ir al Paso 7 Si δ=δ2ir al Paso 4 Si δ=δ3ir al Paso 5 Si δ=δ4ir al Paso 6 Paso 4: [Crecimiendo de un arbol alternante.] Sea uel nodo que define δy sea (u, v)∈MA. Entonces el árbol crece asignando un T-etiquetado a uy un Setiquetado a v. Establecemos d+ v:= δy escaneamos v, es decir, d− k:= m´ın{d− k, d+ j+c0 vk} ∀k6=u, v p(k) := vsi d− k=d+ v+c0 vk Ir al Paso 3 Paso 5: [Contracción de un blossom.] Sea uel nodo que define δy sea v:= p(u). Introduciendo la arista (u, v) al arbol alternante genera un blossom B. Este blossom se contrae en un pseudono vB. Establecemos: A:= A− {B}, d+ B:= δ, yB:= 0 y asignamos una S-etiqueda a vB Para todo i∈B S-etiqueado yi:= yi+ (δ−d+ i) c0 ij := c0 ij −(δ−d+ i)∀j6=i Para todo i∈B T-etiquedado yi:= yi+ (d− i−δ) c0 ij := c0 ij −(d− i−δ)∀j6=i Escaneamos vB(como en el Paso 4 para el nodo v) y vamos al Paso 3. Capítulo 2. Algoritmo SAP 35 Paso 6: [Expansión de un pseudonodo.] Sea vBel pseudonodo que define δ. Establecemos c0 uv := c0 uv +yB∀(u, v)∈δ(B) yB:= 0 y extendemos el matching MA. El blossom Bse particiona en dos caminos, uno tiene un número par de ejes y el otro un número impar Todos los nodos vdel camino par se etiquetan con SyTalternativamente yd+ v:= d− v:= δ. Los nodos S-etiqueados son escaneados. Los nodos u del camino impar no se etiquetan. Establecemos d− v= m´ın{d+ i+c0 iu |iestá Setiquetado} yp(v) := isi el mínimo se alcanza en i. Si no existe, d− v:= ∞. Finalmente, A:= A\ {B}. Vamos al Paso 3. Paso 7: [Aumento.] Hemos detectado un camino de aumento PAcon respecto a MA. Actualizamos MA=MA⊕PAy definimos y definimos para todos los (pseudo)nodos S-etiquetados v: yv:= yv+ (δ−d+ v) c0 uv := c0 uv −(δ−d+ i)u6=v y para los T-etiquedados yv:= yv−(δ−d− i) c0 uv =c0 uv + (δ−d− i) Eliminamos todas las etiquetas y vamos al Paso 1 Paso 8: Expandimos todos los pseudonodos y extendemos el matching MAa un matching perfecto Men G. Este matching es óptimo. 36 2.4. Método de reoptimización 2.4. Método de reoptimización En esta sección vamos a presentar un método de reoptimización para el problema del MCPM y cómo aplicarlo al problema biobjetivo. Antes de presentar este método hemos revisado diversos análsis de postoptimalidad como el de Ball y Taberna [5] o el de Weber [12], pero el que presenta mejores resultados computacionales y supone una simplificación importante en el número de casos a considerar es el que presenta Derigs. Nuestro objetivo en esta sección es, suponiendo que tenemos un matching óptimo Mpara el MCPM, vamos a alterar el coste de ciertos ejes. Denotemos cnew uv los costes nuevos para estas aristas y sea c0new uv los correspondientes costes reducidos respecto de una solución factible dual y. Antes de proseguir, precisemos algunas definiciones. Definición 2.4.1. Sea yuna solución factible del dual descrito gracias al Teorema 2.1.3. Consideremos c0 uv := c0 uv(y), es decir, los costes reducidos asociados a y. Consideremos G0el grafo de admisibilidad asociado. Diremos que yes fuértemente dual factible si yS>0implica que Ses contractible en G0. Definición 2.4.2. En las mismas condiciones de la definición anterior, diremos que (M, y)es un par compatible si yes fuertemente dual factible y se verifica (u, v)∈M⇒c0 uv = 0 yS>0⇒ |M∩γ(S)|=1 2(|S| − 1) Se puede probar que el método que hemos establecido anteriormente comienza con un par compatible (M, y)y encuentra el menor camino de aumento Py una solución fuertemente dual factible y0de manera que (M⊕P, y0)es un par compatible. Proposición 2.4.1. Sea Mun matching óptimo y sea u∈Vun nodo arbitrario. Entonces existe una solución yfuertamente dual factible tal que (M, y)es un par compatible y además yS= 0 ∀S∈A(G)tal que u∈S. Demostración. La demostración es constructiva. Podemos suponer que tenemos (M, y) un par compatible obtenido por el algoritmo SAP. Si yverifica la propiedad adicional, no hay nada que demostrar. En caso contrario, sea vtal que (u, v)∈M. Vamos a construir un grafo Gintroduciendo dos nuevos nodos artificiales r, t /∈Vy aristas (r, u)y(v, t)con costes cru := yu+X u∈S yScvt =yv+X u∈S yS+X v∈S yS Además, extendemos yal nuevo grafo Gde manera que yr:= yt:= 0 eyS= 0 ∀S∈ A(G)tal que S∩{r, t} 6=∅. Esta extensión hace que (M, y)sea un par compatible en G. Aplicamos ahora el algoritmo del SAP. Como Mera óptimo en G, el único camino Capítulo 2. Algoritmo SAP 37 de aumento que puede generarse es {(r, u),(u, v),(v, t)}. Dado que cvt >Pu∈SyS, antes de que se detecte este camino la solución dual será alterada a y0de manera que y0 S= 0 ∀S∈A(G)tal que u∈S. Obtenemos por tanto un par (M, y0)que verifica la propiedad deseada. Es importante hacer una consideración para entender por qué es importante el teorema anterior. Sea yuna solución fuertemente dual factible, definimos A(y) = {S∈ A(G)|yS>0}. Naturalmente, si existe algún S∈A(G)con v∈Stal que yS>0, el nodo vno estará -estará contraído en algún pseudonodo- en el grafo G×A(y). Por tanto, siguiendo el proceso de la demostración anterior, siempre podemos encontrar un par compatible (M, y)de manera que un determinado vértice vesté en el grafo G×A(y). Pasamos a ver el algoritmo de repotimización. Supongamos que todas las aristas cuyos costes son alteradas están en δ(v). Naturalmente, si c0new ij ≥0∀(u, v)∈δ(v), nuestro matching sigue siendo óptimo, luego podemos suponer que hay alguno que es negativo. Inicialización: Sea (M, y)un par compatible óptimo para MCPM Paso 1: Si vno está contraído en G×A(y)vamos al Paso 2. En otro caso, aplicamos el procedimiento anterior y denotamos y:= y0. Vamos al Paso 2. Paso 2: Sea (u, v)∈Mestablecemos M:= M\ {(u, v)} δ:= m´ın{c0new vk |(v, k)∈δ(v)}<0 cvk := cnew vk ∀(v, k)∈δ(v) yv:= yv+δ Ir al Paso 3. Paso 3: Aplicar el algoritmo SAP a (M, y)y determinar el menor camino de aumento Py una solución dual y0tal que el par (M⊕P, y0)sea un par compatible A continuación explicamos cómo aplicar el oráculo anterior para reoptimizar en general un problema de matching y añadimos además una forma algorítima del método. El proceso de reoptimización que muestra Derigs solo nos permite alterar los costes de las aristas relativas a un vértice en concreto. Por tanto, tenemos que invocarlo secuencialmente. Partimos de un par compatible (M, y)óptimo. Calculamos V0⊂V tal que δ(V0)∪γ(V0) = E. Enumeramos los vértices de V0={u1, . . . , us}. En el primer paso, llamamos al método de repotimización inicializandolo mediante (M, y), el vértice u1y los costes cnew uv ∀(u, v)∈δ(u1). El oráculo devuelve un par compatible (M1, y1). En el paso general, partimos del par compatble (Mi, yi)y tomamos el vértice un+1 e inicializamos 38 2.4. Método de reoptimización el proceso de repotimización junto con los costes cnew uv ∀(u, v)∈δ(ui+1). De manera algoritmica, podemos escribir: Inicialización: Sea (M, y)un par compatible óptimo para MCPM y sea V0={u1, . . . , us}un conjunto de vértices tales que δ(V0)∪γ(V0) = E. Paso 1: Invocamos el proceso de repotimización con (M,y), el vértice u1 y la función de costes k1, que mantiene los mismos costes k1 uv =cuv ∀(u, v)/∈δ(u)yk1 uv =cnew uv ∀(u, v)∈δ(u). Obtenemos un par compatible óptimo (M1, y1). Avanzamos al Paso 2 Paso 2: Para todo n= 2, . . . , s Invocamos el proceso de repotimización partiendo de el par compatible (Mn−1, yn−1), el vértice uny la función de costes kndefinida como kn uv =kn−1 uv si (u, v)/∈δ(un)y kn uv =cnew uv si (u, v)∈δ(un). Al final de cada iteración obtenemos un par compatible (Mn, yn)con el cuál se inicializa la siguiente iteración, hasta que n=s. En este caso avanzamos al Paso 3. Paso 3: El par compatible (Ms, ys)es óptimo para Gcon la función de coste cnew Notemos que es posible que algunas o todas las aristas cuyos costes alteramos en un paso general pueden haber sido alteras en otros anteriores, pero tiene un efecto inocuo en el desarrollo del algoritmo. Los conjuntos Wque verifican que δ(W)∪γ(W) = Ese denominan coberturas por vértices de G. En general, el problema de decidir si un grafo Gtiene una cobertura por vértices de tamaño kes NP-completo. A pesar de ello, por la observación anterior, tampoco es estrictamente necesario que obtengamos una cobertura de tamaño mínimo para inicializar el algoritmo. 2.4.1. Problema biobjetivo En el matching paramétrico, tal y como vimos en la Subsección 1.1.1, la función de coste es cuv(λ) = c1 uv +λc2 uv. Por tanto, para aplicar el procedimiento anterior, tenemos que considerar cnew uv =c1 uv + (λ+ ∆)c2 uv =cuv(λ+ ∆). Notemos que ∆no tendría por qué ser un parámetro real, sino que podríamos considerar ∆∈R|E|. En ese caso, cnew uv =c1 uv + (λ+ ∆uv)c2 uv Capítulo 2. Algoritmo SAP 39 2.5. Manteniendo optimalidad Estamos especialmente interesados en el caso del matching biobjetivo en saber cuánto puede variar ∆de manera que la base óptima que representa a nuestro matching Mcontinúe siendo óptima. Recordemos que aunque la base deje de ser óptima, en principio la solución podría seguir siéndolo, pues una solución puede estar representada por más de una base. Para realizar este análisis, consideremos el problema de Programación Lineal: m´ın xX (u,v)∈E xuv(c1 uv +λc2 uv) s.a. X (u,v)∈δ(u) xuv ≤1,∀u∈V X (u,v)∈δ(S) xuv ≥1∀S∈A(G) xuv ≥0∀(u, v)∈E Sea Ala matriz de restricciones sin considerar la introducción de las variables de holgura. Si |V|= 2n, una base Bóptima tendrá, en primer lugar, ncolumnas de A correspondientes a las aristas del matching según el indexado que hayamos hecho de las mismas; y otras scolumnas de la base canónica para las holguras cuyas desigualdades se verifiquen de manera estricta. Sean i1, . . . , inlos índices asociados a las aristas del matching, sean j1, . . . , jklos índices asociados a dichas variables de holgura, sea Nel número total de restricciones del problema tras la introducción de las holguras y denotemos por epel p-ésimo vector de la base canónica de RN, entonces la matriz B es precisamente B=A·i1· · · A·inej1· · · ejk Si denotamos por bla columna de términos independientes y por Nlas columnas no básicas de A, sabemos que la condiciones que certifican que Bes una base óptima son B−1b≥0 cRt=ct N−ct BB−1N≥0 2.5.1. Problema biobjetivo Realizamos el análisis para el caso en que λes un único parámetro escalar común. Alterar λno modifica en ningún caso la primera condición, luego nos centraremos solo en la segunda. cR(λ+ ∆)t=cN(λ+ ∆)t−cB(λ+ ∆)tB−1N =cN(λ)t−cB(λ)tB−1N+ ∆c2 N t−∆c2 B tB−1N =cR(λ)t+ ∆c2 R t≥0 40 2.5. Manteniendo optimalidad Si simplemente utilizamos la aproximación clásica a la reoptimización de un problema de Programación Lineal, sabemos que Bse mantiene como base óptima siempre que m´ax (u,v):c2uv >0 −c(λ)uv c2uv ≤∆≤m´ın j:c2uv <0 −c(λ)uv c2uv Sea yuna solución del problema dual asociado para λ, entonces sabemos que para cualquier arista (u, v)∈E, en particular para una arista representada por una variable no básica, el coste reducido es precisamente c(λ)uv =cuv +λ−yu−yv−X S∈B,(u,v)∈δ(S) yS Una vez estudiado el caso anterior, consideremos un caso más genérico en el que ∆∈R|E|. Para simplificar la notación, sea λ∈Ry∆∈R|E|podemos escribir c(λ+ ∆), donde c(λ+ ∆)uv =c1 uv + (λ+ ∆uv)c2 uv y análogamente para c2. Podemos considerar entonces cR(λ+ ∆)t=cN(λ+ ∆N)t−cB(λ+ ∆B)tB−1N =cN(λ)t−c2 B(λ)tB−1N+cN(∆N)t−c2 B(∆B)tB−1N =cR(λ)t+cR2(∆)t≥0 Nuevamente, podemos escribir los costes reducidos con respecto a λusando una solución dual asociada. Cada variable no básica nos da una desigualdad sobre ∆, de manera que tenemos un poliedro en R|E|de valores admisibles. En ambos casos, al llegar a uno de los vértices del poliedro (el extremo del intervalo en el caso más simple), podemos escoger una nueva base y repetir el análisis. Capítulo 3 Análisis computacional Presentamos a continuación la parte computacional de nuestro trabajo. En primer lugar incluimos las distintas partes de la implementaciones de los algoritmos para pasar después a los resultados obtenidos. 3.1. Implementación en Python Todas las implementaciones están hechas utilizando la distribución de Python 3.7 proporcionada por Anaconda. Además, para la resolución de problemas de Programación Lineal utilizamos el software Gurobi 8.1.0, a cuyo oráculo accedemos a través de Python. Antes de comenzar presentamos un codigo con algunas funciones auxiliares que utilizamos. 1# En este codigo pasamos a escribir funciones auxilares para el resto del trabajo. 2 3import networkx as nx 4import numpy as np 5import operator 6import random 7import itertools 8 9# Construimos una funcion que comprueba si un vector es de coordenadas enteras 10 def entero(s): 11 a = map(lambda x: x == int(x),s) 12 return(all(a)) 13 14 # Dado un grafo G y un subcojunto de aristas F con capacidad , construimos una funcion que nos de la suma de las capacidades de F. 41 48 3.1. Implementación en Python 91 card = (len(nodos)−1)/2 92 return(suma <= card) 93 94 # Generamos las componentes adecuadas 95 conE = [a for ain list(nx.connected_components(GE)) if len(a) %2== 1and comprueba(a)] 96 97 # Si las lista con es no vacia , hemos encontrado planos de cortes , los añadimos a nuestro modelo y volvemos a comprobar. 98 if (conN != []): 99 for ain conN: 100 S = list(GN.edges(a)) 101 card = len(a)−1 102 m.addConstr(suma(S) <= card/2) 103 m.optimize () 104 print("Hemos usado H1") 105 elif (conE != []): 106 for ain conE: 107 S = list(G.edges(a)) 108 card = len(a)−1 109 m.addConstr(suma(S) <= card/2) 110 m.optimize () 111 print("Hemos usado H2") 112 # Paso 7.3: En el caso de que las heuristicas no hayan funcionado , construimos un plano de corte mediante el procedimiento de Padberg y Rao. 113 else: 114 W,F = padraomod(nx.Graph(G),sols) 115 m.addConstr(suma(W.difference(F))−suma(F) >= 1−len (F)) 116 m.optimize () 117 print("Hemos usado Padberg−Rao") 118 # Actualizamos el vector solucion y las iteraciones. 119 sols = [m.getVars()[i].X for iin range(E)] 120 itera = itera + m.IterCount 121 end = time.time() 122 return(m, itera , end−start) Capítulo 3. Análisis computacional 49 3.1.3. Método del simplex sobre el poliedro del matching El siguiente código veremos una forma sencilla de programar en Python (y resolver con Gurobi) el problema del emparajamiento usando el método del simplex sobre el poliedro del matching. 1# En este codigo resolvemos el problema biobjetivo del emparejamiento con funciones de coste c1, c2y parametro lambda , utilizando el simplex sobre el poliedro del matching tal y como lo describio Edmonds. 2 3from gurobipy import ∗ 4from funcionaux import ∗ 5import time 6 7 8# Pasamos a definir la funcion entero , que resuelve el problema del emparejamiento con parametros 9# G: Un grafo. 10 # lamb: Valor de lambda. 11 # c1,c2: Funciones de coste 12 13 def matching(G, c1, c2, lamb = 0): 14 15 start = time.time() 16 17 # Creamos el modelo. 18 m = Model("biobjetivo"); 19 20 21 # Algunos parametros del grafo 22 N = len(G) 23 E = len(G.edges()) 24 25 # Creamos las variables. 26 x = m.addVars(list(G.edges()), vtype=’C’, name = "x"); 27 28 # De esta forma , generamos las dos funciones objetivo y procedemos de manera analoga a los codigos anteriores. 29 cos1= objind(c1) 30 cos2= objind(c2) 31 32 # Definimos el sentido de la optimizacion. 33 m.ModelSense = GRB.MINIMIZE 34 35 # Añadimos la funcion objetivo 50 3.1. Implementación en Python 36 m.setObjective((1−lamb)∗cos1+lamb∗cos2); 37 38 # Generamos las restricciones del problema entero 39 for vin G: 40 m.addConstr(suma(delta(G,[v])) == 1, name = "delta" + str(v)); 41 42 for iin range(3,len(G)+1,2): 43 conj1=set(itertools.combinations(set(G.nodes()), i)) 44 for ain conj1: 45 S = G.subgraph(a) 46 m.addConstr(suma(set(S.edges)) <= (len(a)−1)/2) 47 48 # Resolvemos el prolbema. 49 m.optimize () 50 51 # Guardamos el numero de iteraciones y el tiempo en el que termina el algoritmo 52 itera = m.IterCount 53 end = time.time() 54 return(m, itera , end−start) 3.1.4. Resolución como problema de Programación Entera Para tener un tercer método con el que comparar, también realizamos una implementación del problema del matching en su formulación como problema de Programación Entera. 1# En este codigo resolvemos con Gurobi el problema biobjetivo del emparejamiento con funciones de coste c1, c2y parametro lambda utlizando la formulacion del problema como problema de Progrmacion Entera. 2 3from gurobipy import ∗ 4from funcionaux import ∗ 5import time 6 7# Pasamos a definir la funcion entero , que resuelve el problema del emparejamiento parametro con parametros 8# G: Un grafo. 9# lamb: Valor de lambda. 10 # c1,c2: Funciones de coste 11 12 def integer(G, c1, c2, lamb = 0): Capítulo 3. Análisis computacional 51 13 # Comenzamos capturando el tiempo de inicio 14 start = time.time() 15 16 # Creamos el modelo 17 m = Model("biobjetivo"); 18 19 # Algunos parametros del grafo 20 N = len(G) 21 E = len(G.edges()) 22 23 # Creamos las variables. 24 x = m.addVars(list(G.edges()), vtype=’B’, name = "x"); 25 26 # De esta forma , generamos las dos funciones objetivo. Notemos que escribimos directamente los costes sobre las variables aristas del modelo , pues si una arista no esta en el grafo , no tiene sentido generar una variable para ella. 27 cos1= objind(c1) 28 cos2= objind(c2) 29 30 # Definimos la funcion objetivo y queremos maximizar/ minimizar. 31 m.ModelSense = GRB.MINIMIZE 32 33 # Añadimos la funcion objetivo 34 m.setObjective((1−lamb)∗cos1+lamb∗cos2); 35 36 # Generamos las restricciones del problema entero 37 for vin G: 38 m.addConstr(suma(delta(G,[v])) == 1, name = "delta" + str(v)); 39 40 # Resolvemos el prolbema. 41 m.optimize () 42 43 # Guardamos el numero de iteraciones y de exploraciones de branch and bound 44 itera = m.IterCount 45 bac = m.NodeCount 46 47 end = time.time() 48 return(m, itera , bac , end−start) 52 3.1. Implementación en Python 3.1.5. Grafos de Erdös-Rényi En el siguiente código vamos a escribir una función que nos compare sobre grafos generados aleatoriamente con una cantidad fija de vértices los métodos que hemos programados anteriormente. Para generar los grafos aletoriamente utilizamos el comando de la librería networkx 1networkx.gnp_random_graph(M,p) Esta librería de grafos sobre la que hemos apoyado parte del código utiliza el método de Erdös-Rényi, introducido por dichos matemáticos húngaros en 1959, para generar grafos aleatorios con dos parámetros. El concepto es el siguiente: Supongamos que tenemos una red con Mnodos distribuidos aleatoriamente y desconectados. Ahora cada uno de las posibles aristas del grafo completo KMse añade a nuestra red de manera independiente con probabilidad p. Naturalmente cuanto mayor sea pmás probable es que el grafo resultante sea denso. Si pes demasiado pequeño puede ocurrir que nuestro grafo sea diconexo. Nuestra función tiene como entrada el número Nde nodos de los grafos que queremos explorar. Es decir, todos los grafos generados tendrán Nvértices. El parámetro Ideterminar la cantidad de grafos (y por tanto, de problemas) que se van a generar. Finalmente, utilizamos dos parámetros para establecer si queremos o no utilizar métodos alternativos al algoritmo de Grötschel-Holland. 1# A continuacion generamos grafos aleatorios y ponemos a prueba los metodos y comparamos resultados. 2 3from holland import ∗ 4from integer import ∗ 5from matching import ∗ 6import networkx as nx 7import numpy as np 8import time 9import pandas as pd 10 11 # La siguiente funcion nos devuelve el tiempo e iteraciones medias de los distintos metodos que hemos programado para funciones y grafos aleatorios con un numero de nodos N y utilizando I casos. 12 13 def analisis(N, I, par1=0, par2=0): 14 iter0,iter1,iter2=0,0,0 15 t0, t1, t2=0,0,0 16 bcut = 0 17 aristas = 0 18 for iin range(I): Capítulo 3. Análisis computacional 53 19 #Construimos un grafo aleatorio 20 H = nx.gnp_random_graph(N,0.5) 21 E = len(H.edges()) 22 23 # Construimos funciones de coste aleatorias 24 c1= np.random.rand(E) 25 c2= np.random.rand(E) 26 27 # Optimizamos utilizando Groetschel−Holland 28 modelo0, simplex0, tiempo0= optimiza(H,c1,c2,0) 29 30 # Actualizamos los parametros 31 iter0=iter0+ simplex0 32 t0= t0+ tiempo0 33 aristas = aristas + E 34 # Si hemos activado la Programacio Entera 35 if (par1== 1): 36 # Optimizamos usando Programacion Entera 37 38 modelo1, simplex1, bcut1, tiempo1= integer(H,c1,c 2,lamb) 39 40 # Actualizamos los parametros 41 iter1=iter1+ simplex1 42 t1= t1+ tiempo1 43 bcut = bcut + bcut1 44 45 if (par2== 1): 46 # Utilizamos el simplex normal 47 modelo2, simplex2, tiempo2= matching(H,c1,c2,lamb ) 48 49 # Actualizamos los parametros 50 iter2=iter2+ simplex2 51 t2= t2+ tiempo2 52 53 return([N, I, aristas/I, iter0/I, iter1/I, iter2/I, 54 t0/I, t1/I, t2/I, bcut/I]) 54 3.1. Implementación en Python 3.1.6. Reoptimización Pasamos ahora a presentar el codigo con el que vamos a realizar el análisis de postoptimalidad. En primer lugar, debemos saber que cuando optimizamos un modelo de Programación Lineal en Gurobi, éste posee una gran variedad de parámetros y atributos que podemos consultar. Dos de ellos son SAObjLow ySAObjUp, que para cada componente de la función objetivo nos indica hasta dónde podemos variar hacia abajo y hacia arriba respectivamente de manera que la base continúe siendo óptima. Supongamos que hemos resuelto el modelo para un cierto valor de λ∗, y denotemos por lyua los vectores definidos anteriormente para este valor λ∗. Nuestro objetivo es encontrar el máximo valor de λde manera que la base obtenida siga siendo óptima, es decir, queremos resolver el problema, dada By funciones objetivos ci, m´ax λ s.a. l ≤c1+λc2 u≥c1+λc2 λ≥λ∗ Una vez obtenido este valor de λpodemos alterar el modelo ya construido de Gurobi y cambiar la función objetivo actualizando el valor de λ, a la que sumamos un pequeño εpara que realmente optimice de nuevo, y comprobamos si aun tenemos la misma solución -pues una misma solución puede estar representada por distintas bases- o ésta ha variado. 1# Para el estudio de postoptimalidad tenemos primero que definir la siguiente funcion que resuelve un modelo lineal con ciertos parametros tal y como hemos comentado anteriormente. 2def buscalamb(lam ,c1,c2,u,l): 3m = Model("lambda"); 4m.ModelSense = GRB.MAXIMIZE 5x = m.addVars(1,lb = lam , ub = 1, vtype=’C’, name = "x") 6m.setObjective(x[0]) 7for iin range(len(c1)): 8m.addConstr(l[i] <= (1−x[0])∗c1[i]+(x[0])∗c2[i]); 9m.addConstr(u[i] >= (1−x[0])∗c1[i]+(x[0])∗c2[i]) 10 m.optimize () 11 return(m) 12 13 14 # Construimos una funcion que dado un grafo y dos funciones de coste resuelva el problema biobjetivo para lambda = 0y encuentre los valores lambda donde deja de ser optima y una nueva solucion , hasta llegar a lambda = 1. Capítulo 3. Análisis computacional 55 15 def damelam(G,c1,c2): 16 17 # Inicializamos parametros 18 todos = [0] 19 cambian = [0] 20 lamact = 0 21 sols = [] 22 23 # Resolvemos el modelo para lambda = 0. 24 a1, b1, h1= optimiza(G,c1,c2,0) 25 26 # Desactivamos el presolve y gardamos la primera solcion 27 a1.Params.presolve = 0 28 sols.append([a.X for ain a1.getVars()]) 29 30 # Utilizamos un parametro para guardar 31 q=0 32 while(lamact <1and q == 0): 33 if (lamact == 1): 34 q=1 35 36 # Capturamos parametros del modelo 37 l=a1.getAttr("SAObjLow") 38 u=a1.getAttr("SAObjUp") 39 H1= [a.X for ain a1.getVars()] 40 lamact = buscalamb(lamact ,c1,c2,u,l).ObjVal 41 for i,a in enumerate(a1.getVars()): 42 a.Obj = c1[i] + (1.03)∗lamact∗c2[i] 43 a1.optimize() 44 H2= [a.X for ain a1.getVars()] 45 if (H1!= H2): 46 cambian.append(lamact) 47 sols.append([a.X for ain a1.getVars()]) 48 todos.append(lamact) 49 print() 50 return(todos ,cambian , sols) 56 3.1. Implementación en Python 3.1.7. Dibujando matchings En L A T EX es común utilizar, tal y como hemos hecho nosotros, el paquete tikz para dibujar grafos (entre otras muchas funcionalidades). Sin embargo, no existe una forma sencilla de dibujar grafos relativamente grandes a partir, por ejemplo, de su matriz de adyacencia. La librería networkx tiene distintas formas de dibujar grafos apoyándose en otra librería conocida de Python, matplotlib. Vamos a utilizar estas librerías para representar algunos grafos y respectivos matching perfectos de coste mínimo obtenido por el algoritmo de Grötschel-Holland. En general todos los comandos para dibujar grafos pueden colorear de manera independiente cada nodo o eje, además de una gran cantidad de parámetros. 1# Este comando nos permite dibujar el grafo situando los nodos de manera dispersa sobre una elipse. 2networx.draw_circular(G) Figura 3.1 Figura 3.2 Capítulo 3. Análisis computacional 57 1# Si simplemente utilizamos draw , el dibujo es algo mas caotico. 2networx.draw(G) Figura 3.3 1# Tambien tenemos una opcion para utilziar el algoritmo de Kamada−Kawai para dibujar redes. 2networx.draw_kamada_kawai(G) Figura 3.4 64 3.3. Reoptimización 0 1 2 3 4 5 6 7 Figura 3.9: Solución para λ= 0 0 1 2 3 4 5 6 7 Figura 3.10: Solución para λ= 0,5814219436621 0 1 2 3 4 5 6 7 Figura 3.11: Solución para λ= 0,6626459613176332 Capítulo 3. Análisis computacional 65 Reportamos también el siguiente análisis. Para un cierto número Nde vértices, generamos 5 grafos y sendas funciones de costo. A continuación usando las mismas técnicas que usamos anteriormente, calculamos el número de soluciones soportadas, es decir, aquellas que están en la frontera del poliedro y que conseguimos obteniendo al variar λ∈[0,1]. Mostramos el número promedio de soluciones para cada N considerado. N Sols 20 10.00 40 19.80 60 31.00 80 45.60 100 51.60 120 63.80 140 68.00 Por desgracia es realmente costoso tratar de buscar las soluciones interiores para el problema biobjetivo, de manera que para grafos de 30 nodos ya resulta relativamente costoso en un ordenador usual. En nuestro caso, para ciertos valores de Nentre 16 y30, vamos a generar 5 problemas para estudiar el comportamiento medio de las soluciones soportadas frente a no soportadas. La función para generar grafos aleatorios toma en este caso el parámetro p= 0,7. En la siguiente tabla encontramos el número Nde vértices de los grafos considerados, la Media Total de soluciones (soportadas y no soportadas), la media de soluciones no soportadas y la media de las soluciones soportadas frente al total. N Media Total Media No Sop Media Porcentaje 14 9.40 4.60 45.39 16 15.80 10.20 63.72 18 18.80 10.40 55.90 20 22.60 14.80 65.31 22 30.00 19.20 63.55 24 39.20 27.00 68.81 26 40.00 29.40 72.81 28 53.80 38.60 73.05 30 63.00 48.40 76.91 Adjuntamos además algunas gráficas de la tabla anterior. 66 3.3. Reoptimización ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● 20 40 60 15 20 25 30 N Valor Leyenda ● ● Media.No.Sop Media.Total Figura 3.12: Número medio de soluciones frente a nodos ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● 60 80 100 15 20 25 30 N Media.Porcentaje Figura 3.13: Media de porcentajes frente a nodos Capítulo 4 Conclusiones A lo largo de este trabajo hemos analizado ampliamente el problema del emparejamiento, tanto en su formulación usual como en su formulación biobjetivo. De esta forma, hemos comprobado su versatilidad y utilidad tanto en el marco de otros problemas matemáticos como en insertados dentro de aplicaciones prácticas, destacando el problema de la asignación. En nuestra descripción del algoritmo del SAP hemos analizado conceptos teóricos que desarrollan la estructura intrínseca del emparejamiento en grafos, gracias a la cual se vertebra dicho método. Este es especialmente interesante porque permite realizar un proceso de reoptimización arista a arista invocando al oráculo de la optimización. Finalmente, al decidirnos establecer una primera aproximación a una implementación computacional del algoritmo de Grötschel-Holland, programamos distintos algoritmos en un entorno de Python. Nos ha servido para explorar las características de algunas librerías de Teoría de Grafos como networkx, comprender la generación de planos de corte gracias la cuál se resuelve el problema de separación mediante el algoritmo de Padberg-Rao o acercarnos a la distribución, y por tanto al comportamiento, de las soluciones soportadas y no soportadas en el poliedro entero del matching para el caso biobjetivo. 67 68 Bibliografía [1] Edmonds, J. (1965). Maximum matching and a polyhedron with 0, 1-vertices. Journal of Research of the National Bureau of Standards B, 69, 125–130. [2] Kuhn, H. W. (1955), The Hungarian method for the assignment problem. Naval Research Logistics, 2: 83-97. doi:10.1002/nav.3800020109 [3] Grötschel, M.; Holland, O. Mathematical Programming (1985) 33: 243. https://doi.org/10.1007/BF01584376 [4] Schrijver, Alexander. (1983). Short proofs on the matching polyhedron. Journal of Combinatorial Theory, Series B. 34. 104-108. 10.1016/0095-8956(83)90011-4. [5] Ball, M. O., Taverna, R. (1985). Sensitivity analysis for the matching problem and its use in solving matching problems with a single side constraint. Annals of Operations Research, 4(1), 25-56. [6] Papadimitriou, C.H.; Steiglitz, K. (1998), Combinatorial optimization: algorithms and complexity, Mineola, NY: Dover, pp.308-309. [7] Nicos Christofides, Worst-case analysis of a new heuristic for the travelling salesman problem, Report 388, Graduate School of Industrial Administration, CMU, 1976 [8] George L. Nemhauser and Laurence A. Wolsey. 1988. Integer and Combinatorial Optimization. Wiley-Interscience, New York, NY, USA. [9] J. Edmonds, “Paths, Trees, and Flowers,” Can. J. Math., 17,449-467 (1965). [10] Pulleyblank, William and Edmonds, Jack. (1970). Facets of 1-Matching Polyhedra. 10.1007/BFb0066196. [11] Derigs, U. (1981), A shortest augmenting path method for solving minimal perfect matching problems. Networks, 11: 379-390. doi:10.1002/net.3230110407 [12] Weber, G. M. (1981), Sensitivity analysis of optimal matchings. Networks, 11: 41-56. doi:10.1002/net.3230110105 69 70 Bibliografía [13] Padberg, Manfred W., and M. R. Rao. “Odd Minimum Cut-Sets and b- Matchings.” Mathematics of Operations Research, vol. 7, no. 1, 1982, pp. 67–80. JSTOR, www.jstor.org/stable/3689360. [14] Letchford, Adam; Reinelt, Gerhard; Theis, Dirk Oliver. (2004). A Faster Exact Separation Algorithm for Blossom Inequalities. 10.1007/978-3-540-25960-215. [15] M. Grötschel, L. Lovász, and A. Schrijver. The ellipsoid method and its consequences in combinatorial optimization. Combinatorica, 1(2):169-197, 1981 [16] R. Borndörfer, S. Schenker, M. Skutella, T. Strunk: PolySCIP. Mathematical Software - Proceedings of ICMS 2016, G.-M. Greuel, T. Koch, P. Paule, A. Sommese (Eds.), Lecture Notes in Computer Science Vol. 9725, ISBN: 978-3-319-42431-6 [17] Gurobi Optimization, Inc. Gurobi Optimizer Reference Manual Version 3.0. Houston, Texas: Gurobi Optimization, April 2010.