Metodología Matemática-Algorítmica de Programación de Operaciones Aplicada a un Caso de Estudio de Escala Industrial

Algorithmic-Mathematical Scheduling Methodology Applied to an Industrial Size Test Case

Presentación: 25/01/2020

Aprobación: 02/11/2020

Santiago Zuffiaurre

Universidad Tecnológica Nacional, Facultad Regional Santa Fe - Lavaisse 610, (3000) Santa Fe - Argentina

santiago.zuffiaurre@gmail.com

Pablo A. Marchetti

Instituto de Desarrollo Tecnológico para la Industria Química (INTEC), Universidad Nacional del Litoral (UNL) - Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) - Güemes 3450, (3000) Santa Fe - Argentina

Universidad Tecnológica Nacional, Facultad Regional Santa Fe - Lavaisse 610, (3000) Santa Fe - Argentina

pmarchet@intec.unl.edu.ar

Resumen

En este trabajo se presenta una metodología para la programación de operaciones de procesos “batch” en instalaciones multiproducto multietapa. El modelo matemático empleado es de tipo mixto-entero lineal (MILP) y utiliza una representación de ranuras de tiempo (“time slots”). El modelo se complementa con un algoritmo iterativo, basado en la resolución de una secuencia de subproblemas, que permite identificar y fijar la programación de la etapa cuello de botella en cada paso. La metodología propuesta apunta a obtener soluciones de buena calidad para problemas de escala industrial en tiempos de cómputo razonables. Fue aplicada a un caso de estudio real de la industria farmacéutica, que involucra la programación de 30 productos en una planta de 6 etapas y 17 equipos. Si bien no garantiza la optimalidad de la solución hallada, a diferencia de otros aportes de tipo heurístico provee una cota inferior rigurosa que permite medir la calidad de la solución.

Palabras claves: Programación de operaciones, Planta multiproducto, Optimización.

Abstract

This paper presents a batch process scheduling methodology for multiproduct multistage facilities. A mixed-integer linear programming (MILP) mathematical model based on a time-slot representation is used. The model is complemented with an iterative algorithm that solves a sequence of subproblems, through which the bottleneck stage is identified and fixed at each step. The proposed methodology seeks to attain good quality solutions for industrial size problems in reasonable computational times. It has been applied to a real case study from the pharmaceutical industry, comprising the scheduling of 30 products in a plant with 17 units and 6 processing stages. Even though the proposed method does not guarantee the optimality of the best solution found, in contrast to other heuristic approaches it provides a rigorous lower bound from which its quality can be determined.

Keywords: Scheduling, Multiproduct facility, Optimization.

Introducción

Entre las actividades de planificación de la producción, la programación de operaciones en piso de planta (o “scheduling”) es fundamental para garantizar la competitividad de las organizaciones. De hecho, la rentabilidad de una planta industrial está estrechamente vinculada a la eficiente utilización de los recursos de producción disponibles. Esta actividad se caracteriza por un horizonte temporal de corto plazo, realizándose con una frecuencia semanal o quincenal.

En las últimas décadas se han realizado importantes esfuerzos desde las disciplinas de Ingeniería de Procesos e Investigación de Operaciones para el desarrollo de algoritmos y métodos de optimización para el problema de “scheduling”. Dado que permite garantizar la optimalidad de la solución cuando el problema planteado es resuelto satisfactoriamente, la formulación de modelos matemáticos de tipo mixto-entero lineal (MILP) es una de las metodologías más utilizadas. Las formulaciones de este tipo de problemas se pueden clasificar, según su representación de la variable tiempo, en aquellas de tiempo discreto y de tiempo continuo (Méndez et al., 2006).

Las representaciones de tiempo discreto se caracterizan por dividir el horizonte de programación en un número finito de intervalos de duración predeterminada, y restringir los tiempos de inicio o finalización de tareas a los extremos de estos períodos. Generan representaciones matemáticas más simples que otras formulaciones, al costo de perder exactitud si la grilla de tiempo no tiene la granularidad suficiente, y disminuir su desempeño computacional cuando el número de intervalos de tiempo aumenta (Kondili et al., 1993; Rodrigues et al., 2000).

Por otra parte, las representaciones de tiempo continuo asocian el tiempo a variables continuas, lo que permite modelar los instantes exactos en los que suceden los eventos. Sin embargo, estas formulaciones muchas veces requieren el uso de restricciones más complejas e incluso restricciones de tipo “big-M”. Las metodologías para procesos secuenciales utilizan representaciones continuas de la variable tiempo tales como: ranuras de tiempo o “time-slots” (Pinto y Grossmann, 1995; Lim y Karimi, 2003) y variantes del concepto de precedencia, incluyendo precedencia inmediata (Cerdá et al., 1997; Gupta y Karimi, 2003) y general (Méndez et al., 2001; Marchetti y Cerdá, 2009; Marchetti et al., 2012).

Harjunkoski et al. (2014) presentan una extensa revisión donde analizan las fortalezas y debilidades de los diferentes modelos de programación de operaciones colocando el foco en la implementación de aplicaciones industriales concretas. Dichos autores sostienen que las herramientas de optimización actuales pueden aplicarse satisfactoriamente para aumentar la eficiencia de la producción en piso de planta.

En este trabajo se desarrolla una metodología matemática-algorítmica para la programación de operaciones de una planta “batch” multiproducto multietapa, utilizando una representación basada en ranuras de tiempo (“time-slots”). El modelo se complementa con un algoritmo iterativo, basado en la resolución de una secuencia de subproblemas, que permite identificar y fijar la programación de la etapa cuello de botella en cada paso. Una característica fundamental del método propuesto consiste en su capacidad de resolver problemas reales de escala industrial, y en la consideración de una cota rigurosa que permite evaluar la calidad de la solución obtenida. La metodología fue aplicada a un caso de estudio real de la industria farmacéutica, introducido y estudiado previamente por Castro et al. (2009) y Kopanos et al. (2010).

La estructura del trabajo es la siguiente: en la sección Metodología se formaliza el problema abordado, se describe la formulación matemática basada en ranuras de tiempo utilizada, y se presenta en forma resumida la metodología iterativa propuesta en sus dos niveles. Luego, en la sección Resultados se introduce el caso de estudio tomado como base para este trabajo, y se presenta la solución obtenida junto con el requerimiento computacional requerido para su obtención. En la sección Discusión se analizan los resultados y presentan comparaciones y, para finalizar, se presentan las Conclusiones obtenidas.

Metodología

Descripción del problema

Dados: (i) una planta multiproducto multietapa de producción por lotes (planta “batch”) con numerosos equipos operando en paralelo en cada etapa del proceso, (ii) los lotes o “batches” a ser procesados y la secuencia de etapas requeridas para cada lote, (iii) el subconjunto de equipos disponibles para cada tarea (i, l), (iv) los tiempos de procesamiento ptij del “batch” i en el equipo j, (v) los tiempos de “setup” (o alistamiento) τi' ij dependientes de la secuencia (es decir, el tiempo de limpieza o preparación del equipo j cuando el lote i' se procesa inmediatamente antes del lote i), y (vi) la duración H del horizonte de tiempo disponible para la programación. El objetivo del problema es encontrar un programa de operaciones factible, que contemple la asignación, secuenciación y temporización detallada de las tareas, completando todos los lotes dentro del horizonte temporal especificado y, simultáneamente, minimizando el tiempo de finalización máximo o “makespan.”

En este trabajo no se consideran fechas de entrega para los lotes y se asume que existe disponibilidad ilimitada o suficiente de otros recursos de fabricación (por ejemplo: tanques de almacenamiento para productos intermedios y finales, mano de obra, servicios, etc.)

Formulación de ranuras de tiempo para la programación de operaciones multietapa

El foco de estudio del presente trabajo es un caso real y frecuente de la industria, aquel de una planta multiproducto en el que la cantidad de productos es ampliamente mayor a la cantidad de equipos en cada etapa. Se eligió representarlo mediante una formulación de ranuras de tiempo debido a su capacidad para manejar de manera eficiente la secuenciación de tareas consecutivas en un equipo dado.

El modelo de ranuras de tiempo (“time slots”) especifica un número fijo de períodos de tiempo para cada equipo, durante los cuales los lotes pueden ser procesados. Como formulación de base se tomó aquella propuesta por Pinto y Grossmann (1995), con modificaciones para representar los tiempos de transición dependientes de la secuencia. Esta representación considera intervalos de tiempo con duración variable, de modo que los instantes de comienzo y finalización de las tareas son modelados con variables continuas. A continuación se describe la formulación.

La decisión de asignar un determinado lote i a una ranura k asociada a un equipo j de una etapa l, se representa mediante la variable binaria Wijkl , la cual tomará el valor uno cuando se realiza la asignación; o cero en caso contrario. La ecuación (1) indica que cada lote debe ser asignado a una ranura de alguno de los equipos, para cada una de las etapas l que le correspondan. La ecuación (2) limita la cantidad de lotes que se procesan por ranura a 1. Seguidamente, la ecuación (3) asegura que las ranuras sean utilizadas en orden, evitando que se dejen ranuras sin utilizar entre otras asignadas (esto evita soluciones simétricas).

(1)

(2)

(3)

Las ecuaciones (4) y (5) sincronizan los tiempos de comienzo para el lote y la ranura asociados (y , respectivamente) cuando Wijkl = 1. Dado que se utilizan restricciones de tipo “big-M”, esta sincronización entre los tiempos de comienzo de “batches” y ranuras añade dificultad a la resolución del problema. En cambio, la relación entre los tiempos de comienzo mencionados y los tiempos de finalización (y ) se obtiene en forma directa considerando el tiempo de procesamiento ptij asignado al lote o ranura, respectivamente, en las ecuaciones (6) y (7).

(4)

(5)

(6)

(7)

 

Los tiempos de transición entre “batches” consecutivos procesados en el mismo equipo se introducen en la ecuación (8), donde COjk es una variable continua positiva que representa el tiempo de alistamiento efectivamente necesario entre los lotes asignados a las ranuras (k-1) y k. Su valor se determina en la ecuación (9), como la diferencia entre el máximo tiempo de alistamiento posible para el lote i en el equipo j , y la magnitud que corresponde restar según el lote precedente i'. Nótese que el lado izquierdo de la ecuación (9) sólo tendrá un valor positivo cuando Wijkl=1, en cuyo caso el “changeover” efectivamente requerido será el límite inferior de la variable COjk.

(8)

(9)

 

La ecuación (10) determina que el tiempo de comienzo de un “batch” en una etapa debe ser mayor o igual al tiempo de finalización en la etapa anterior. Finalmente, las ecuaciones (11) y (12) relacionan los tiempos de finalización de los lotes y los equipos como límite inferior del “makespan,” representado por la variable continua MK, la cual se minimiza en la función objetivo.

(10)

(11)

(12)

A fin de obtener una región factible más compacta y mejorar la velocidad de convergencia del resolvedor en la búsqueda de la solución óptima, se incorporan al modelo algunas restricciones de ajuste (“tightening constraints”). Se describen a continuación dos de estas restricciones, basadas en parámetros adicionales que representan el tiempo de comienzo más temprano (ESTil) y el tiempo remanente mínimo (MRTil) para la tarea (i, l). Estos parámetros son calculados inicialmente mediante la ecuación (13). Como se verá más adelante, a medida que avanza en las iteraciones el algoritmo propuesto actualiza los valores de estos parámetros para incluir información de las variables de decisión fijadas. Las ecuaciones (14) y (15) utilizan los parámetros ESTil y MRTil para establecer límites inferiores y superiores para el tiempo de comienzo y finalización, respectivamente, de cada ranura (j, k). En el caso de la ecuación (15), permite estimar un tiempo máximo luego del cual un lote no puede ser finalizado a fin de asegurar el cumplimiento del “makespan.” Otras ecuaciones de ajuste relevantes son aquellas diseñadas para proveer de cotas inferiores a los tiempos de transición, para lo cual se establecieron parámetros cuyo valor representa el mínimo tiempo de alistamiento para cada “batch.”

(13)

(14)

(15)

Algoritmo iterativo propuesto basado en la identificación de cuellos de botella

La aplicación a casos reales de la formulación matemática presentada se ve limitada a problemas de tamaño mediano/pequeño, ya que sólo allí el tiempo de cómputo será razonable. Para abordar la complejidad resultante de la combinatoria en problemas más grandes, se propone una metodología iterativa basada en los cuellos de botella. Comprende una iteración principal (main iteration) y otra interna (inner iteration).

En términos generales, la iteración principal selecciona en cada paso una etapa cuello de botella o crítica, fijando las decisiones 0-1 para dicha etapa. El proceso de selección consiste en resolver subproblemas (iteración interna) para cada etapa pendiente a fin de identificar aquella en la que se produce un cuello de botella, obteniendo una solución óptima o de buena calidad para dicha etapa. En base a esta solución se fijan las variables de decisión Wijkl correspondientes (en un vector solución x). Luego, el proceso se repite incorporando la información de las etapas previamente fijadas, de manera que en sucesivas iteraciones la cantidad de decisiones pendientes sea menor.

A continuación, se presenta en pseudocódigo el algoritmo principal utilizado para resolver el problema:

MAIN ITERATION

El algoritmo comienza fijando valores iniciales para z y zCUT, parámetros que representan el mayor valor mínimo de la función objetivo y un límite inferior para interrumpir la búsqueda, respectivamente. También, se inicializa la lista de etapas aún no fijadas Lpending con todas las etapas, preordenadas decrecientemente en base a una estimación del tiempo de procesamiento requerido. Seguidamente, mientras la lista mencionada no esté vacía se repite la iteración principal. Esta consiste en los siguientes pasos: Se actualizan las estimaciones de los parámetros ESTil (earliest start time) y MRTil (minimum remaining time). Se resetea la variable lsel, que especifica la siguiente etapa a fijar. Luego, para cada l pendiente, se resuelve la iteración interna (inner iteration) y procesa su resultado. Si se encuentra una solución factible del subproblema donde la función objetivo (FO) es menor que zCUT, el “solver” es interrumpido y se continúa analizando otra etapa. Caso contrario, y cuando no hubiera ninguna etapa seleccionada o el valor de la FO obtenido en la iteración interna zl fuera mayor que z, se define la etapa considerada como la siguiente a fijar, y se asigna z=zl. Seguidamente, se asigna a zCUT el valor máximo entre el valor de corte actual y el “lower bound” obtenido para la solución de la etapa l considerada. Una vez inspeccionadas todas las etapas pendientes, lsel será la siguiente etapa a fijar, es decir aquella con el mayor mínimo valor de la FO (“makespan”). La solución de esta etapa se incorpora el vector solución x, y la etapa lsel es quitada del conjunto Lpending. Mediante este proceso, en cada iteración del ciclo While se fijan las decisiones correspondientes a las variables binarias de una de las etapas, construyendo gradualmente una solución completa del problema.

Por su parte, la iteración interna resuelve una secuencia de modelos que representan el subproblema obtenido al considerar una única etapa l, y tiene por objeto encontrar una solución de la mejor calidad posible para dicha etapa. Para esto, se exploran diferentes configuraciones para el número de ranuras en los equipos de l. En esta iteración se resuelven en forma alternada un subproblema aproximado, para determinar un número candidato de ranuras posible, y un subproblema exacto, para encontrar la mejor solución para dicho número de ranuras. En los subproblemas aproximados se considera el último “slot” de cada equipo como una ranura “comodín.” Este tipo de ranuras admitirá la asignación de múltiples lotes, no considerándose las decisiones de secuenciación ni tiempos de comienzo y finalización asociados. Por lo tanto, para la ranura la ecuación (2) se omite y las demás restricciones se modifican consistentemente. Como no considera todas las decisiones de la etapa l, la solución del subproblema aproximado será incompleta. Sin embargo, el número de ranuras seleccionadas en dicha solución se utiliza como configuración de ranuras candidata para el subproblema exacto posterior. En este último, todas las ranuras de la etapa l serán utilizadas, por lo tanto la ecuación (2) se convierte en igualdad y la restricción (3) se omite. En la siguiente iteración interna, la última configuración de ranuras candidata se descarta mediante la introducción de restricciones de corte de tipo entero. A medida que se avanza en las iteraciones se registra la mejor solución encontrada para un subproblema exacto, que constituye la mejor solución factible para el subproblema de la etapa l.

En resumen, la estrategia planteada permite analizar diferentes configuraciones para el número de ranuras de los equipos en l, encontrando una solución exacta xl y manteniendo información rigurosa de las cotas durante el proceso. La variable zl es la función objetivo obtenida por la iteración interna al resolver el subproblema de la etapa l, mientras que constituye la mejor solución posible (o lower bound) de dicho problema. El vector xl representa las variables de decisión, e incluye las variables Wijkl de asignación de lotes a ranuras. El proceso termina cuando la diferencia entre las cotas es suficientemente pequeña o cuando se alcanza un número máximo de iteraciones. Si algún subproblema exacto obtiene una solución entera menor al parámetro zCUT la iteración interna termina anticipadamente. En este caso, la etapa l no será considerada crítica (es decir, cuello de botella) porque hay otra etapa cuya mejor solución hallada es mayor.

Utilizando la información de las variables de asignación fijadas es posible actualizar las estimaciones de ESTil y MRTil en la iteración principal. En las ecuaciones (16) y (17), cuando los índices js, ip e in representan el equipo seleccionado, el “batch” previo y el “batch” siguiente en dicho equipo, respectivamente, para el lote i en la etapa l. Los parámetros con índices (l - 1) o (l + 1) no definidos toman el valor cero, y se define en la ecuación (18). Nótese que la ecuación (13) es la estimación inicial de estos parámetros cuando ninguna de las etapas ha sido fijada.

(16)

(17)

(18)

Resultados

La metodología propuesta se aplicó a un caso de estudio real de la industria farmacéutica, estudiado previamente por Castro et al. (2009) y Kopanos et al. (2010). El ejemplo considerado involucra la producción de 30 “batches”. La planta de producción tiene 6 etapas y 17 equipos en total, pertenecientes los dos primeros a la etapa inicial y los restantes en grupos de tres a las etapas sucesivas. Algunos productos no requieren su procesamiento en la etapa 3, resultando un total de 162 tareas de procesamiento a ser programadas. El problema involucra tiempos de “changeover” dependientes de la secuencia de procesamiento, que incluso en algunas etapas son superiores a los tiempos de procesamiento de las tareas. Los datos detallados del problema pueden encontrarse en el material suplementario de Kopanos et al. (2010). A diferencia del problema original, en este estudio no se consideran los “due dates” introducidos por dichos autores.

El algoritmo y formulación matemática se implementaron con la herramienta de modelado GAMS 25.0, utilizando para hallar la solución de los modelos MILP el resolvedor o “solver” GUROBI 7.5. Este resolvedor fue seleccionado por su mejor desempeño y porque dispone del parámetro bestobjstop, el cual permite detener la búsqueda cuando se ha encontrado una solución de cierta calidad (zCUT). Se utilizó una PC genérica con procesador Intel Core i7 3.2 GHz y 16 GB de RAM, y se consideró un tiempo límite de 5 minutos para la resolución de los subproblemas de la iteración interna. Al ejecutar el algoritmo principal, fueron seleccionadas sucesivamente las etapas E4, E6, E2, E3, E5 y E1, hasta que el conjunto Lpending estuvo vacío. En cuanto al tamaño del modelo, éste presentó variaciones en las distintas iteraciones del algoritmo entre (1556, 716, 8051) y (423, 581, 2360) para el número de variables binarias, variables continuas y restricciones, respectivamente. En general, los subproblemas de las primeras iteraciones son los de mayor tamaño. La Figura 1 muestra el diagrama Gantt de la mejor solución obtenida. El tiempo de computación total de la metodología propuesta para alcanzar esta solución fue de aproximadamente 2 horas y 50 min.

Discusión

La solución obtenida y reportada en la Figura 1 presenta un “makespan” de 26,31 horas, obteniéndose con la metodología propuesta una cota inferior rigurosa de 21,94 h. Estos valores permiten determinar que el problema planteado aún admite una solución un 16,6% mejor, lo cual es el indicador de calidad de la solución (“relative gap”). Cabe mencionar que esta solución es menor a aquella reportada por Kopanos et al. (2010), con un valor de “makespan” de 26,56 h. Sin embargo, ambas soluciones no son directamente comparables porque Kopanos et al. (2010) incorporan fechas de entrega para cada lote, característica que requiere modificaciones del método que serán abordadas a futuro.

El caso de estudio planteado es de una dificultad considerable debido a la explosión combinatoria que generan las múltiples alternativas de asignación de tareas a equipos y de secuenciación de tareas asignadas a un mismo equipo. Para ilustrar esta dificultad, se resolvió una formulación más cercana al problema completo, utilizando también ranuras de tiempo. En dicha formulación se incluyeron todas las variables de decisión y ecuaciones para todas las etapas y equipos pero, para simplificar la elección del número de ranuras por equipo, se seleccionó la configuración de “slots” que aparece en la Figura 1. Por lo tanto, aún considerando todas las etapas, al prefijar el número de ranuras por equipo no se está considerando el espacio de soluciones completo del problema planteado. Este modelo con todas las etapas requirió 4066 variables binarias, 502 variables continuas, y 16479 restricciones, y permitió obtener como mejor solución un “makespan” de 33,57 h luego de 5 h de tiempo de CPU (con gap relativo 30%). En contraste, la metodología matemática-algorítmica propuesta permitió obtener una solución de mejor calidad debido al proceso de construcción incremental de la solución planteado, basado en resolver una secuencia de subproblemas enfocados cada uno en una etapa.