Estudio de las características y técnicas de resolución de los problemas de asignación multiperiodo

Universidad de Sevilla Facultad de Matem´ aticas Estudio de las caracter ´ ısticas y t´ ecnicas de resoluci´ on de los problemas de asignaci´ on multiperiodo Trabajo Fin de M´aster M´aster Universitario en Matem´aticas Autor: ´ Alvaro Bl´azquez Molino Resumen En este Trabajo Fin de M´aster trabajaremos con problemas de optimizaci´on combinatoria con variables enteras que pueden ser dif´ıciles de resolver con exactitud debido a que tienen un gran n´umero de variables o un conjunto de restricciones complicadas. Para hacer frente a estos problemas, se aplicar´a un m´etodo de descomposici´on conocido como “Relajaci´on Lagrangiana”. El objetivo de este trabajo involucra diferentes aspectos de la investigaci´on en optimizaci´on combinatoria. En primer lugar, se definir´an los principales conceptos asociados a los problemas de programaci´on matem´atica y algunos resultados. En segundo lugar, tambi´en se definir´a el concepto de Relajaci´on Lagrangiana de un problema y c´omo nos puede ayudar a acotar o encontrar la soluci´on del problema original. Finalmente, se estudiar´an los Problemas de Asignaci´on Multiper´ıodo, y se desarrollar´an y detallar´an algunos m´etodos y procedimientos de soluci´on para mostrar a trav´es de casos pr´acticos de estos problemas c´omo los conocidos m´etodos iterativos que utilizan la Relajaci´on Lagrangiana pueden ayudarnos a aproximarnos a la soluci´on del problema original en un l´ımite de tiempo. iii iv Abstract In this Master Thesis we will work with combinatorial optimization problems with integer variables that can be difficult to solve exactly due to having a large number of variables or a set of complicating restrictions. To deal with these problems, a decomposition method known as “Lagrangian Relaxation”will be applied. The aim of this work involves different aspects of research in combinatorial optimization. Firstly, the main concepts associated with mathematical programming problems and some results will be recalled. Secondly, the Lagrangian Relaxation concept of a problem will also be defined and how it can help us to bound or find the solution of the original problem. Finally, Multiperiod Assignment Problems will be studied, and some solution methods and procedures will be developed and detailed to show through practical cases of these problems how the well-known iterative methods that use Lagrangian Relaxation can help us to approximate the solution of the original problem in a reasonable time limit. v vi ´ Indice general Resumen III Abstract V ´ Indice VII Prefacio XI 1. Preliminares 1 1.1. Problemas y formulaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Cotas.......................................... 4 1.3. Relajaciones...................................... 6 1.4. CotasPrimales .................................... 7 1.5. Algoritmo del simplex y de punto interior . . . . . . . . . . . . . . . . . . . . . . 7 1.6. Algoritmos de ramificaci´on y poda . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.7. Algoritmos de planos de corte . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2. Dualidad Lagrangiana 11 2.1. Relajaci´on Lagrangiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2. La fuerza del Dual Lagrangiano . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3. Resoluci´on del Dual Lagrangiano . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4. Algunas consideraciones para el algoritmo de subgradiente . . . . . . . . . . . . 19 2.4.1. Una forma de seleccionar las longitudes de paso . . . . . . . . . . . . . . 19 2.4.2. Condiciones de parada . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.3. Divisi´on en subproblemas enteros o continuos . . . . . . . . . . . . . . . 21 3. Problemas de Asignaci´on Multiperiodo 23 3.1. Permutaciones y el problema de asignaci´on . . . . . . . . . . . . . . . . . . . . . 23 3.1.1. Aplicaciones de problemas de asignaci´on . . . . . . . . . . . . . . . . . . 25 3.2. Problema de asignaci´on Planar de 3 ´ındices . . . . . . . . . . . . . . . . . . . . . 26 3.3. Problema de asignaci´on Axial de 3 ´ındices . . . . . . . . . . . . . . . . . . . . . 29 3.4. El Dual Lagrangiano de los problemas de asignaci´on multiperiodo Planar y Axial 30 4. Estudio computacional 33 vii 4.1. Python y el solver Gurobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2. Instancias ....................................... 34 4.3. Cotas Primales iniciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.4. Cotas primales y duales del algoritmo exacto . . . . . . . . . . . . . . . . . . . . 36 4.5. Par´ametros iniciales del algoritmo de subgradiente . . . . . . . . . . . . . . . . . 37 4.6. Algoritmo de subgradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.7. Algoritmo de reparaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.8. Variantes de resoluci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.9. Resultados....................................... 40 4.9.1. Resoluci´on con algoritmo exacto . . . . . . . . . . . . . . . . . . . . . . . 40 4.9.2. Resoluci´on con algoritmo de subgradiente . . . . . . . . . . . . . . . . . . 43 Conclusiones 49 Lista de figuras 50 Lista de tablas 51 Bibliograf´ıa 54 viii Cap´ıtulo 1. Preliminares 3 mantiene para conjuntos no acotados enteros X={x:Ax ≤by entero}y para conjuntos enteros mixtos X={(x, y) : Ax +Gy ≤b, x ≥0 y entero, y≥0}con A, G, b racionales. En cualquier caso se trata de una soluci´on te´orica, ya que en la mayor parte de los casos no hay una caracterizaci´on sencilla para todas las desigualdades que describen conv(X). Veamos ahora en qu´e caso podemos decir que una formulaci´on P1es mejor que otra formulaci´on P2. Como la soluci´on ideal conv(X) satisface X⊆conv(X)⊆Ppara cualquier formulaci´on P, esto sugiere que la mejor formulaci´on ser´a aquella que m´as se ajuste a conv(X) lo que formalmente se expresar´ıa as´ı: Definici´on 5. Dado un conjunto X⊆Zn, y dos formulaciones P1yP2para X,P1es mejor formulaci´on que P2si P1⊆P2. Ejemplo 1. Consideremos el conjunto S={(0,0),(0,1),(1,0),(1,1)}. Dos formulaciones v´alidas para este conjunto son: P1 = {x∈R2:−x1≤0 −x2≤0 x1≤1 x2≤1} P2 = {x∈R2:−x1≤0 −x2≤0 x1≤1,5 x2≤1,5} Como P1⊆P2, entonces P1es mejor formulaci´on que P2para el conjunto S, v´ease la Figura 1.1 y la Figura 1.2. Para m´as informaci´on sobre estos problemas consultar [17] o [10]. 4 Secci´on 1.2. Cotas Figura 1.1: Formulaci´on P1para el conjunto S Figura 1.2: Formulaci´on P2para el conjunto S 1.2. Cotas Consideramos el problema (COP) (COP) Z= max {c(x) : x∈X⊆Zn}. Nos preguntamos c´omo se podr´ıa probar que un punto x∗es ´optimo. Trataremos de encontrar una cota inferior Z≤Zy una cota superior ¯ Z≤Zde forma que Z=¯ Z=Z. En la pr´actica, Cap´ıtulo 1. Preliminares 5 necesitaremos emplear un algoritmo que encuentre una sucesi´on decreciente de cotas superiores ¯ Z1>¯ Z2> ... > ZU≥Z y una sucesi´on creciente de cotas inferiores Z1< Z2< ... < ZL≤Z y detenerse cuando ZU−ZL≤ϵ, donde ϵes un peque˜no valor positivo que nosotros elegiremos. Definici´on 6. Sean ZUyZLuna cota superior y una cota inferior del problema (COP), respectivamente. Se define el Gap Relativo (como se define en [14]) entre dichas cotas como: gapUL =ZU−ZL ZU ×100 % Cada soluci´on factible x∗∈Xotorga una cota inferior o primal Z=c(x∗)≤Z. Por otra parte, para encontrar cotas superiores para un problema de maximizaci´on (o cotas inferiores para un problema de minimizaci´on), que llamaremos cotas duales, se reemplaza el problema original por un problema de optimizaci´on m´as simple cuyo valor ´optimo sea mayor o igual que el del problema original Z. Para que este problema, al que denotaremos problema relajado, cumpla esta propiedad, hay dos opciones: 1) Aumentar el conjunto de soluciones factibles. 2) Reemplazar la funci´on objetivo por una funci´on cuyo valor sea igual o superior en cada punto factible. Definici´on 7. Un problema (RP) ZRP =max{f(x) : x∈T⊆Rn}es una relajaci´on del problema (IP) Z=max{c(x) : x∈X⊆Rn}si: 1) X⊆T, y 6 Secci´on 1.3. Relajaciones 2) f(x)≥c(x)∀x∈X. Proposici´on 3. Si (RP) es una relajaci´on de (IP), entonces ZRP ≥Z. Demostraci´on. Si x∗es soluci´on ´optima de (IP), x∗∈X⊆TyZ=c(x∗)≤f(x∗). Como x∗∈T,f(x∗) es una cota inferior de ZRP yZ≤f(x∗)≤ZRP . 1.3. Relajaciones Empezamos esta secci´on recordando los conceptos definidos anteriormente de cota primal, que es una cota inferior proporcionada por una soluci´on factible del problema de maximizaci´on, y el concepto de cota dual, que es una cota superior del problema de maximizaci´on. En el caso de tener un problema de minimizaci´on, las cotas primales ser´an l´ogicamnte cotas superiores, y las cotas duales ser´an cotas inferiores. Definici´on 8. Dado un problema de programaci´on entera Z=max{cx :x∈P∩Zn}con formulaci´on P={x∈Rn +:Ax ≤b}, la relajaci´on de la programaci´on lineal es ZLP =max{cx :x∈P}. Proposici´on 4. Sean P1, P2formulaciones para el problema de programaci´on entera max{cx : x∈X⊆Zn}de forma que P1es mejor formulaci´on que P2, es decir, P1⊂P2. Si ZLP i=max{cx :x∈Pi}para i= 1,2son los valores asociados a cada relajaci´on lineal, entonces ZLP 1≤ZLP 2para cualquier c. Las relajaciones no solo proporcionan cotas duales, sino que a veces tambi´en nos permiten probar optimalidad. Proposici´on 5. Consideremos un problema (IP) y una relajaci´on de este problema (RP): 1) Si una relajaci´on (RP) es infactible, el problema original (IP) es infactible. Cap´ıtulo 1. Preliminares 7 2) Sea x∗soluci´on ´optima de (RP). Si x∗∈Xyf(x∗) = c(x∗), entonces x∗es soluci´on ´optima de (IP). Demostraci´on. 1) Al ser (RP) infactible, T=∅y por tanto X=∅. 2) Como x∗∈X,Z≥c(x∗) = f((x∗) = ZRP ). Como Z≤ZRP , obtenemos c(x∗) = Z=ZRP . 1.4. Cotas Primales Consideramos algunas formas simples de obtener soluciones factibles o cotas primales. Definici´on 9. Un problema (RE) ZRE =max{f(x) : x∈T⊆Rn}es una restricci´on del problema entero Z=max{c(x) : x∈X⊆Rn}si: T⊆X, y f(x)≤c(x)∀x∈X. Proposici´on 6. Si (RE) es una restricci´on de (IP) con soluci´on ´optima xRE,xRE es factible en (IP) y c(xRE)≤Z. Algunas formas de construir restricciones son fijar valores de algunas variables, reemplazar desigualdades por igualdades o a˜nadir restricciones adicionales. Esto se suele hacer de forma que la restricci´on resultante sea un problema m´as f´acil que se pueda resolver r´apidamente. 1.5. Algoritmo del simplex y de punto interior Consideremos el problema: max cx s.a Ax ≤b x≥0 8 Secci´on 1.6. Algoritmos de ramificaci´on y poda El m´etodo simplex es un algoritmo iterativo aplicado a problemas de programaci´on lineal. El m´etodo comienza a˜nadiendo nuevas variables, llamadas variables de holgura, de forma que al a˜nadirlas a las restricciones, est´as pasan a ser de igualdad. El algoritmo divide en las iteraciones las variables en b´asicas y no b´asicas, y va mejorando la soluci´on en el sentido de hacer su valor objetivo mayor hasta que no se puede continuar. Por otro lado, los algoritmos de punto interior buscan el punto ´optimo por el interior del dominio, mientras que el simplex busca por los v´ertices. Adem´as, el algoritmo de punto interior encuentra el ´optimo en tiempo polinomial, aunque en la pr´actica el algoritmo del simplex converge m´as r´apido. Para m´as informaci´on de estos algoritmos, ver [15]. 1.6. Algoritmos de ramificaci´on y poda Consideramos el problema: Z=max{cx :x∈S} Nos preguntamos c´omo dividir el problema en una serie de subproblemas que sean m´as f´aciles de resolver. Proposici´on 7. Sea S=S1∪ · · · ∪ Skuna descomposici´on de S en conjuntos m´as peque˜nos, y sea Zk=max{cx :x∈Sk}para k= 1,· · · , K. Entonces, Z=maxkZk. Proposici´on 8. Sea S=S1∪ · · · ∪ Skuna descomposici´on de S en conjuntos m´as peque˜nos, y sea Zk=max{cx :x∈Sk}para k= 1,· · · , K,¯ Zkuna cota superior de ZkyZkuna cota inferior de Zk. Entonces, ¯ Z=maxk¯ Zkes una cota superior de Z y Z=maxkZkes una cota superior de Z. Ahora veremos c´omo realizar un algoritmo de ramificaci´on y poda basado en la relajaci´on lineal. Se inicia el algoritmo con el conjunto Scon formulaci´on Pen una lista L. La cota inferior inicial es Z=−∞ yx∗inicialmente ser´a vac´ıo. Despu´es se usa un algoritmo heur´ıstico para encontrar una soluci´on factible xHcon valor objetivo ZH, y se actualiza Z=ZHyx∗=xH. A continuaci´on se realizar´a un paso de comprobaci´on que llamaremos inicio del bucle, y al que volveremos de forma recurrente cada vez que se realice una poda. En este paso se comprueba si la lista L es vac´ıa. Si lo es el algoritmo se detendr´a. En caso contrario, se elimina Sicon formulaci´on Pide L para un cierto i. Despu´es se resuelve la relajaci´on lineal de la formulaci´on Pi, con soluci´on xi(LP) y valor objetivo ¯ Zi. en caso de Pisea vac´ıo se poda o descarta la ramificaci´on por Sipor infactibilidad, y si ¯ Zi≤Zse poda por acotaci´on. Si no se poda, se Cap´ıtulo 1. Preliminares 9 comprueba si xi(LP) es entera o no. Si lo es, se actualiza la cota primal Z=¯ Zi, se hace x∗=xi(LP) y se poda por optimalidad. Recordemos que cada vez que se realiza una poda se vuelve al inicio del bucle. Si por el contrario xi(LP) no es de coordenadas enteras, se devuelven dos subproblemas Si 1ySi 2a la lista L con formulaciones Pi 1yPi 2y cotas superiores ¯ Zi, y se vuelve al inicio del bucle. A partir de la primera vez que se vuelve al bucle en la que contemos con una cota superior ¯ Zi, antes de resolver la relajaci´on lineal, se comprobar´a si ¯ Zi≤Z, y en tal caso se podar´a por acotaci´on. 1.7. Algoritmos de planos de corte Consideramos el problema: (IP)max{cx :x∈X} donde X={Ax ≤b, x ∈Zn +}. Definici´on 10. Una desigualdad πx ≤π0es una desigualdad v´alida para X⊆Rnsi πx ≤π0∀x∈X. Supongamos que X=P∩Zny que conocemos una familia Fde desigualdades v´alidas πx ≤π0, (π, π0)∈ P para X. En muchos casos, Fcontiene demasiadas desigualdades (2no m´as) como para ser agregadas a priori. Adem´as, dado una funci´on objetivo espec´ıfica, el objetivo no es encontrar la envolvente convexa completa, sino encontrar una buena aproximaci´on de conv(X). A continuaci´on se describe un algoritmo b´asico de planos de corte para (IP), max{cx :x∈X}, que genera desigualdades ´utiles a partir de F. El algoritmo de plano de corte se inicia estableciendo t=0yP0=P. En cada iteraci´on t, se resuelve el problema lineal Zt=max{cx :x∈Pt}. Sea xtla soluci´on ´optima. Si xtes de coordenadas enteras, se detiene el algoritmo y xtes ´optimo de (IP). En caso contrario, se resuelve el problema de separaci´on de xtyF. Si se encuentra una desigualdad (πt, πt 0)∈ F con πtxt> πt 0, se establece Pt+1 =Pt∩{x:πtx≤πt 0}y se aumenta t. En caso contrario se detiene. Si el algoritmo termina sin encontrar una soluci´on entera para IP, Pt=P∩ {x:πix≤πi 0, i = 1,· · · , t} es una formulaci´on mejorada que se puede dar como input a un solver que utilice un algoritmo de ramificaci´on y corte o branch-and-cut. En la pr´actica suele ser mejor a˜nadir varios cortes no 10 Secci´on 1.7. Algoritmos de planos de corte verificados en una sola iteraci´on que a˜nadir un solo corte. Cap´ıtulo 2 Dualidad Lagrangiana En este segundo cap´ıtulo profundizaremos en la Relajaci´on Lagrangiana citando el contenido de [17]. En concreto, veremos que esta clase de relajaciones busca eliminar restricciones que complican el problema pas´andolas a un t´ermino de penalizaci´on en la funci´on objetivo. Tambi´en se describir´a el algoritmo de subgradiente, que en cada iteraci´on nos proporcionar´a cotas que busquen la convergencia al valor ´optimo del problema original, y se analizar´a c´omo seleccionar los par´ametros del algoritmo para asegurar la convergencia. 2.1. Relajaci´on Lagrangiana Consideremos el problema IP: (IP)          Z= max cx s.a Ax ≤b Dx ≤d x∈Zn + Supongamos que las restricciones Dx ≤dson buenas en el sentido de que si el problema solo tuviera esas restricciones ser´ıa f´acil de resolver. Entonces, si eliminamos las restricciones complicantes o complicadas Ax ≤b, la relajaci´on resultante es m´as f´acil de resolver que el problema original (IP). Sin embargo, la cota resultante puede ser d´ebil porque algunas restricciones retiradas sean importantes y est´en siendo ignoradas. Una manera de abordar esta dificultad ser´a utilizar la Relajaci´on Lagrangiana. 11 12 Secci´on 2.1. Relajaci´on Lagrangiana Consideremos el problema IP de una forma algo m´as general: (IP)      Z= max cx s.a Ax ≤b x∈X donde Ax ≤bson mrestricciones complicadas. Para cualquier valor de u= (u1, ..., um)≥0, definimos el problema: (IP(u)) (Z(u) = max cx +u(b−Ax) x∈X Proposici´on 9. El problema (IP(u)) es una relajaci´on de (IP) ∀u≥0. Demostraci´on. Recordemos que (IP(u)) es una relajaci´on de (IP) ∀u≥0 si: La regi´on factible contiene a la regi´on factible original. Esto se cumple pues S={x: Ax ≤b, x ∈X} ⊆ X. El valor objetivo es como m´ınimo igual de grande en (IP(u)) que en (IP) para cada una de las soluciones factibles de (IP). Si x∈Syu≥0, entonces Ax ≤b, u(b−Ax)≥0 y entonces cx +u(b−Ax)≥cx ∀x∈S. En el problema (IP(u)) las restricciones complicadas se tratan a˜nadi´endolas a la funci´on objetivo con un t´ermino de penalizaci´on u(b−Ax). Se llama a umultiplicador de Lagrange asociado a las restricciones Ax ≤b. El problema (IP(u)) se llama Relajaci´on Lagrangiana de (IP) con par´ametro u. Al ser (IP(u)) una relajaci´on de (IP), Z(u)≥Zy tenemos as´ı una cota superior del valor ´optimo de (IP). Para encontrar la mejor cota superior, es decir, la m´as baja, de entre todos los posibles valores de u, hay que resolver el Problema Dual Lagrangiano: (LD) WLD = m´ın{Z(u) : u≥0} Veamos un ejemplo num´erico: Cap´ıtulo 2. Dualidad Lagrangiana 19 2.4. Algunas consideraciones para el algoritmo de subgradiente 2.4.1. Una forma de seleccionar las longitudes de paso Vamos a ver una forma de seleccionar las longitudes de paso que se detalla en [2], basada en el apartado c) del Teorema 2 para que la convergencia sea lo m´as r´apida posible. Sea ξkel subgradiente (b−Ax(uk)). Nos preguntamos ahora qu´e longitudes de paso µkdebemos escoger en cada iteraci´on. El valor ideal ser´a el que minimice ∥u∗−(uk+µξk)∥2, donde u∗es el uque optimiza WLD. Este valor ideal viene dado por µ∗ k=ξt k(u∗−uk) ∥ξk∥2 Adem´as, por la concavidad de Z(u), tenemos que µ∗ k=ξt k(u∗−uk) ∥ξk∥2≥Z∗−Z(uk) ∥ξk∥2 donde Z∗= sup{Z(u) : u≥0}. Una posible aproximaci´on ser´ıa µk=¯ Z−Z(uk) ∥ξk∥2 donde ¯ Z≥Z∗. Este valor lo podemos obtener con el valor objetivo de una soluci´on factible del problema de minimizaci´on original. Pero en lugar de escoger un valor fijo para ¯ Z, nos gustar´ıa actualizar este valor peri´odicamente de forma que ¯ Z−Z∗ ∥ξk∥2est´e en relaci´on directa al gap de µ∗ k=ξt k(u∗−uk) ∥ξk∥2≥Z∗−Z(uk) ∥ξk∥2. As´ı, ¯ Zvendr´a dado por ¯ Z=αrZ0+ (1 −αr)Zc, donde Z0es un n´umero fijo tal que Z0≥Z∗,Zces el mejor valor objetivo hasta la iteraci´on k, y{αr}es una sucesi´on mon´otona decreciente tal que α0= 1 y αr→ϵ0>0. Relacionaremos la sucesi´on {αr}con {uk}de la siguiente manera: si tras una cantidad de iteraciones el valor Zcno mejora, probablemente estemos tomando longitudes de paso demasiado 20 Secci´on 2.4. Algunas consideraciones para el algoritmo de subgradiente grandes, por lo que necesitamos reducir ¯ Z. Entonces, incrementamos el valor de r en 1 unidad, haciendo as´ı que αry¯ Zdisminuyan. Por otra parte, nos interesa evitar que ¯ Zest´e por debajo del valor de Z∗, cosa que puede ocurrir por ejemplo en las primeras iteraciones cuando Zcest´a demasiado por debajo de Z∗. Entonces debemos seleccionar {αr}de forma que otorgue m´as peso a Z0en las etapas iniciales, ya que el valor ´optimo Zcen esas etapas estar´a alejado del ´optimo Z∗, y que vaya dando, gradualmente al principio y mas r´apido despu´es, m´as valor a Zcconforme este valor se aproxima a Z∗. Una sucesi´on {αr}que optimiza estas exigencias es: αr=(e−0,6933(r/r1)3,26 si r < r2 ϵ0si r≥2 Para la elecci´on de los par´ametros r1yϵ0conviene probar varios valores y ver cual ofrece mejores resultados. Normalmente si Z0est´a pr´oximo a Z∗, lo deseable ser´a un valor alto para r1 y uno bajo para ϵ0, y al contrario, si Z0est´a alejado a Z∗interesar´a un r1bajo y un ϵ0alto. Por su parte, r2es el menor entero que satisface ϵ0≥e−0,6933(r2/r1)3,26 El procedimiento por el que se llega a esta sucesi´on {αr}se puede consultar en la referencia indicada al principio de esta secci´on. Con todo esto, lo primero que se realiza en el algoritmo es elegir unos multiplicadores de Lagrange iniciales u0. Se suelen seleccionar todos los multiplicadores nulos. En general, dado un uk, se obtiene la soluci´on x(uk) y el subgradiente ξk= (b−Ax(uk)) correspondiente. Acto seguido se calcula µkde la expresi´on µk=¯ Z−Z(uk) ∥ξk∥2y obtenemos los multiplicadores de la siguiente iteraci´on uk+1 mediante uk+1 i=max{uk i−µk(b−Ax(uk)i),0}para cada i= 1, ..., m. Por el bien de la convergencia, se llamar´a a esta iteraci´on un ´exito o un fracaso en funci´on de si Z(uk+1)−Zc≥ϵ > 0 se cumple o no, donde ϵes un valor de tolerancia previamente fijado. Si se experimentan una cantidad νde fracasos consecutivos, donde νtambi´en es un valor entero que fijamos al inicio del experimento, entonces incrementaremos el contador ren 1 unidad, actualizamos αry retornamos al multiplicador previo que otorgaba la mejor soluci´on hasta el momento. Esto evita que la sucesi´on {uk}se pierda durante demasiadas iteraciones sin mejorar. Continuamos as´ı hasta que rllegue a r2. En este momento, αrse hace igual a ϵ0y se congela ese valor, de forma que ¯ Z=ϵ0Z0+ (1 −ϵ)Zc, y se pasar´an a seleccionar las longitudes de paso como: µk=1 βk[¯ Z−Z(uk) ∥ξk∥2] Cap´ıtulo 2. Dualidad Lagrangiana 21 donde βk= 1 inicialmente, y se incrementa en 2 unidades cada νiteraciones sin importar si es un ´exito o un fracaso. A partir de esta segunda fase tras llegar rar2, el retorno a la mejor soluci´on se realiza mientras βksea menor que un entero ¯ β≥1. Una vez se supera este valor no se realiza el reseteo a la mejor soluci´on hasta el momento. 2.4.2. Condiciones de parada Antes de comenzar a iterar, se deben establecer una serie de condiciones para que el algoritmo cese. Algunas de las m´as utilizadas son: 1. Condici´on sobre las longitudes de paso: Si en alguna iteraci´on k la longitud µkes menor que un ϵpreviamente fijado se detendr´an las iteraciones. 2. Condici´on sobre el gap relativo: En cada iteraci´on, se tienen una cota superior, procedente de una soluci´on factible, y la mejor cota inferior conseguida hasta el momento. Si el Gap Relativo entre estas cotas ¯ ZU−ZL ¯ ZU×100 %, donde ¯ ZUes la mejor cota superior hasta el momento y ZLes la mejor cota inferior hasta el momento, es menor que un valor determinado que fijamos previamente, se detendr´a el algoritmo. 3. Si el algoritmo realiza un n´umero de iteraciones que hemos elegido previamente como el n´umero m´aximo de iteraciones que se permitir´an, entonces parar´a. 4. Se puede establecer un tiempo l´ımite tras el que se detenga el proceso iterativo. 2.4.3. Divisi´on en subproblemas enteros o continuos En el paso del algoritmo en el que se resuelve el Problema Lagrangiano IP(uk) se pueden realizar modificaciones. Con el fin de hacer el algoritmo m´as r´apido, se puede resolver la relajaci´on lineal IP(uk), que otogar´a una cota inferior menor, es decir, menos buena, que la que proporcionar´ıa IP(uk) con las variables enteras, pero la calcular´a m´as r´apido. Las opciones de combinar la resoluci´on de los problemas IP(uk) con variables continuas o enteras son muy diversas, se puede hacer de la forma que deseemos. Por ejemplo, se puede realizar un n´umero fijado de iteraciones resolviendo la relajaci´on lineal, y realizar las posteriores iteraciones con variables enteras. Tambi´en se puede iterar con la relajaci´on lineal hasta que el Gap de Optimalidad sea menor que una valor fijado, y a partir de ah´ı usar variables enteras. 22 Secci´on 2.4. Algunas consideraciones para el algoritmo de subgradiente Cap´ıtulo 3 Problemas de Asignaci´on Multiperiodo En este cap´ıtulo se definir´an los problemas de asignaci´on multiperiodo, y se desarrollar´an los problemas Axial y Planar de 3´ındices, que son los problemas a los que se le aplicar´a la Relajaci´on Lagrangiana y el algoritmo de subgradiente vistos en el Cap´ıtulo 2. Por ´ultimo se definir´a el Dual Lagrangiano de estos problemas y su descomposici´on en subproblemas m´as sencillos. 3.1. Permutaciones y el problema de asignaci´on Los problemas de asignaci´on consisten en asignar nelementos, como pueden ser por ejemplo agentes, a otros nelementos, como por ejemplo tareas. Una forma en matem´aticas de describir las asignaciones es mediante las permutaciones. Definici´on 13. Consideremos un conjunto Xde nelementos. Una permutaci´on φes una aplicaci´on biyectiva que asigna a cada ´ındice i= 1, ..., n un elemento del conjunto X. Se puede representar como 1 2 ··· n φ(1) φ(2) ··· φ(n). Tal representaci´on significa que ies asignado a φ(i) para i= 1, ..., n. Cada permutaci´on φdel conjunto {1,2,· · · , n}corresponde de forma ´unica a una matriz n×n de permutaci´on Xφ= (xij) donde: xij =(1 si j=φ(i) 0 en caso contrario. 23 24 Secci´on 3.1. Permutaciones y el problema de asignaci´on Veamos un ejemplo de problema de asignaci´on. Supongamos que nagentes deben ser asignados antareas, de forma que cada agente solo puede ser asignado a una tarea, y cada tarea ha de ser realizada por un ´unico agente. Consideremos una matriz de costes n×n,C= (cij), donde cij es el coste de asignar el agente i a la tarea j, y supongamos que queremos minimizar el coste total de las asignaciones, es decir, queremos encontrar la permutaci´on que encuentre m´ın φ∈Sn n X i=1 ciφ(i). Siendo X= (xij) una matriz de permutaci´on, se puede formular el problema de asignaci´on como m´ın n X i=1 n X j=1 cijxij s.t. n X j=1 xij = 1 i= 1,2, ..., n n X i=1 xij = 1 j= 1,2, ..., n xij ∈ {0,1}i, j = 1,2, ..., n. Un problema de asignaci´on multiperiodo es un problema consistente en asignar nagentes a n tareas en ninstantes de tiempo. Sea cijk el coste de asignar el agente ia la tarea jen el instante k. Buscamos encontrar una asignaci´on φde los agentes a las tareas y una asignaci´on ψde los agentes a los instantes de tiempo. Esto nos lleva a Problema de Asignaci´on Axial de 3 ´ındices: m´ın φ,ψ∈Sn n X i=1 ciφ(i)ψ(i), donde Sndenota el conjunto de todas las permutaciones de los enteros {1,2, ..., n}. Este Problema Axial, junto a otro problema multiperiodo conocido como el Problema Planar, se formular´an en las posteriores secciones y ser´an los modelos que se implementar´an en el estudio computacional. Sin embargo, estos nos son los ´unicos problemas de asignaci´on multiperiodo que existen. Esta clase de problemas es muy diversa y tiene m´ultiples aplicaciones pr´acticas. Veamos a continuaci´on algunos modelos para ciertos problemas y la motivaci´on de su planteamiento y resoluci´on. Cap´ıtulo 3. Problemas de Asignaci´on Multiperiodo 25 3.1.1. Aplicaciones de problemas de asignaci´on Los problemas de asignaci´on tienen diversas aplicaciones en diferentes campos. En primer lugar se plantea el problema de horarios en la universidad. En este problema se busca compatibilizar los horarios de alumnos, profesores y dem´as personal de acuerdo a una serie de normas o restricciones. Algunas de ellas son: No se permiten colisiones. Una colisi´on ocurre cuando dos o m´as cursos se programan a la vez para el mismo profesor, en el mismo aula o para el mismo grupo de alumnos. Tambi´en cuando dos o m´as profesores se asignan al mismo grupo de alumnos para impartir la misma asignatura; y por ´ultimo, cuando dos o m´as aulas se asignan al mismo grupo de alumnos y al mismo curso. Un horario debe estar completo. Un horario se considera completo cuando aparecen todos los cursos planificados para cada grupo de estudiantes en el horario, con la cantidad correcta de per´ıodos de tiempo para cada curso y cada parte de cada curso. Adem´as, cuando a cada miembro del personal docente se le asigna el n´umero total de per´ıodos de ense˜nanza que requiere la instituci´on. Un horario debe poder adaptarse a solicitudes de sesiones de per´ıodos de ense˜nanza consecutivos. Dependiendo del curso y el n´umero de per´ıodos de ense˜nanza asignados por semana, un profesor puede optar por impartir el curso en sesiones de un solo per´ıodo o de varios per´ıodos. El horario debe ser capaz de programar un curso determinado en cualquier esquema que el profesor responsable elija seguir. Un horario debe poder adaptarse a solicitudes de sesiones repetidas de un curso dado o parte de un curso. Esta regla se refiere a la necesidad presentada por algunos cursos de sesiones de repaso o trabajos de laboratorio dise˜nados para grupos peque˜nos. En este caso, la misma sesi´on debe repetirse varias veces para acomodar el n´umero total de estudiantes registrados. El problema se puede consultar en detalle en [5], art´ıculo en el que se desarrolla la formulaci´on con m´as variables y restricciones que las aqu´ı presentadas. Otro problema es el de realizar una planificaci´on eficiente de actividades culturales y eventos, teniendo en cuenta diferentes partes involucradas como ayuntamientos, empresas y administraci´on, entre otras cosas. Este problema se aborda en [12], motivado por el deseo de los municipios de brindar una oferta cultural atractiva con respecto a una combinaci´on de preferencias, principalmente relacionadas con la elecci´on de empresas y del calendario. Por su 26 Secci´on 3.2. Problema de asignaci´on Planar de 3 ´ındices parte, las compa˜n´ıas art´ısticas tratan de maximizar sus beneficios actuando el mayor tiempo posible en los municipios m´as atractivos. La administraci´on tiene como objetivo maximizar el beneficio social. As´ı, se propone medir la calidad de las soluciones relacionando esos factores. Algunas restricciones que se modelan en el art´ıculo son: Una misma compa˜n´ıa no puede actuar el mismo d´ıa. Ning´un agente puede actuar por j-´esima ocasi´on en un lugar m´as de 1 vez. Se proh´ıbe que dos actuaciones coincidan el mismo d´ıa y en el mismo lugar. Ning´un sitio puede exceder el presupuesto permitido para una modalidad, y se otorga un descuento en el coste por la j-´esima repetici´on de un agente en el sitio. Tambi´en se puede encontrar aplicaciones dirigidas a servicios de cuidado de hogares, como se aborda en [13]. El problema consiste en asignar personal m´edico a hogares de clientes para llevar a cabo determinados cuidados. En este problema se establecen restricciones como la obligaci´on de que todos los pacientes sean atendidos alg´un d´ıa, establecer la dependencia de la realizaci´on de unos cuidados de otros cuidados previos, o la prohibici´on de que un equipo vuelva al hospital y luego salga de este el mismo d´ıa. Otras aplicaciones est´an relacionadas con los modelos Planar y Axial que se definir´an en las siguientes secciones. En [4] se recogen las aplicaciones que mencionaremos a continuaci´on, y recoge referencias a estudios de estas aplicaciones. El modelo Planar tiene aplicaciones directas en modelos de horarios y en la asignaci´on de intervalos de tiempo para sistemas de telecomunicaciones, que usan sat´elites que deben comunicarse con la Tierra. El modelo Axial tiene aplicaciones en diversos campos, siendo ´util por ejemplo para invertir capital en diferentes localizaciones durante una cantidad de tiempo, para modelar el ensamblaje de placas de circuitos o para la laminaci´on de lingotes. 3.2. Problema de asignaci´on Planar de 3 ´ındices Definimos el Problema de Asignaci´on Planar de 3 ´ındices (Planar 3AP) en lo que sigue. Se dice que n permutaciones φ1, φ2, ..., φnson mutuamente distintas si φr(i)=φs(i) para cualquier i= 1,2, ..., n yr=s. Dados n3coeficientes de costes cijk (i, j, k = 1,2, ...n), el problema consiste en encontrar n permutaciones mutuamente distintas tal que: n X k=1 n X i=1 ciφk(i)k Cap´ıtulo 3. Problemas de Asignaci´on Multiperiodo 27 sea m´ınimo. Su correspondiente programaci´on lineal entera es: m´ın n X i=1 n X j=1 n X k=1 cijkxijk (3.1) s.t. n X k=1 xijk = 1 i, j = 1,2, ..., n (3.2) n X j=1 xijk = 1 i, k = 1,2, ..., n (3.3) n X i=1 xijk = 1 j, k = 1,2, ..., n (3.4) xijk ∈ {0,1}i, j, k = 1,2, ..., n (3.5) La variable xijk tomar´a el valor 1 si se asigna el agente ia la tarea jen el instante k. De esta forma, la funci´on objetivo (3.1) nos indica que hay que minimizar los costes de las combinaciones que hemos escogido en las asignaciones. La restricci´on (3.2) nos indica que cada agente debe realizar cada tarea exactamente en un ´unico instante de tiempo. De forma an´aloga, (3.3) nos fuerza a que cada agente debe estar realizando una y solo una tarea en cada instante de tiempo. Por ´ultimo, (3.4) nos dice que cada tarea se debe estar llevando a cabo en cada instante de tiempo por un ´unico agente. El Planar 3AP est´a muy relacionado con los cuadrados Latinos. Un cuadrado Latino es un conjunto de n×nentradas lij que toman los valores de 1 a n. Cada fila y cada columna de un cuadrado Latino contiene exactamente una entrada con el valor k, con 1 ≤k≤n. Cada soluci´on factible del Planar 3AP se puede representar como un cuadrado Latino. Sea L= (lij) un cuadrado Latino de tama˜no n. Entonces, para i, j = 1,2, ..., n,lij es el ´unico valor de ´ındice k tal que xijk = 1 en una soluci´on factible de Planar 3AP. Un ejemplo de cuadrado latino ser´ıa:    1 2 3 2 3 1 3 1 2    Este cuadrado latino corresponde a la asignaci´on: agente 1, tarea 1 e instante de tiempo 1. agente 1, tarea 2 e instante de tiempo 2. 28 Secci´on 3.2. Problema de asignaci´on Planar de 3 ´ındices agente 1, tarea 3 e instante de tiempo 3. agente 2, tarea 1 e instante de tiempo 2. agente 2, tarea 2 e instante de tiempo 3. agente 2, tarea 3 e instante de tiempo 1. agente 3, tarea 1 e instante de tiempo 3. agente 3, tarea 2 e instante de tiempo 1. agente 3, tarea 3 e instante de tiempo 2. En concreto, este cuadrado latino es un cuadrado latino reducido, es decir, un cuadrado latino que tiene la permutaci´on identidad en la primera fila y en la primera columna. El problema Planar 3AP tiene muchas aplicaciones a la hora de organizar y repartir tareas en un horario para distintos empleados, pero es un problema NP-duro complicado de resolver (ver [17]). Y es que el n´umero de soluciones factibles, es decir, de cuadrados latinos existentes, aumenta muy r´apido en funci´on de n. Por ejemplo, si n= 9, el n´umero de soluciones factibles ser´ıa aproximadamente de 551026 de acuerdo con Bammel y Rothstein en [1]. En [9] se muestra el n´umero de cuadrados latinos reducidos hasta orden 11, y se afirma que el n´umero de cuadrados latinos hasta orden n es igual al n´umero de cuadrados latinos reducidos hasta ese orden multiplicado por n!(n−1)!. Se comprueba en la siguiente tabla que el n´umero de cuadrados latinos reducidos (y por tanto tambi´en los totales) aumenta a una velocidad muy grande conforme aumenta n. Por este motivo usaremos la relajaci´on lagrangiana, que puede ayudarnos a obtener soluciones cercanas a la ´optima en un tiempo aceptable. Tabla 3.1: N´umero de cuadrados latinos reducidos de orden n ncuadrados latinos reducidos 1 1 2 1 3 1 4 4 5 56 6 9408 7 16942080 8 535281401856 9 377597570964258816 10 7580721483160132811489280 11 5363937773277371298119673540771840 No se conocen muchos algoritmos para resolver el Planar 3AP. El primer algoritmo de ramificaci´on y acotaci´on se remonta a [16], quien calcul´o cotas inferiores mediante reducci´on Cap´ıtulo 4. Estudio computacional 35 Algoritmo 2: Algoritmo Greedy para el Axial 3AP input : Se solicita una matriz de costes n×n×n,C= (cijk)n×n×n output: Una cota superior o primal ZUB 1Se fija un valor Mtal que M > cijk ∀i, j, k = 1,· · · , n 2for cada agente i= 1,· · · , n do 3Asigna el agente ia la tarea jy el tiempo kque cumplen que cijk es el menor valor entre los cij∗k∗∀j∗, k∗= 1,· · · , n 4for cada agente i∗= 1,· · · , n do 5for cada tiempo k∗= 1,· · · , n do 6ci∗jk∗=M 7for cada agente i∗= 1,· · · , n do 8for cada tarea j∗= 1,· · · , n do 9ci∗j∗k=M Algoritmo 3: Algoritmo para el Planar 3AP input : Se solicita una matriz de costes n×n×n,C= (cijk)n×n×n output: Una cota superior o primal ZUB 1contador=0 2for cada agente i= 1,· · · , n do 3for cada tiempo k= 1,· · · , n do 4if k+contador < n then 5Asigna la tarea ie instante ka la tarea k+contador 6if k+contador ≥nthen 7Asigna la tarea ie instante ka la tarea k+contador −n 8contador =contador + 1 36 Secci´on 4.4. Cotas primales y duales del algoritmo exacto En el caso del Planar 3AP, no resulta evidente obtener una cota primal. Una opci´on que implementamos en Python es la descrita en el Algoritmo 3, que selecciona un cuadrado latino como soluci´on factible y eval´ua su valor objetivo. El cuadrado latino que selecciona el algoritmo es el siguiente: 1. Para el agente 1: asigna la tarea 1 al tiempo 1, la tarea 2 al tiempo 2, ... 2. Para el agente 2: asigna la tarea 2 al tiempo 1, la tarea 3 al tiempo 2, ... 3. Para el agente 3: asigna la tarea 3 al tiempo 1, la tarea 4 al tiempo 2, ... Existen otros m´etodos m´as complicados para el Planar 3AP, como el descrito en [8]. En este art´ıculo se expone un algoritmo de ramificaci´on y acotaci´on para el Planar 3AP, de forma que, en lugar de ramificar una sola variable, ramifica un conjunto de variables que se encuentran en una misma restricci´on. El procedimiento para encontrar una cota superior encuentra primero una soluci´on inicial y posteriormente la mejora, mientras que para las cotas inferiores combina heur´ısticas duales y la relajaci´on lagrangiana. En particular, en el procedimiento de mejora de la cota superior, dada una soluci´on factible, es decir, un cuadrado latino, se le aplica una mejora local que encuentra una soluci´on vecina con menor valor objetivo. Se representa la soluci´on factible como un cuadrado latino, es decir, en forma de tabla n×nen el que en la entrada (i, j) se encuentra el ´ındice kal que se asignan el agente iy la tarea j, y se cambia el valor de alguna entrada. Eso convierte la soluci´on en infactible. Entonces, se desplaza desde la entrada modificada mediante filas y columnas hasta que se encuentra el ´ındice kque est´a repetido, y se cambia por el que fue modificado. Al acabar, se tiene un cuadrado latino diferente al original, y se eval´ua su valor objetivo para ver si es menor que el del cuadrado latino original. 4.4. Cotas primales y duales del algoritmo exacto Para cada problema de asignaci´on multiperiodo, se definieron funciones PLANAR yAXIAL que, mediante el solver Gurobi, y al darle como datos de entrada la matriz de costes (cijk)n×n×n, resuelve los problemas Planar 3AP y Axial 3AP con un algoritmo exacto de ramificaci´on y corte. Adem´as del valor ´optimo de la relajaci´on lineal del problema, Gurobi obtiene cotas superiores e inferiores que cuando coinciden nos garantizan haber encontrado el valor ´optimo del problema. A estas funciones se les puede poner un tiempo l´ımite, de forma que si llegado ese tiempo todav´ıa no han coincidido las cotas superior e inferior, la funci´on detendr´a los c´alculos y devolver´a el valor de estas cotas inferior y superior que ten´ıa en ese instante. Cap´ıtulo 4. Estudio computacional 37 4.5. Par´ametros iniciales del algoritmo de subgradiente Ahora que ya tenemos las herramientas para hallar el ´optimo, la relajaci´on lineal y cotas primales y duales de los problemas Planar 3AP y Axial 3AP, pasamos a implementar el algoritmo de subgradiente. Debemos fijar una serie de par´ametros ya descritos en el Cap´ıtulo 2 antes de comenzar a iterar. Estos son: Multiplicadores de Lagrange iniciales u0. Los elegiremos nulos de inicio, u0= 0. Max Iter: El n´umero de iteraciones m´aximo que se realizar´an hasta detener el algoritmo. Scale: es un par´ametro proporcional a la longitud de paso en cada iteraci´on. Desempe˜na un papel similar al de ηken el apartado c) del Teorema 2 del Cap´ıtulo 2, o al de 1 βk descrito en la secci´on 2.3. Su valor debe estar entre 0 y 2. SameLimit: n´umero m´aximo de iteraciones seguidas en las que la cota inferior no se actualiza. De llegar a esta cantidad de iteraciones seguidas en las que ocurra esto, se dividir´a el Scale entre 2 y el SameLimit se resetear´a a 0. Tiene el papel de νde la secci´on 2.3. ϵ: valor positivo, de forma que si el Optimality Gap es menor que esta cantidad, entonces se detendr´a el algoritmo. En nuestro caso elegiremos ϵ= 10−4. ϵstep: valor positivo, de forma que si la longitud de paso calculada en alguna iteraci´on es inferior a este valor, se detendr´a el algoritmo. 4.6. Algoritmo de subgradiente A continuaci´on se explica c´omo se aplicar´a el algoritmo de subgradiente visto en el Cap´ıtulo 2. PASO 0: Fijamos los par´ametros iniciales tal y como se establecen en la Secci´on 4.5. PASO 1: Nos situamos en la iteraci´on k. Contamos con los multiplicadores uk, con los que se resolver´a la relajaci´on lagrangiana con la funci´on lower boundP olower boundA, seg´un se trate el modelo Planar o Axial, y se obtendr´a una cota dual o inferior. Las funciones lower boundP ylower boundA tienen como dato de entrada los multiplicadores de Lagrange u, junto a la matriz de costes (cijk), y devuelve la soluci´on de la relajaci´on lagrangiana de los problemas 38 Secci´on 4.6. Algoritmo de subgradiente Planar 3AP y Axial 3AP, as´ı como su valor objetivo. Este supone una cota dual o inferior. Antes de empezar las iteraciones tambi´en es necesaria una cota inferior inicial. En nuestro experimento se escoger´a una cota inferior trivial, que es el valor 0. En otros casos tambi´en servir´ıa el valor objetivo ´optimo de la relajaci´on lineal del problema original aunque en algunas ocasiones, como en los problemas que nos ocupan, la fuerza del dual lagrangiano no permite que el valor de la relajaci´on lineal sea superado. Adem´as, este trabajo tambi´en pretende comprobar la viabilidad del algoritmo de subgradiente como alternativa al c´alculo de la relajaci´on lineal mediante el solver Gurobi. En caso de que la cota inferior obtenida de la relajaci´on lagrangiana sea mayor que la cota inferior que se pose´ıa al iniciar la iteraci´on, se actualizar´a como la mejor cota inferior hasta el momento. Si la cota inferior no se actualiza, se sumar´a una unidad al valor same, que es una variable auxiliar que definimos antes de comenzar a iterar y tiene como valor inicial 0. Esta variable servir´a para llevar la cuenta del n´umero de veces que no se actualiza la cota inferior. Si el valor de same iguala a SameLimit, entonces, tal y como explicamos antes, se dividir´a el valor de scale entre 2 para calcular la longitud de paso y se resetear´a la cuenta de same a 0. PASO 2: Con la soluci´on x(uk) del problema relajado calcularemos los slacks o subgradientes ξkcomo: b−Ax(uk) en el caso general. Si el subgradiente es nulo, entonces significa que la soluci´on x(uk) de la relajaci´on lagrangiana en esta iteraci´on cumple las restricciones relajadas y por lo tanto es factible en el problema original, por lo que su valor objetivo en el problema original se evaluar´a con la funci´on que definimos para encontrar cotas primales: upper boundP oupper boundA. Estas funciones tienen como argumento la matriz de costes (cijk) y una soluci´on factible x, y devuelven el valor objetivo en el problema original de la soluci´on factible x. Estas funciones sirven para evaluar el valor objetivo de una soluci´on factible en el caso de que el algoritmo de subgradiente encuentre alguna en alguna iteraci´on. Si la cota superior de la soluci´on factible es mejor, i.e., inferior a la que ten´ıamos hasta el momento, la sustituir´a como mejor cota superior hasta el momento. Si por el contrario el subgradiente no es nulo, la longitud de paso µkse calcular´a como µk=scale ∗ZUB −ZD ∥ξk∥2 donde ZUB es la cota superior y ZDes el valor objetivo de la relajaci´on lagrangiana con multiplicadores ukresuelto anteriormente. Cap´ıtulo 4. Estudio computacional 39 PASO 3: Con las longitudes de paso ya calculadas, que son las descritas en el apartado c) del teorema 2 de la secci´on 2.3, se obtienen los multiplicadores de Lagrange de la siguiente iteraci´on: uk+1 =uk+µkξk Estos multiplicadores son irrestrictos en signo debido a que las restricciones relajadas tanto en el Planar 3AP como en el Axial 3AP son de igualdad. Con estos multiplicadores uk+1 se dar´a paso al Paso 1 de la siguiente iteraci´on, siempre que no se cumpla alguna condici´on de parada: se ver´a si el paso µkes menor que 10−8y si el Optimality Gap es menor que el ϵ= 10−4que fijamos al principio, y de cumplirse alguna de estas condiciones, se detendr´a el algoritmo. En caso contrario, se continuar´a iterando. El proceso tambi´en terminar´a si se alcanza el n´umero m´aximo de iteraciones fijado. 4.7. Algoritmo de reparaci´on En esta secci´on se explica c´omo se ha implementado un algoritmo que toma la mejor soluci´on que ha encontrado el algoritmo de subgradiente aplicado al modelo Axial 3AP y la transforma en una soluci´on factible lo m´as parecida posible. PASO 1: El algoritmo empieza con una soluci´on inicial con todas las entradas nulas. Para cada agente i, obtiene el valor m´as grande de la soluci´on iterada en ese agente y lo llama xmax. Este valor generalmente ser´a 1. PASO 2: Para cada tarea y cada instante de tiempo, compara los costes de las entradas en las que la soluci´on iterada vale xmax, y selecciona y guarda los ´ındices de la tarea y el instante de tiempo (j, k) de menor coste, en caso de que se pueda. Esto es debido a que puede que no se seleccione ninguna tarea ni instante de tiempo por el siguiente motivo: en el caso en el que no se seleccione ning´un (j, k), significar´a que ninguna de las entradas con valor xmax en la soluci´on iterada se puede seleccionar por factibilidad, pues han sido bloqueados por los anteriores ´ındices seleccionados en iteraciones previas. En este caso, se tendr´an que seleccionar los ´ındices de la entrada (j, k) de menor coste en el agente i, aunque tenga un valor 0 en la soluci´on iterada. 40 Secci´on 4.8. Variantes de resoluci´on PASO 3: Una vez tenemos el par de ´ındices (j, k) para el agente i, se establece el valor de la soluci´on reparada xijk = 1, se bloquean todas las entradas de los agentes restantes correspondientes a la tarea jy el instante de tiempo k, y se pasa al siguiente agente. Al finalizar con la asignaci´on del ´ultimo agente, el algoritmo ha elaborado una soluci´on parecida a la iterada que adem´as es factible. En la implementaci´on, el algoritmo va recorriendo los agentes por orden, de modo que al haber menos entradas bloqueadas en los primeros, la soluci´on se parecer´a m´as a la iterada en los primeros agentes y en los ´ultimos se parecer´a menos. Se podr´ıa recorrer los agentes en otro orden en funci´on de si se quisiera priorizar m´as que se parezca en seg´un qu´e agentes. Otra posible variaci´on del algoritmo ser´ıa recorrer las tareas o los instantes de tiempo e ir seleccionando y bloqueando los otros ´ındices. 4.8. Variantes de resoluci´on Para resolver la relajaci´on lagrangiana, no solo se pueden usar las funciones lower boundP o lower boundA. Tambi´en se programaron las funciones lower boundPk olower boundAk que resuelven los ksubproblemas en los que se puede dividir la relajaci´on lagrangiana y devuelven el valor objetivo del subproblema. Una vez resueltos los subproblemas, se agregan sus valores para obtener el mismo valor que se obten´ıa usando las anteriores funciones pero de forma m´as eficiente. Tambi´en se puede alternar la resoluci´on de los problemas con las variables enteras o continuas como se explic´o en la Secci´on 2.4.4. 4.9. Resultados Los resultados que se detallan a continuaci´on se han obtenido un ordenador con un procesador Intel Core i7 de 2GHz y 8GB de memoria RAM. La versi´on de Gurobi instalada es la 10.0.1. 4.9.1. Resoluci´on con algoritmo exacto Vamos a ver los resultados obtenidos al resolver diversas instancias de los modelos Axial 3AP y Planar 3AP con Gurobi. Veremos para cada modelo una tabla cuyas columnas corresponden a los siguientes valores: Cap´ıtulo 4. Estudio computacional 41 n: el tama˜no de la instancia. g¯ UR: el gap entre la cota superior que encuentra Gurobi y el valor de la relajaci´on lineal. gUL: el gap entre la cota superior y la cota inferior que encuentra Gurobi. Tiempo: el tiempo en segundos que emplea Gurobi hasta encontrar la soluci´on ´optima. Est´a limitado a 300 segundos para el Axial y 60 segundos en el Planar, de modo que si llega a este tiempo, acabar´a la iteraci´on que est´e realizando y devolver´a las cotas superior e inferior que ten´ıa en ese momento. TiempoRL: el tiempo en segundos empleado en calcular la relajaci´on lineal. Est´a limitado a 120 segundos en el Axial y a 60 segundos en el Planar. z∗: ´optimo calculado por Gurobi. zRL: valor de la relajaci´on lineal. n g ¯ UR gUL Tiempo TiempoRL z∗zRL 5 0 0 0.01 0.01 44.62 44.62 10 0 0 0.08 0.11 19.29 19.29 15 14.2 0 0.64 0.22 35.64 30.58 20 5 0 1.05 0.46 36.39 34.57 55 2.48 1.8 ≥310.31 10.06 - 59.76 125 - - - ≥120 - - Tabla 4.1: Resultados del algoritmo exacto del Axial 3AP n g ¯ UR gUL Tiempo TiempoRL z∗zRL 5 0 0 0.02 0.14 712.66 712.66 6 2.88 0 0.33 0.09 946.2 918.98 8 2.64 0 0.57 0.05 1552.17 1511.04 11 4.15 0 17.44 0.11 2429.74 2328.89 15 12.63 12.05 ≥60 0.44 - 3421.4 60 - - - ≥60 - - Tabla 4.2: Resultados del algoritmo exacto del Planar 3AP En ambos modelos se observan buenos resultados de tiempo y gap peque˜nos en las instancias de menor tama˜no n, y a medida que este tama˜no aumenta, se observa que el tiempo y el gap tambi´en lo hacen. Para n= 55 en el Axial 3AP y para n= 15 en el Planar 3AP, se observa que el tiempo empleado para obtener el ´optimo supera a los tiempos que hab´ıamos fijado previamente, pero sigue siendo capaz de resolver la relajaci´on lineal. Adem´as, en estas instancias los gap son de 2.48 en el Axial y de 12.63 en el Planar, por lo que Gurobi sigue dando una buena acotaci´on 42 Secci´on 4.9. Resultados Figura 4.1: Modelo Axial 3AP. En azul, tiempo del ´optimo. En rojo, tiempo de la relajaci´on lineal. Figura 4.2: Modelo Planar 3AP. En azul, tiempo del ´optimo. En rojo, tiempo de la relajaci´on lineal. con gaps bajos. En cambio, para n= 125 en el Axial 3AP y para n= 60 en el Planar 3AP, los tiempos de la relajaci´on lineal que establecimos como m´aximos se superan, y el problema con variables enteras tarda demasiado en terminar su iteraci´on, por lo que no se recogen las cotas o gaps. En las siguientes gr´aficas se puede ver la evoluci´on de los tiempos de resoluci´on en funci´on del tama˜no n. Cap´ıtulo 4. Estudio computacional 43 4.9.2. Resoluci´on con algoritmo de subgradiente A continuaci´on se exponen los resultados obtenidos al aplicar el algoritmo de subgradiente a las instancias del Axial 3AP. Para el modelo Planar 3AP resulta muy complejo encontrar un algoritmo de reparaci´on de soluci´on factible, y la cota superior implementada que corresponde al cuadrado latino explicada en 4.3 no es eficiente. Esto es debido a que encontrar cuadrados latinos no es nada trivial. Por eso la implementaci´on del algoritmo de subgradiente del Planar 3AP en nuestro caso no proporciona resultados interesantes y no se ver´a. El algoritmo de subgradiente ha sido aplicado dividiendo la relajaci´on lagrangiana de la forma mostrada en la Secci´on 3.4. Veamos qu´e valores recoge la tabla de resultados: n: el tama˜no de la instancia. ZUB: la mejor cota superior obtenida tras aplicar el algoritmo de subgradiente, es decir, la menor cota superior de las obtenidas mediante el algoritmo voraz o mediante el algoritmo de reparaci´on de la soluci´on iterada. ZLB: la mejor cota inferior obtenida tras aplicar el algoritmo de subgradiente, es decir, la cota inferior m´as alta iterada. gUR: el gap entre la mejor cota superior ZUB y el valor de la relajaci´on lineal, que es la mejor cota inferior de todas las instancias. gLR: el gap entre la cota inferior ZLB y la cota inferior que encuentra la relajaci´on lineal. Se calcula para comprobar c´omo de cerca se queda la cota iterada de la cota de la relajaci´on lineal Tiempo: el tiempo en segundos que emplea el algoritmo de subgradiente hasta que se cumpla una condici´on de parada. Max Iter: n´umero m´aximo de iteraciones permitidas. Iteraciones: la cantidad de iteraciones realizadas hasta la finalizaci´on del algoritmo. ϵstep: valor del par´ametro. 44 Secci´on 4.9. Resultados n ZUB ZLB gUR gLR Tiempo max iter Iteraciones ϵstep 5 44.62 44.62 0 0 0.56 500 64 10−8 10 19.29 19.29 0 0 1.38 500 20 10−8 15 126.97 30.5 75.91 0.26 43.7 500 217 10−8 20 83.84 34.46 58.77 0.32 176.86 500 384 10−8 55 170.97 56.82 65.05 4.92 812.8 500 87 10−8 125 218.59 125 - - 108.56 1 1 10−8 Tabla 4.3: Resultados del algoritmo de subgradiente del Axial 3AP Figura 4.3: Modelo Axial 3AP. En azul, tiempo del ´optimo. En rojo, tiempo de la relajaci´on lineal. En negro, tiempo del algoritmo de subgradiente. Nos fijamos en primer lugar en gUR, el gap entre la mejor cota superior ZUB y el valor de la relajaci´on lineal. Para instancias de menor tama˜no, el algoritmo encontr´o el ´optimo, por lo que el gap fue nulo. Con instancias mayores, el gap aument´o a valores grandes, superando el 50 %. Sin embargo, en el caso de la instancia de tama˜no n= 125 se opt´o por realizar una sola iteraci´on, dado que el gran tama˜no de la instancia provocaba que cada iteraci´on requiriese mucho tiempo. Para esta instancia, la relajaci´on lineal no pudo ser obtenida en un tiempo menor a 120 segundos, pero gracias a la iteraci´on del algoritmo, en 108.56 segundos se pudo obtener una cota inferior. No solo eso, sino que el gap entre la mejor cota superior del algoritmo de reparaci´on y la cota inferior obtenidas es de un 41.89 %, que es un gap mejor (ya que es menor) que los gUR obtenidos en las instancias medianas de tama˜no 15, 20 o 55. ´ Indice de figuras 1.1. Formulaci´on P1paraelconjuntoS.......................... 4 1.2. Formulaci´on P2paraelconjuntoS.......................... 4 4.1. Modelo Axial 3AP. En azul, tiempo del ´optimo. En rojo, tiempo de la relajaci´on lineal. ......................................... 42 4.2. Modelo Planar 3AP. En azul, tiempo del ´optimo. En rojo, tiempo de la relajaci´on lineal. ......................................... 42 4.3. Modelo Axial 3AP. En azul, tiempo del ´optimo. En rojo, tiempo de la relajaci´on lineal. En negro, tiempo del algoritmo de subgradiente. . . . . . . . . . . . . . . 44 4.4. Modelo Axial 3AP. En negro, tiempo del algoritmo con ϵstep = 10−8. En azul, tiempo del algoritmo con ϵstep = 10−4. En rojo, tiempo del algoritmo con ϵstep = 10−2....................................... 46 51 52 ´ Indice de tablas 3.1. N´umero de cuadrados latinos reducidos de orden n................ 28 4.1. Resultados del algoritmo exacto del Axial 3AP . . . . . . . . . . . . . . . . . . . 41 4.2. Resultados del algoritmo exacto del Planar 3AP . . . . . . . . . . . . . . . . . . 41 4.3. Resultados del algoritmo de subgradiente del Axial 3AP . . . . . . . . . . . . . 44 4.4. Resultados del algoritmo de subgradiente del Axial 3AP . . . . . . . . . . . . . 45 4.5. Resultados del algoritmo de subgradiente del Axial 3AP . . . . . . . . . . . . . 45 4.6. Resultados del algoritmo de subgradiente del Axial 3AP sin dividir la relajaci´on ensubproblemas ................................... 46 4.7. Resultados del algoritmo de subgradiente del Axial 3AP al resolver la relajaci´on convariablescontinuas................................ 47 53 54 Bibliograf´ıa [1] S.E. Bammel and J. Rothstein. The number of 9 x 9 latin squares. Discr. Math, 1995. [2] Mokhtar S. Bazaraa and Hanif D. Sherali. On the choice of step size in subgradient optimization. European Journal of Operational Research, 7(4):380–388, 1981. [3] R. E. Burkard and R. Rudolf. Computational investigations on 3-dimensional axial assignment problems. JORBEL - Belgian Journal of Operations Research, Statistics, and Computer Science, 32(1, 2):85?98, Jan. 2020. [4] R.E. Burkard, M. Dell’Amico, and S. Martello. Assignment Problems. SIAM e-books. Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), 2009. [5] S Daskalaki, T Birbas, and E Housos. An integer programming formulation for a case study in university timetabling. European Journal of Operational Research, 153(1):117–135, 2004. Timetabling and Rostering. [6] S.A. Dichkovskaya and M.K Kravtsov. Investigation of polynomial algorithms for solving the three-index planar assignment problem. 2006. [7] D. Magos. Tabu search for the planar three-index assignment problem. J. Glob.Optim., 1996. [8] D. Magos and P. Miliotis. An algorithm for the planar three-indexassignment problem. European Journal of Operational Research, 1994. [9] Brendan D. McKay, Alison Meynert, and Wendy Myrvold. Small latin squares, quasigroups, and loops. Journal of Combinatorial Designs, 2007. [10] George L. Nemhauser and Laurence A. Wolsey. Integer and combinatorial optimization. Wiley- Interscience, New York, NY, USA, 1988. [11] S. Nickel, J. Puerto, and A. Rodriguez-Chia. Location problems with multiple criteria. In Gilbert Laporte, Stefan Nickel, and Francisco Saldanha da Gama, editors, Location Science, chapter 9, pages 205–247. Springer International Publishing Switzerland, 2015. [12] Francisco A. Ortega, Miguel A. Pozo, and Justo Puerto. Modelling and planning public cultural schedules for efficient use of resources. Computers and Operations Research, 58:9–23, 2015. [13] Dilson Lucas Pereira, Julio Cesar Alves, and Mayron Cesar de Oliveira Moreira. A multiperiod workforce scheduling and routing problem with dependent tasks. Computers and Operations Research, 118:104930, 2020. [14] Miguel A. Pozo, Justo Puerto, and Antonio M. Rodr´ıguez-Ch´ıa. The ordered median tree of hubs location problem. TOP, 29:78–105, 2021. [15] Robert J. Vanderbei. Linear ProgrammingFoundations and Extensions. 2014. [16] M. Vlach. Branch-and-bound method for the three-index assignment problem. Ekonomicko- Matematicky Obzor, 1967. 55 [17] L.A. Wolsey. Integer Programming. Wiley Series in Discrete Mathematics and Optimization. Wiley, 1998. 56