Implementaciones paralelas para un problema de control de inventarios de productos perecederos

Implementaciones paralelas para un problema de control de inventarios de productos perecederos Alejandro G. Alcoba, Eligius M.T. Hendrix, Inmaculada Garc´ıa1, Gloria Ortega 2 Resumen—En este trabajo se analizan y eval´uan dos implementaciones de un algoritmo de optimizaci´on para un problema de control de inventarios de productos perecederos. Las implementaciones se han llevado a cabo utilizando una arquitectura heterog´enea donde cada nodo est´a compuesto por varios multicores y varias GPUs. Las versiones paralelas que se han desarrollado son: (1) una versi´on MPI-PTHREADS en la que se extrae el paralelismo tanto a nivel de proceso MPI como a nivel de hilo y (2) una versi´on multiGPU en la que se obtiene el paralelismo a nivel de proceso MPI y a nivel de cores de GPU. Este algoritmo puede ser descompuesto f´acilmente en un conjunto de tareas que no presentan ninguna dependencia entre s´ı. Sin embargo, la carga computacional asociada a cada una de las tareas es diferente y el problema del reparto de las tareas entre los elementos de proceso se puede modelar como un problema de Bin Packing. Ello implica que la selecci´on del conjunto de tareas asociadas a cada una de las unidades de computaci´on requiere del dise˜no de heur´ısticas que sean capaces de balancear la carga eficientemente y de forma est´atica. En este trabajo hemos analizado y evaluado varias heur´ısticas. Finalmente, la mejor heur´ıstica ha sido la utilizada en la implementaci´on paralela del algoritmo de control de inventarios que ha sido evaluado en la versi´on MPI-PTHREADS y en la versi´on multiGPU. Para la implementaci´on MPI-PTHREADS los resultados obtenidos muestran una buena escalabilidad mientras que las versi´on MultiGPU para el ejemplo que se ha evaluado deja de ser eficiente cuando se usan mas de 2 GPUs. Palabras clave—Multihilo, Multi-GPU, Bin Packing, Monte-Carlo, inventarios, productos perecederos. I. Introducci´ on EL objetivo de este trabajo consiste en determinar hasta que punto el uso de una arquitectura heterog´enea (multicore-multiGPU) facilita la resoluci´on de un problema de optimizaci´on del control de inventarios de productos perecederos. Pretendemos aprovechar la capacidad computacional de estas arquitecturas para obtener soluciones mas exactas, para ejemplos mas pesados desde el punto de vista de la computaci´on, manteniendo tiempos de respuesta aceptables. El problema del control de inventarios queda definido a lo largo de una serie finita de Tperiodos de tiempo en los que se ha de satisfacer la demanda (estoc´astica) de un determinado producto perecedero que desde que se produce tiene una vida ´util de Jperiodos. En el modelado de este problema se supone que la distribuci´on se realiza siguiendo la pol´ıtica de distribuci´on FIFO, entregando 1Computer Architecture, Universidad de M´alaga. Campus de Excelencia Internacional Andaluc´ıa Tech, e-mail: {agutierreza,eligius,igarciaf}@uma.es 2Informatics, Univ. of Almer´ıa, Agrifood Campus of Int. Excell., ceiA3, e-mail: [email protected] el producto demandado con mayor antig¨uedad. Se supone, adem´as, que la demanda que no se satisfaga en un periodo queda perdida, no pudi´endose acumular al periodo siguiente. La soluci´on a este problema consiste en encontrar que cantidades de pedido a lo largo de todos los periodos resulta ´optima, en el sentido de minimizar el coste asociado a la producci´on, distribuci´on, almacenamiento y desecho de los productos que sobrepasen su vida ´util. Actualmente, las arquitecturas de computaci´on de altas prestaciones m´as extendidas son las plataformas heterog´eneas basadas en sistemas de memoria distribuida, donde cada nodo tiene una arquitectura multicore que podr´ıa albergar un n´umero distinto de cores [1]. Por lo tanto, las implementaciones paralelas tienen que ser adaptadas para poder ser ejecutadas en dichas arquitecturas heterog´eneas. En este contexto, es necesario tener un conocimiento detallado tanto del algoritmo a paralelizar como de los recursos computacionales que se van a utilizar para la implementaci´on [2]. Adem´as, a estas arquitecturas se les pueden incorporar aceleradores, como son FPGAS, GPUs, coprocesadores Intel Xeon Phi, etc. En concreto, en el problema del control de inventario para productos perecederos se ha optado por la combinaci´on de cl´usteres de Multi-GPUs. De este modo, el uso de plataformas masivamente paralelas (GPUs) permite la aceleraci´on de las tareas computacionalmente m´as costosas, porque estas unidades tienen mucha potencia de c´alculo para los esquemas de computaci´on vectorial. De forma adicional, el uso de plataformas de memoria distribuida permite obtener unos resultados m´as precisos debido a que el uso de computaci´on paralela permite incrementar el n´umero de simulaciones realizadas para resolver un caso particular sin que el tiempo de ejecuci´on se incremente. El modelo de computaci´on paralela asociado a este problema se puede describir en t´erminos de un conjunto de conjunto de tareas que no presentan dependencias entre s´ı. Sin embargo, la carga computacional de cada una de estas tareas es variable y por lo tanto pueden aparecer problemas de desbalanceo de la carga si se hace un reparto de la carga a ciegas. Este problema de asignaci´on de tareas a elementos de procesamiento se conoce en la literatura de complejidad como problema de Bin packing [3]. Dado que es un problema NP-Completo, se han desarrollado varias heur´ısticas que permiten tener una soluci´on en un tiempo razonable. El resto del trabajo se organiza de este modo. La Secci´on II estudia y analiza el modelo propuesto para el control de inventario de productos perecederos. La Secci´on III describe el algoritmo secuencial implementado. En la Secci´on IV se explican los detalles de las versiones paralelas implementadas. Adem´as, se describe el problema de Bin packing que hay que resolver inicialmente para repartir la carga de trabajo entre las diferentes plataformas. La Secci´on V ofrece algunos resultados experimentales computacionales al evaluar el modelo utilizando un cl´uster de Multi-GPUs. Finalmente, la Secci´on VI expone las conclusiones y las principales l´ıneas de actuaci´on futuras. II. Descripci´ on del modelo La base de las implementaciones que se presentan en este trabajo es un algoritmo desarrollado en Matlab para resolver un problema MINLP (Mixed Integer NonLinear Programming). Se trata de planificar, a lo largo de un n´umero finito de periodos T, las cantidades que se deben proveer de cierto producto perecedero para satisfacer la demanda bajo una restricci´on que establece un nivel de servicio βque necesariamente se debe satisfacer. En concreto, esta restricci´on establece que para cada periodo (siempre hablando en t´erminos de esperanza matem´atica) a lo sumo una fracci´on βde la demanda no pueda ser satisfecha y sea perdida por falta de stock, ya que se supone que esta no puede ser servida en un periodo posterior. Esta condici´on es equivalente a que al menos una fracci´on (1 −β) de la demanda sea cubierta en cada periodo. La duraci´on de cada item producido desde que est´a disponible para el consumidor hasta que ha de ser retirado es de J < T periodos. Adem´as, se supone que los productos se distribuyen siguiendo la regla FIFO: los productos son expedidos comenzando por los m´as antiguos. El problema de optimizaci´on que se plantea es el de encontrar la cantidad de producto perecedero que hay que producir en cada periodo de forma que se satisfgan todas las restricciones del problema y que adem´as se minimice una funci´on coste. A continuaci´on se detallan las principales variables del modelo: Indices t´ındice del periodo, t= 1, . . . , T, siendo Tel n´umero total de periodos j´ındice de edad, j= 1, . . . , J, siendo Jla vida ´util de cada unidad Data dtDemanda en cada periodo con distribuci´on normal dada por su media µt>0 y varianza (cv ×µt)2dado por un coeficiente de variaci´on cv, id´entico en cada periodo. kCoste por periodo en el que se decide realizar un pedido, k > 0 cCoste unitario de producto, c > 0 hCoste por almacenamiento, h > 0 wCoste unitario de desecho, puede ser negativo con la condici´on, w > −c βNivel de servicio, 0 < β < 1 Variables Qt≥0 Cantidad de producto producido y disponible en el periodo t. Denotamos por Qal vector completo (Q1, . . . , QT) Yt∈ {0,1}Indica si se produce un pedido en el periodo t. Es 1 si y solo si Qt>0. Denotamos por Yal vector completo (Y1, . . . , YT) XtVentas perdidas en el periodo t Ijt Inventario de edad jal final del periodo t, considerando un periodo inicial fijo, Ij0= 0, Ijt ≥0 para j= 1, . . . , J. Adem´as, se usar´a la notaci´on (·)+=max(·,0). La funci´on coste que se pretende minimizar depende del vector Q= (Q1, . . . , QT) y se puede definir como: f(Q) = T X t=1  C(Qt) + E h J−1 X j=1 Ijt +wIJt  , (1) siendo C(x) = k+cx, if x > 0,and C(0) = 0.(2) El nivel de inventario para cada periodo t= 1, . . . , T y cada edad jsiguiendo la regla FIFO puede calcularse como sigue: Ijt =          Qt−(dt−PJ−1 j=1 Ij,t−1)++ j= 1, (IJ−1,t−1−dt)+j=J, Ij−1,t−1−(dt−PJ−1 i=jIi,t−1)++ otro j (3) Por otra parte, la restricci´on del nivel de servicio puede expresarse como: E(Xt)≤(1 −β)µt, t = 1, . . . , T. (4) Para controlar el cumplimiento de esta restricci´on es necesario calcular las ventas perdidas que se producen en cada periodo t, lo cual viene dado por: Xt= dt− J−1 X j=1 Ij,t−1−Qt  + .(5) El valor esperado de las ventas perdidas es una funci´on conocida como loss-function que en general no admite una expresi´on en t´erminos elementales. Algunas aproximaciones factibles pueden verse en [4], [5], [6], [7]. Para nuestro modelo hemos decidido utilizar la simulaci´on Monte-Carlo para obtener una estimaci´on de la loss-function. Con las condiciones impuestas, el problema de encontrar las cantidades de producto perecedero que se deben producir en cada periodo y que minimizan la funci´on coste f(Q) dada en (1), y con ello la pol´ıtica Y∈ {0,1}Tde periodos de pedido ´optima, es un problema MINLP (Mixed Integer NonLinear Programming). Su resoluci´on, mediante t´ecnicas de programaci´on din´amica, puede consultarse en [8]. Como veremos, la t´ecnica usada presenta caracter´ısticas adecuadas para su implementaci´on en computadores de alto rendimiento. III. Algoritmo secuencial En esta secci´on se detalla el algoritmo que resuelve el problema MINLP mediante programaci´on din´amica, tal como se ha planteado en la secci´on anterior. El Algoritmo 1 presenta el nivel m´as general de la resoluci´on del problema. En primer lugar, se generan en forma de vectores todas las pol´ıticas de pedido v´alidas {0,1}Tpara unos valores TyJdados. Una pol´ıtica de pedido Y∈ {0,1}Tes considerada no v´alida si contiene Jo m´as ceros consecutivos, ya que se considera que en ning´un periodo la demanda puede ser determinista e igual a 0 y, por lo tanto, no se cumplir´ıan los requerimientos de nivel de servicio. Para una determinada pol´ıtica de pedidos Y, el Algoritmo 2 calcula, para cada periodo, cuales son las cantidades de producto a producir (Q(Y)) que minimizan la funci´on coste (1). El Algoritmo 2 solo se ejecuta para los casos en los que la funci´on LB(Y)<mincost, siendo LB(·) un limite inferior de la estimaci´on de la funci´on coste (1) (en [9] puede encontrarse informaci´on sobre la estimaci´on de LB(.)). Por lo tanto, el Algoritmo 1 realiza una evaluaci´on exhaustiva (salvando una cantidad relativa de casos con la funci´on LB), para encontrar tanto la pol´ıtica ´optima Y∗como el vector de producci´on ´optimo Q∗ y su coste asociado f(Q∗). Algorithm 1 AllY (): C´alculo de la pol´ıtica de pedidos ´optima (Y∗) a partir de la evaluaci´on de todas las politicas factibles Y 1: Generate all feasible Y 2: mincost=∞ 3: for all Ydo 4: if LB(Y) <mincost then 5: QY=MINQ(Y); # Algorithm 2 6: Determine f(QY) 7: if f(QY)<mincost then 8: mincost=f(QY) 9: Y∗=Y;Q∗=QY 10: end if 11: end if 12: end for 13: return Y∗;Q∗;f(Q∗) = mincost Algorithm 2 MinQ(Y): Optimizaci´on de Q(Y) 1: Generar para Ylos vectores AyR 2: for i= 1 to lenght(A)do 3: if Ai= 1 o Ai−Ai−1=Jthen # No Invent. 4: Qi=QRi,Ai; 5: floss(Qi, Ai, Ai+Ri−1) #Alg 4 6: else 7: Qi=OrdV al(Ai, Ai+Ri) # Alg 3 8: end if 9: end for 10: return Q; Para cada periodo ten el que Yt= 1, el Algoritmo 2 calcula la cantidad de producto producido que permita cubrir la demanda de un ciclo, entendiendo por ciclo los periodos comprendidos entre dos elementos consecutivos del vector Yque valen 1. El vector A almacena el ´ındice del inicio de cada ciclo, mientras que Rguarda su duraci´on en periodos. Para cada ciclo, la cantidad de producto necesaria es aquella cantidad que ajusta el valor de las ventas perdidas al l´ımite inferior permitido por (4). La forma de calcular Qi(Y) para los ciclos en los que el inventario es nulo es diferente de la usada para los ciclos en los que existe alg´un inventario. En caso de no existir inventario anterior el valor esperado de las ventas perdidas para este problema concreto puede ser calculado resolviendo una sencilla ecuaci´on [8]. Para los ciclos en los que existe alg´un inventario es necesario recurrir a la simulaci´on reiterada del problema, haciendo uso de (3) y (5) para obtener aproximaciones del valor de ventas perdidas que permitan encontrar el valor ´optimo de Q(Y). De eso se encarga la funci´on OrdV al(t1, t2), descrita en el Algoritmo 3, que, por el m´etodo de la secante, calcula la cantidad ´optima para el ciclo. La funci´on floss(q, t1, t2) (Algoritmo 4 realiza la simulaci´on y devuelve la cantidad de ventas perdidas en el ´ultimo periodo del ciclo, t2, suponiendo una llegada de qunidades al comienzo del ciclo (periodo t1). Algorithm 3 OrdV al(t1, t2): Calculo de la cantidad de pedido m´ınima que satisface (4) para un ciclo entre t1hasta t2 1: K= (1 −β)∗µt2; 2: q1= t2 P i=t1 µi; 3: q2=q1(1 + cvˆx); 4: f1←floss(q1, t1, t2); #Alg 4 5: f2←floss(q2, t1, t2); #Alg 4 6: repeat 7: q←q1+(K−f1)(q2−q1) f2−f1; #secant method 8: f1←f2; 9: f2←floss(q, t1, t2); #Alg 4 10: q1←q2; 11: q2←q; 12: until |f2 µt2 | − (1 −β)<  13: return q Algorithm 4 floss(q, t1, t2): M´etodo Monte-Carlo para aproximar E(X) 1: for n= 1 to Ndo 2: for t=t1to t2do 3: for j= 1 to j=Jdo 4: Update Ij,t,n using (3) 5: end for 6: end for 7: end for 8: X=1 N N P n=1 dt2,n −PJ−1 j=1 Ij,t2−1,n −q+; 9: return X; IV. Implementaciones paralelas El orden de complejidad del problema, partiendo del Algoritmo 1, est´a relacionado con el n´umero de pol´ıticas factibles Y, lo cual depende de los valores de JyT. Independientemente del valor de J, el n´umero de casos posibles a tratar aumenta de forma exponencial con el valor de T, es decir, O(eT). M´as a´un, para cada Ytratado por el Algoritmo 2, su complejidad depende del n´umero de veces que es necesario utilizar el Algoritmo 3, limitada a Ten cada caso. Por ´ultimo, cada ejecuci´on del Algoritmo 3 requiere el c´omputo de la funci´on floss que a su vez ejecuta N simulaciones de Monte-Carlo en las que se tiene que calcular el inventario y las ventas perdidas mediante (3) y (5). Por tanto, el orden de complejidad para el m´etodo completo, es decir, hallar el vector de pedidos Y´optimo y las cantidades ´optimas de pedido, es aproximadamente del orden de O(N·T·eT). En la secci´on anterior se ha puesto de manifiesto la necesidad de realizar simulaciones para obtener aproximaciones de la funci´on floss. Esta funci´on es la que consume la mayor parte del tiempo computacional de la ejecuci´on del problema de inventarios. Es importante destacar la necesidad de realizar un elevado n´umero de simulaciones para que las aproximaciones que lleva a cabo la funci´on floss sean suficientemente exactas. Como ejemplo, en la Figura 1 se muestra el histograma de 20000 ejecuciones de la funci´on floss, realizando en cada una de ellas N= 30000 simulaciones de Monte-Carlo, para un caso en el que anal´ıticamente, su valor exacto es 0.05. Las aproximaciones de la funci´on floss siguen una distribuci´on normal (pasa un test chi-cuadrado (χ2) al 5% de nivel de significaci´on) y, aunque la media se ajusta a la realidad, 0.05, la dispersi´on es alta. Por tanto, para realizar aproximaciones relativamente precisas del valor de las ventas perdidas (X), es necesario realizar un n´umero Nde simulaciones del problema suficientemente alto, lo que constituye la verdadera carga computacional del problema. Como caso de estudio, se ha considerado un ejemplo del problema de control de inventarios en el que T= 12 y J= 3, que es bastante realista. Consideraremos siempre N= 50000 simulaciones de Monte-Carlo. Revisando el Algoritmo 2 se observa que los ciclos en los que existe inventario previo son aquellos en los que hay que proceder a la simulaci´on de Monte-Carlo de forma reiterada mediante el Algoritmo 3 (funci´on OrdV al), que es lo que representa la mayor carga computacional de cada optimizaci´on de un vector Y. Podemos hacer una estimaci´on de la carga computacional asociada con cada vector de pedidos (Y) (Algoritmo 2) como la suma de los valores de los elementos del vector Rpara los casos en los que se ejecuta el Algoritmo 3 (funci´on OrdV al). Para el caso que vamos a analizar, con T= 12 y J= 3, el n´umero total de vectores Y(diferentes politicas de pedidos) que se generan es 927, y la carga computacional estimada asociada a cada vector Yvalorada entre 0 y 11 presenta una distribuci´on como la mostrada en Fig. 1 Histograma de 20000 aproximaciones de floss mediante el m´ etodo Monte-Carlo con N= 30000 simulaciones en cada caso. El valor exacto, calculado anal´ ıticamente para ese caso, es 0.05 la Tabla I, donde, por ejemplo, existen 24 configuraciones diferentes del vector Yque tendr´ıan una carga estimada de 5 unidades. Experimentalmente, se ha podido comprobar que esta estimaci´on de la carga para cada vector Yest´a directamente relacionada con el tiempo de ejecuci´on del Algoritmo 2. Una vez analizada la estructura algor´ıtmica del problema, pasamos a describir los detalles de las implementaciones que hemos llevado a cabo sobre una arquitectura heterog´enea formada por un cl´uster de Multi-GPUs (multicores y dispositivos GPUs). El hecho de explotar una plataforma heterog´enea de un cl´uster tiene dos ventajas fundamentales: poder abordar la resoluci´on de problemas de mayor tama˜no y reducir el tiempo de ejecuci´on de un caso concreto. Las implementaciones consideradas en este trabajo han sido: •MPI-PTHREADS: Esta implementaci´on obtiene el paralelismo de los procesadores multicore y de los nodos disponibles en el cl´uster. Para ello, se utiliza programaci´on basada en hebras [10] y MPI [11]. •Multi-GPU: Esta implementaci´on est´a basada en el uso de GPUs para realizar las simulaciones de Monte-Carlo, las cuales son la parte computacionalmente m´as costosa del problema a resolver. Para ello, la interfaz de programaci´on que se utiliza es CUDA [12]. En las siguientes subsecciones se describen ambas implementaciones, as´ı como el reparto inicial de la carga entre las diferentes unidades de proceso, en base a la carga estimada asociada a cada vector Y tal y como se ha descrito en la Tabla I, y que se puede modelar como un problema de Bin packing [3]. A. Implementaci´on MPI-PTHREADS Centrando nuestra atenci´on en la implementaci´on MPI-PTHREADS, se ha explotado el paralelismo en TABLA I N´ umero de vectores Y(Casos posibles) en funci´ on de la carga computacional estimada, para T= 12 yJ= 3. Carga estimada 0 1 2 3 4 5 6 7 8 9 10 11 Casos posibles 1 1 2 6 14 24 50 86 120 185 260 178 dos niveles: a nivel de nodo (memoria distribuida) y a nivel de multicore (memoria compartida). Por un lado, existen m´ultiples formas de paralelizar rutinas en modelos de memoria compartida, aunque la librer´ıa est´andar es Pthreads (POSIX threads). Pthreads provee un conjunto unificado de rutinas en una librer´ıa de C cuyo principal objetivo es facilitar la implementaci´on de threads o hilos en el programa. Por otro lado, debido a su portabilidad, MPI ha sido el interfaz considerado para explotar el paralelismo a nivel de nodo. Partiendo del algoritmo de optimizaci´on del problema de inventarios, se ha realizado una paralelizaci´on h´ıbrida (MPI y Pthreads), en la cual el conjunto de vectores Yque se van a evaluar en el problema de optimizaci´on son repartidos entre los procesadores de acuerdo a las heur´ısticas de balanceo de la carga que se describen en la Secci´on IV-B. La evaluaci´on de esta implementaci´on se ha realizado en un cl´uster Bullx y los resultados se describen en la Secci´on V. B. Heur´ısticas para el problema de Bin packing El reparto inicial de la carga de trabajo entre los procesadores disponibles puede considerarse como un problema de Bin packing con algunas restricciones. El problema de Bin packing se enmarca dentro de la optimizaci´on combinatoria (NP-completo), y en nuestro caso se puede modelar de la siguiente forma: Dado un conjunto de Eejecuciones independientes del Algoritmo 2 (items), cada una de ellas con una carga computacional 0 < wi< B y dado un conjunto de Pprocesadores (Bins), repartir las ejecuciones del algoritmo entre los procesadores de forma que la carga computacional m´axima asignada a un procesador sea m´ınima (ver [3] para una formulaci´on general del problema de Bin packing). Debido a la dificultad de encontrar soluciones ´optimas para este tipo de problemas, habitualmente se utilizan t´ecnicas heur´ısticas y metaheur´ısticas, que son capaces de encontrar una soluci´on aceptable en un tiempo razonable. Algunas de estas heur´ısticas est´an inspiradas en computaci´on evolutiva [13]. Para resolver el problema de balanceo de la carga que se ha modelado como un problema de tipo Bin packing, proponemos tres algoritmos heur´ısticos (H1, H2 y H3) para repartir la carga de trabajo (en nuestro caso, los posibles vectores Y) entre todos los elementos de procesamiento disponibles de forma que se minimice el tiempo de ejecuci´on del problema de optimizaci´on. 1. (H1): Heur´ıstica basada en Round Robin: Ordenando previamente, de mayor a menor, el peso de las tareas a asignar, estas se reparten entre los Pprocesadores siguiendo el patr´on (1, . . . , P, P, P −1,...,1,1, . . .) 2. (H2): Heur´ıstica basada en asignar sucesivamente los items wial procesador que menos carga de trabajo haya acumulado. 3. (H3): Similar a la heur´ıstica H2, pero previamente ordenando los items de mayor a menor carga. Para valorar las heur´ısticas se han utilizado tres instancias del problema denominadas γ1,γ2y U, en las que la carga computacional estimada que se asocia a cada vector Yes diferente. (γ1)y(γ2) est´an basadas en distribuciones gamma con par´ametros de forma y escala (10,4) y (1,25), respectivamente, y (U) sigue una distribuci´on uniforme con valores entre 0 y 100. De cada una se han tomado 4000 muestras distintas, entendiendo que cada una de ellas es una tarea cuya carga computacional est´a asociada a su valor. En la Figura 2 aparecen representadas las muestras de las distribuciones gamma(10,4), gamma(1,25) y la uniforme. Para medir el grado de balanceo de la carga asignada a cada procesador, se ha utilizado el coeficiente de Gini (G), ampliamente utilizado en el campo de la econom´ıa para medir el grado de desigualdad de la distribuci´on de la riqueza en poblaciones [14], [15], [16]. Este ´ındice var´ıa entre 0 (equidad absoluta) y 1 (un solo individuo (procesador) posea toda la riqueza de la poblaci´on (carga computacional)). G se define como la media de la diferencia entre cada posible par de procesadores, divididos por su carga media. Para un n´umero de ejecuciones Easignadas a Pelementos de proceso, siendo wila carga computacional asignada al procesador i, ordenadas de forma ascendente, Gse calcula como sigue: G= 2 P P i=1 i·wi P P P i=1 wi −P+ 1 P(6) Gr´aficamente, Grepresenta el ratio entre la diferencia del ´area rodeada por la l´ınea de uniformidad y la curva de Lorenz de la distribuci´on, y el ´area triangular que hay debajo de la l´ınea de uniformidad. G toma valores entre un m´ınimo de 0, cuando todos los procesadores tienen la misma carga, a un m´aximo de 1, cuando todos los procesadores (excepto uno) tienen una carga de cero. Por lo tanto, cuando Gse acerca a 0 la carga est´a bien balanceada, y cuando se acerca a 1 est´a desbalanceada. En la Tabla II se resume el comportamiento de las diferentes heur´ısticas a trav´es del valor del coeficiente de Gini (G) para los ejemplos planteados (compar´andolos con un reparto a ciegas (columna Fig. 2 Distribuciones γ1: Gamma(10,4) (arriba), γ2: Gamma(1,25) (centro) y U Uniforme[0,100] (abajo) de una muestra de 4000 casos. (HR)). Claramente se demuestra en esta tabla que la heur´ıstica H3 es, al menos, un orden de magnitud mejor que las heur´ısticas H1 y H2 y que el reparto aleatorio de las tareas entre los procesadores (HR) es al menos dos ordenes de magnitud peor que H3 y un orden de magnitud peor que H1 y H2. De los datos de la Tabla II, se concluye que la heur´ıstica H3 es la que presenta mejores resultados, consiguiendo balancear la carga de forma casi exacta, por lo tanto esta es la heur´ıstica que produce mejores tiempos de ejecuci´on en la evaluaci´on de las implementaciones paralelas del problema del control de inventarios que se muestra en la Secci´on V. TABLA II Comportamiento de las heur´ ısticas H1, H2 y H3 y su comparativa con un reparto aleatorio de tareas HR, para P= 8,16,32,64 y las distribuciones γ1,γ2y U. El comportamiento se mide por el coeficiente de Gini G. γ1 PH1 H2 H3 HR 8 1,3 10−44,8 10−41,2 10−53,7 10−3 16 2,5 10−41,2 10−35,2 10−57,6 10−3 32 5,8 10−42,0 10−31,4 10−41,3 10−2 64 3,2 10−33,9 10−31,5 10−32,1 10−2 γ2 PH1 H2 H3 HR 8 6,9 10−41,0 10−31,9 10−52,3 10−2 16 1,2 10−32,1 10−31,9 10−53,9 10−2 32 3,3 10−34,1 10−37,8 10−56,1 10−2 64 8,0 10−37,6 10−31,0 10−47,9 10−2 U PH1 H2 H3 HR 8 1,3 10−47,5 10−47,4 10−61,2 10−2 16 1,8 10−41,4 10−38,6 10−61,8 10−2 32 2,2 10−42,4 10−33,9 10−52,6 10−2 64 3,7 10−44,1 10−35,4 10−54,0 10−2 C. Implementaci´on Multi-GPU La versi´on Multi-GPU se ha basado en la explotaci´on de diversas GPUs para la paralelizaci´on de las simulaciones del m´etodo de Monte-Carlo, realizadas por la funci´on flossGPU (ver Algoritmo 5). Por una parte, cada una de las Nsimulaciones son independientes entre s´ı. Al mismo tiempo, para el c´alculo de todo el inventario de cada una de las posibles edades (3), este solo depende del inventario del periodo anterior. Por tanto, separando por periodos, un kernel de CUDA puede realizar en paralelo el c´alculo de las Nsimulaciones y, al mismo tiempo, la actualizaci´on del inventario de Jedades diferentes. Por tanto, la computaci´on que se realiza con la GPU es la simulaci´on de Monte-Carlo. Al mismo tiempo, el modelo se ha implementado de forma que cada proceso MPI abre uno o dos hilos, que a su vez abren una o dos GPUs del nodo en el que se encuentra. Las tareas se reparten usando la heur´ıstica (H3) entre las GPUs que se vayan a utilizar en cada caso. El Algoritmo 5 resume el cambio realizado en el Algoritmo 4 para adaptarlo a su ejecuci´on en una o varias GPU. Tal y como se ha mencionado, el bucle que recorre los periodos se sit´ua en el primer nivel. En la siguiente secci´on se muestran los resultados experimentales de la ejecuci´on de las implementaciones MPI-PTHREADS y Multi-GPU. V. Resultados Para la evaluaci´on de las implementaciones paralelas hemos utilizado un cl´uster compuesto de ocho nodos Bullx R424-E3 Intel Xeon E5 2650 (cada uno con 16 cores), interconectados por un puerto Infini- Fig. 3 (a) Cuatro nodos del cluster Multi-GPU utilizado para la evaluaci´ on del problema de control de inventarios y (b) caracter´ ısticas de las GPUs. Algorithm 5 flossGPU(q, t1, t2): M´etodo MonteCarlo para aproximar E(X) 1: for t=t1to t2do 2: for n= 1 to Ndo #GPU (l´ıneas 2-6) 3: for j= 1 to j=Jdo 4: Update Ij,t,n using (3) 5: end for 6: end for 7: end for 8: X=1 N N P n=1 dt2,n −PJ−1 j=1 Ij,t2−1,n −q+; 9: return X; Band QDR/FDR embebido en la placa madre, 8-GB RAM y 16-GB SSD) con ocho GPUs TeslaM2075 (de los ocho nodos, cuatro de ellos tienen dos GPUs por nodo). El driver de CUDA que se ha utilizado es CUDA 6.5 [17]. La arquitectura Multi-GPU y las caracter´ısticas de las GPUs se muestran en la Figura 3. Las Tablas III, IV y V muestran los tiempos de ejecuci´on, en segundos, del ejemplo del problema de control de inventarios de productos perecederos descrito anteriormente, en el que T= 12 y J= 3. Adicionalmente, se ha medido el tiempo tomado por cada heur´ıstica para asignar las tareas a los procesadores, siendo en el peor de los casos de aproximadamente mil microsegundos, por lo que todos ellos cumplen su funci´on de balancear la carga en un tiempo que resulte insignificante para el problema. Cada proceso MPI se ejecuta en un nodo diferente del cl´uster (hemos evaluado con hasta ocho nodos), mientras que en cada uno de ellos se ha probado con hasta diecis´eis hilos. La multiplicaci´on del n´umero de procesos MPI por el n´umero de hilos, da como resultado el n´umero de cores entre los que las heur´ısticas distribuiran la carga computacional asociada a los vectores Ydel problema. En la tabla VI se recogen los tiempos cuando se aplica un reparto aleatorio de la carga entre los cores (HR). TABLA III Tiempos de ejecuci´ on del Algoritmo 1 usando la heur´ ıstica H1 (en segundos). MPI\threads 1 2 4 8 16 1 63.10 37.42 18.79 8.07 4.96 2 31.75 18.82 9.56 4.74 2.61 4 16.00 9.60 4.96 2.56 1.42 8 8.17 4.98 2.63 1.43 0.88 TABLA IV Tiempos de ejecuci´ on del Algoritmo 1 usando la heur´ ıstica H2 (en segundos). MPI\threads 1 2 4 8 16 1 63.18 37.35 18.38 8.04 4.94 2 31.76 18.77 9.51 4.83 2.62 4 16.04 9.53 4.86 2.55 1.48 8 8.14 4.89 2.57 1.40 0.92 TABLA V Tiempos de ejecuci´ on del Algoritmo 1 usando la heur´ ıstica H3 (en segundos). MPI\threads 1 2 4 8 16 1 63.14 37.22 18.57 8.04 4.80 2 31.72 18.76 9.51 4.77 2.59 4 16.01 9.49 4.87 2.21 1.42 8 8.14 4.89 2.54 1.40 0.87 TABLA VI Tiempos de ejecuci´ on del Algoritmo 1 usando la heur´ ıstica HR (en segundos). MPI\threads 1 2 4 8 16 1 63.08 38.03 19.49 8.34 4.98 2 32.19 19.51 9.94 4.94 2.83 4 16.49 9.95 5.05 2.60 1.36 8 8.43 5.06 2.64 1.43 0.96 En cuanto a la paralelizaci´on MPI-PTHREADS, se puede apreciar en las tablas un buen nivel de speedup para todas las heur´ısticas consideradas. Aunque en diversas pruebas, simulando el comportamiento de las heur´ısticas con problemas con una gran disparidad de carga computacional de las tareas (distribuciones γ1,γ2y uniforme) se pueden apreciar diferencias entre ellas, para el problema de inventarios los tiempos de ejecuci´on son similares, ligeramente inferiores para la heur´ıstica H3. En t´erminos relativos, tambi´en puede apreciarse una diferencia en tiempos de ejecuci´on con respecto al reparto a ciegas (HR). En la Tabla VII se recogen los tiempos de ejecuci´on del problema completo, tomando hasta las 8 GPUs existentes en el cl´uster para la paralelizaci´on del m´etodo Monte-Carlo en CUDA. En dicha tabla se aprecia como se pueden conseguir buenos resultados paralelizando una escasa porci´on del c´odigo. En cambio, al usar m´ultiples GPUs, la escalabilidad est´a penalizada por el tiempo de inicializaci´on de las GPUs (aproximadamente 5 segundos), lo cual hace que no se obtenga un buen rendimiento con m´as de 2 GPUs para este ejemplo concreto. Para el ejemplo de inventarios considerado, el uso de los diecis´eis cores de un solo nodo resulta m´as beneficioso que el uso de las dos GPUs disponibles. TABLA VII Tiempo en segundos para la paralelizaci´ on del m´ etodo Monte-Carlo para GPUs. Number of GPUs Tiempo 1 13.54 2 9.53 4 7.44 6 6.80 8 6.39 VI. Conclusiones En este trabajo se ha analizado la implementaci´on paralela del problema basado en el control de inventarios de productos perecederos. Dicha implementaci´on se basa en tres ideas principales: 1. La explotaci´on de plataformas multicore, bas´andonos en el uso de Pthreads y MPI. 2. El uso de la computaci´on GPU para acelerar la parte computacionalmente m´as costosa. 3. El uso de computaci´on MPI para distribuir la carga entre los distintos procesadores y poder utilizar GPUs ubicadas en diferentes nodos. La paralelizaci´on MPI-PTHREADS con un balanceo est´atico de la carga basado en heur´ısticas presenta una buena escalabilidad. La rapidez con la que las heur´ısticas se ejecutan las hace apropiadas incluso para problemas en los que el desbalanceo no es muy acusado, mejorando un reparto aleatorio de las tareas. Para casos en los que el desbalanceo es mucho mayor (ejemplos con las distribuciones gamma y uniforme) el beneficio es mucho mayor y se observa que la heur´ıstica (H3) presenta mejores resultados. Con respecto a la paralelizaci´on del m´etodo Monte-Carlo en CUDA, la escalabilidad con el uso de m´ultiples GPUs se ve limitada por el tiempo de inicializaci´on requerido, aunque podr´ıa ser beneficiosa en comparaci´on con la paralelizaci´on MPI-PTHREADS si el problema a tratar requiere m´as precisi´on, siendo necesarias m´as simulaciones del m´etodo de MonteCarlo, de forma que el tiempo de inicializaci´on de las GPUs se hiciera comparativamente irrelevante. Agradecimientos Alejandro G. Alcoba es becario del programa FPI. Este trabajo ha sido financiado por Universidad de M´alaga, Campus de Excelencia Internacional Andaluc´ıa Techel, el Ministerio de Ciencia (TIN201237483) y la Junta de Andaluc´ıa (P11-TIC7176), parcialmente financiados por el Fondo Europeo de Desarrollo Regional (FEDER). Referencias [1] J. L. Hennessy and D. A. Patterson, Computer Architecture: A Quantitative Approach, Morgan Kaufmann, 2011. [2] A.L. Lastovetsky, “Special issue on heterogeneity in parallel and distributed computing,” Journal of Parallel and Distributed Computing., vol. 73, no. 12, 2013. [3] M. R. Garey and D. S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness, W. H. Freeman & Co., New York, NY, USA, 1979. [4] A.A. Kurawarwala and H. Matsuo, “Forecasting and inventory management of short life-cycle products,” Operations Research, vol. 44, pp. 131–150, 1996. [5] R. Rossi, S. Armagan Tarim, S. Prestwich, and B. Hnich, “Piecewise linear lower and upper bounds for the standard normal first order loss function,” Applied Mathematics and Computation, vol. 231, pp. 489–502, 2014. [6] S. K. De Schrijver, El-H. Aghezzaf, and H. Vanmaele, “Double precision rational approximation algorithm for the inverse standard normal first order loss function,” Applied Mathematics and Computation, vol. 219, no. 3, pp. 1375 – 1382, 2012. [7] G. R. Waissi and D. F. Rossin, “A sigmoid approximation of the standard normal integral,” Applied Mathematics and Computation, vol. 77, no. 1, pp. 91 – 95, 1996. [8] A.G. Alcoba, E.M.T. Hendrix, I. Garcia, K.G.J. PaulsWorm, and R. Haijema, “On computing order quantities for perishable inventory control with nonstationary demand,” in Proceedings of MAGO 2014, 2014, pp. 5–8. [9] A.G. Alcoba, E.M.T. Hendrix, I Garcia, G. Ortega, K. G. J. Pauls-Worm, and R. Haijema, “On computing order quantities for perishable inventory control with non-stationary demand,” in Proceedings of ICCSA 2015. 2015, vol. En prensa of Lecture Notes in Computer Science, Springer. [10] D.R. Butenhof, Programming with POSIX Threads, Professional Computing Series. Addison-Wesley, 1997. [11] M. Snir, S. Otto, S. Huss-Lederman, D. Walker, and J. Dongarra, MPI-The Complete Reference, Volume 1: The MPI Core, MIT Press, Cambridge, MA, USA, 1998. [12] “CUDA C Best Practices Guide,” http://docs.nvidia. com/cuda/cuda-c-best-practices-guide/, accessed 23 February 2015. [13] C. Blum and A. Roli, “Metaheuristics in combinatorial optimization: Overview and conceptual comparison,” ACM Comput. Surv., vol. 35, no. 3, pp. 268–308, Sept. 2003. [14] Dror G. Feitelson, “Workload modeling for performance evaluation,” in Performance Evaluation of Complex Systems: Techniques and Tools, Performance 2002, Tutorial Lectures, London, UK, UK, 2002, pp. 114–141, Springer-Verlag. [15] Zujie Ren, Jian Wan, Weisong Shi, Xianghua Xu, and Min Zhou, “Workload analysis, implications, and optimization on a production hadoop cluster: A case study on taobao,” Services Computing, IEEE Transactions on, vol. 7, no. 2, pp. 307–321, April 2014. [16] Andrew Burkimsher, Iain Bate, and LeandroSoares Indrusiak, “Scheduling hpc workflows for responsiveness and fairness with networking delays and inaccurate estimates of execution times,” in Euro-Par 2013 Parallel Processing, Felix Wolf, Bernd Mohr, and Dieter an Mey, Eds., vol. 8097 of Lecture Notes in Computer Science, pp. 126–137. Springer Berlin Heidelberg, 2013. [17] “NVIDIA CUDA TOOLKIT V6.5 (RN-06722001 v6.5),” http://docs.nvidia.com/cuda/pdf/CUDA_ Toolkit_Release_Notes.pdf, August 2014.