ISSN electrónico 2145-9371 ISSN impreso 0122-3461 vol. 31 n.° 1, Enero-Junio de 2013 Fecha de recepción: 13 de agosto de 2012 Fecha de aceptación: 27 de febrero de 2013

ARTÍCULO DE INVESTIGACIÓN / RESEARCH ARTICLE

Supply chain optimization with stochastic parameters by using the sample
average approximation method

John Willmer Escobar^*
Pontificia Universidad Javeriana, Cali (Colombia)

Juan José Bravo^**

Carlos Julio Vidal^***
Universidad del Valle, Cali (Colombia)

Correspondencia: John Willmer Escobar, Calle 18 No. 118-250 Cali (Colombia). Tel-Fax. +572 3218266.

Subvenciones y apoyos: El primer autor agradece el apoyo económico a través de una beca como asistente de docencia de la Universidad del Valle.

Resumen

Este artículo presenta el problema del diseño de una red de distribución de gran escala con parámetros estocásticos. La problemática central radica en la determinación de las decisiones de expansión o contracción de algunos de los eslabones de una red, teniendo en cuenta la variabilidad de la demanda. El problema se ha formulado como un modelo estocástico de dos etapas. Las decisiones de la primera etapa son de tipo estratégico; mientras que las decisiones de la segunda etapa son de tipo táctico. El modelo está basado en el caso de una compañía multinacional de alimentos que abastece a todo el territorio colombiano y varios mercados internacionales, dentro de los que se destacan Venezuela, Ecuador, Chile y algunos países de Centroamérica. La estrategia de solución adoptada es conocida como Sample Average Approximation (SAA). Dicha estrategia usa un esquema de aproximación por promedios muestrales para la solución de problemas estocásticos. Se presentan experimentos computacionales con diferentes tamaños de muestra. Los resultados obtenidos reflejan la importancia y eficiencia de la metodología propuesta como alternativa para el tratamiento de la variabilidad para redes de suministro de gran escala.

Palabras Clave: Diseño de Redes de Suministro, Logística, Programación Estocástica, Sample Average Approximation (SAA), Simulación Montecarlo.

Abstract

This paper presents a supply chain design problem with stochastic parameters for a large-scale company. The main problem consists to determine the decisions of expansion or contraction of some echelons by considering the variability of the demand. The problem is formulated as a two-stage stochastic model. The first-stage decisions are strategic, while the second-stage decisions are tactical. The model is based on a real-world case from amultinational food company, which suppliesthe Colombian territory and different international markets such as: Venezuela, Ecuador, Chile and some Central American Countries. The solution strategy adopted is known as Sample Average Approximation (SAA). This strategy uses an approximation scheme by sample averages for solving stochastic problems. Computational experiments with different sample sizes are presented. The results show the importance and efficiency of the proposed approach as analternative to the treatment of the variability for large-scale supply chains.

Keywords: Logistics, Montecarlo Simulation, Sample Average Approximation (SAA), Stochastic Programming, Supply Chain Design.

La logística comprende la integración de las actividades de abastecimiento, producción, transporte, distribución, inventario, almacenamiento, manipulación de materiales y empaque, así como el flujo de información entre ellas. Ballou [1] ha definido tres niveles de planeación logística dependiendo del horizonte de tiempo: nivel estratégico, nivel táctico y nivel operacional. El nivel estratégico considera el horizonte de tiempo más largo (usualmente más de un año). Este nivel requiere aproximación y agregación de datos. En este nivel generalmente se consideran decisiones de selección y asignación de proveedores y/o clientes, número, tamaño y localización de las instalaciones, tipo de productos a fabricar y/o distribuir y las decisiones de tercerización de alguna de las operaciones logísticas.

El nivel táctico es el intermedio en el horizonte de tiempo y requiere exactitud en los datos. Este nivel comprende decisiones de análisis del comportamiento de la demanda, la selección de técnicas de pronóstico, la administración de inventarios, la determinación de políticas de producción, almacenamiento y distribución y la selección del modo de transporte, entre otras. El nivel operacional involucra decisiones de corto plazo, frecuentemente en términos de días u horas, para las cuales se necesitan datos transaccionales. Generalmente este nivel comprende decisiones de planiicación de recursos, determinación de planes de emergencia, prioridades y asignación de picking, distribución de carga y ruteo de vehículos, entre otras. El problema considerado en este artículo involucra algunas decisiones de orden táctico y estratégico relacionadas con el diseño de la cadena de suministro.

La importancia del diseño de la cadena de suministro fue reconocida a principio de los años 60. Las primeras teorías sobre localización fueron propuestas por economistas agrarios y geógrafos [1]. El tema común de todos los trabajos presentados en dicha época fue el reconocimiento de los costos de transporte para determinar la localización de una instalación. Todos estos trabajos asumieron que los parámetros que influencian el diseño de la red son en general determinísticos. Esta idea a lo largo del tiempo ha evolucionado considerando la aleatoriedad en los parámetros que afectan el diseño de una cadena de suministro.

A través de la vida útil de una red de distribución, una compañía experimenta fluctuaciones en la demanda, precios y tasas de cambio en ambientes competitivos. En contextos de este tipo, la alta gerencia de las compañías debe considerar estrategias que le ayuden a tomar mejores decisiones del diseño de la red. Estos desafíos han incrementado el interés por el estudio de modelos de programación estocástica en la última década. Dichos modelos, consideran de alguna manera, la variabilidad presente en los elementos críticos de decisión (demanda, precios, costos, etc.).

Bajo incertidumbre no toda la información está disponible y algunos parámetros de una red de suministro deben ser modelados como variables aleatorias. En ese sentido la construcción de un modelo de programación entera mixta determinístico, acompañado de la estadística y la simulación, provee métodos eicientes para obtener soluciones al problema de diseño de redes de distribución de productos de consumo masivo con elementos estocásticos.

En el presente artículo se propone un modelo de diseño de redes de distribución de productos de consumo masivo estocástico de dos etapas. El modelo considera un sistema de distribución - inventario de tres eslabones en la cadena de suministro, en la cual se deben tomar decisiones de expansión o repliegue en centros de distribución, permitiendo tomar decisiones tácticas y estratégicas de la red logística. Además del artículo presentado por Santoso et al.[2], son escasos los artículos que incluyen la solución de un caso real de cadena de abastecimiento con consideraciones estocásticas aplicando la metodología SAA, razón por la cual la principal contribución de este artículo es examinar la aplicabilidad y efectividad de un modelo estocástico en un caso real de estudio asociado a una cadena de suministro de productos de consumo masivo, utilizando esta técnica de solución y considerando la variabilidad de la demanda en las zonas de mercado nacionales e internacionales.

Este artículo aborda el caso de una compañía multinacional de alimentos que abastece a todo el territorio colombiano y varios mercados internacionales, dentro de los que se destacan Venezuela, Ecuador, Chile y algunos países de Centroamérica. La empresa cuenta con aproximadamente 40 familias de productos terminados, que representan alrededor de 410 sku's (Stock Keeping Units). La compañía cuenta con más de 1.800 clientes nacionales y más de 2.000 clientes internacionales. Entre los centros de distribución y las plantas existe un flujo continuo de productos por medio del transporte terrestre con un operador logístico. Los centros de distribución tienen gran cantidad de producto en inventario cuyo contenido se subdivide en dos grupos básicos: a) productos para exportación, y b) productos para venta nacional.

En primera instancia se ha diseñado un modelo matemático de programación entera mixta determinístico (MILP) que permite representar la realidad del sistema bajo estudio. Seguidamente, se obtienen muestras de los parámetros aleatorios con los cuales se alimenta el modelo de programación entera mixta estocástico (SILP). De esta manera, se obtiene una solución para un escenario determinado con la cual se trata de validar el modelo mediante el análisis de respuestas y la experimentación. De acuerdo con los resultados de la validación, fue posible retornar al rediseño del modelo MILP y su reinamiento, convirtiéndose así en un proceso iterativo que evolucionaba de acuerdo con la dinámica del sistema estudiado.

Para resolver el modelo SILP se ha utilizado la técnica SAA desarrollada por [2] y Kleywegt et al. [3].Dicha metodología utiliza un esquema de aproximación mediante Simulación Montecarlo, a problemas de optimización de modelos de programación lineal entera mixta estocásticos.

En la sección 2 se presenta la revisión de la literatura del diseño de cadenas de suministro con elementos estocásticos. La sección 3 presenta el modelo determinístico (MILP) para la compañía bajo estudio. La estructura general del modelo estocástico (SILP) y su estrategia de solución se detalla en la sección 4. Finalmente, resultados computacionales y conclusiones son presentados en las secciones 5 y 6 respectivamente.

Dada la naturaleza integral de los modelos de diseño de redes de distribución, y de la incidencia que sus decisiones tienen a lo largo del sistema de abastecimiento, el campo de las áreas del conocimiento relacionadas con esta problemática es bastante amplio. Algunas de las áreas especíicas de estudio que se distinguen, entre otras, son: modelos de diseño de redes determinísticos de una sola instalación, modelos de diseño de redes determinísticos considerando múltiples instalaciones, modelos de diseño de redes utilizando métodos exactos, modelos de diseño de redes considerando métodos de simulación, modelos de diseño de redes utilizando métodos heurísticos, modelos dinámicos de diseño de redes y modelos de diseño de redes estocásticos.

Muchas de estas áreas, se han integrado para el estudio de la problemática por parte de expertos y académicos, logrando importantes desarrollos en la evolución de la teoría de diseño de redes de distribución [4] , [6], [17] , [18]. Igualmente diferentes funciones objetivo han sido consideradas para tratar numerosas aplicaciones [7], [8]. Existe una gran cantidad de aplicaciones exitosas a nivel nacional e internacional del diseño de cadenas de abastecimiento a nivel industrial; entre ellas podemos mencionar el caso Hewlett Packard propuesto por Laval et al. [9], caso IBM propuesto por Fleischmann et al. [10] y caso Elkem propuesto por Ulstein et al. [11], entre otros.

La literatura que considera modelos de diseño de redes estocásticos ha sido divida en dos grandes aproximaciones: aproximación probabilística y aproximación basada en escenarios. En ambos casos, cualquier tipo de parámetro del sistema es considerado como incierto. La aproximación probabilística considera explícitamente la distribución de probabilidad de las variables aleatorias modeladas, mientras que la aproximación basada en escenarios considera un conjunto de posibles valores futuros de la variable o parámetro en el análisis. En la práctica, la solución obtenida mediante la técnica de aproximación probabilística agrega un alto grado de complejidad computacional y matemático para la solución de problemas de diseño de redes de distribución, más aún cuando se consideran una gran cantidad de parámetros aleatorios con distribuciones de probabilidad continuas.

En la aproximación basada en escenarios, la incertidumbre es representada por un conjunto de escenarios discretos que capturan las diversas maneras en las que se puede modelar la variabilidad. Cada escenario tiene asociado un nivel de probabilidad que representa las expectativas del tomador de decisiones respecto a la ocurrencia de unos hechos en particular [12]. La metodología basada en escenarios en la práctica no es tan sencilla de implementar, debido a que es necesario predecir todas las posibles salidas de los parámetros o variables que se consideran aleatorios, además de que en algunos casos no se puede identiicar fácilmente un conjunto inito de escenarios discretos. Por ejemplo, suponga que se desea diseñar una red de distribución de 50 instalaciones (plantas, centros de distribución y zonas de mercado) con un parámetro aleatorio (demanda, tiempo de reposición, etc.). Asuma que este parámetro puede tomar un valor determinado de solamente tres posibilidades, y dicha cifra es independiente para cada una de las instalaciones. Entonces existe un total de 3⁵⁰ ≈7.18x10²³ escenarios para la representación de la incertidumbre, cifra que no podría ser enumerada fácilmente.

Modelos robustos de localización probabilísticos considerando la mayor cantidad de parámetros aleatorios posibles han sido estudiados por Chen et al. [13], Dasci y Laporte [14], Gabor y Van Ommeren [15] y Yu y Li [16]. Modelos de aproximación basada en escenarios han sido estudiados por Vidal y Goetschalckx [18].

Aun cuando las metodologías de aproximación probabilística y planeación por escenarios han sido ampliamente estudiadas por varios autores [13] , [16]; según [17], estas técnicas no son adecuadas para investigar y capturar en un camino sistemático las iteraciones de varios parámetros estocásticos de la cadena de abastecimiento. Para dar una solución a la problemática, los autores han desarrollado un método de descomposición estocástico de dos etapas, el cual explícitamente incorpora estas iteraciones. El método permite generar un número limitado de coniguraciones factibles. Para cada una de ellas se genera N réplicas calculando el valor esperado y la desviación estándar de los beneicios para cada diseño. Finalmente, se selecciona la mejor coniguración con la mejor combinación ponderada del valor esperado y la desviación estándar. Uno de los puntos interesantes de este trabajo es que ha demostrado que la coniguración de la cadena de abastecimiento óptima para el caso determinístico no es la mejor coniguración para el caso estocástico.

Tres de las publicaciones relevantes en el desarrollo de la teoría de optimización estocástica en cadenas de abastecimiento han sido propuestas por Kiya y Davoudpour [19], [2] y Snyder [20]. En la primera de ellas se considera el problema de rediseño de la red de distribución, donde se busca eliminar o integrar centros de distribución en una red doméstica, considerando aleatoriedad de la demanda y de los costos operacionales en un modelo de programación lineal entera mixta. En este trabajo se presenta un caso hipotético probando diferentes distribuciones de probabilidad de los parámetros estocásticos. El segundo artículo se concentra en el diseño de una red de suministro para una cadena de múltiples eslabones. En este trabajo, descomposición primal acelerada de Benders y SAA son utilizadas como técnicas de solución para el modelo estocástico. Resultados computacionales para un caso real son presentados. En [20] se propone un modelo de optimización que evidencia que la solución óptima de los modelos estocásticos es más robusta cuando se involucra toda la variabilidad de todos los parámetros en el modelo.

En Alonso-Ayuso et al. [21] se presenta un modelo de programación estocás-tica binaria de dos etapas. De igual manera, en [12] se propone un modelo de planeación de la demanda bajo incertidumbre utilizando programación estocástica. En [22], se muestra una aplicación interesante de la metodología SAA para el problema ULSP - Uncapacitated Lot Sizing Problem. Una aplicación interesante de la teoría de optimización estocástica a la industria papelera puede ser consultada en Santoso et al. [23]. Escobar [24] y Escbar et al. [25], describen de manera general la metodología utilizada para abordar el problema propuesto.

El desarrollo de la investigación se realizó en tres grandes etapas: formulación del modelo matemático de programación entera mixta determinístico (MILP), formulación del modelo de programación entera mixta estocástico(SILP), y descripción de la estrategia de solución SAA. A continuación se describe en detalle cada etapa de la metodología propuesta.

Formulación del Modelo Matemático Determinístico (MILP)

Características y supuestos generales del modelo

Las características generales y supuestos del modelo matemático determinístico de diseño de la red de distribución de productos de consumo masivo son los siguientes:

Notación Matemática

Conjuntos principales e inducidos

C = Conjunto de zonas geográficas de clientes, indexadas por k;

CD = Conjunto de centros de distribución, indexados por j;

PL = Conjunto de plantas de manufactura, indexadas por i;

PT = Conjuntos de productos terminados, indexados por p;

C(j) = Conjunto de zonas de mercado que pueden ser abastecidas por el centro de distribución j ∈ CD, C(j)⊆ C;

CD(k) = Conjunto de centros de distribución que pueden abastecer al cliente o zona de consumo k ∈ ;CD, cd(k) ⊆ CD;

CDEN(j) = Conjunto de centros de distribución que pueden recibir productos enviados desde el centro de distribución j ∈ cd, cden(j) ⊆cd;

CDRE(j) = Conjunto de centros de distribución que pueden enviar productos al centro de distribución j ∈ CD, CDRE(j)⊆cd;

CD(i) = Conjunto de centros de distribución que puede recibir productos desde la planta i ∈ PL, CD(i) ⊆ CD;

PL(p) = Conjunto de plantas que pueden fabricar el producto terminado p ∈ PT, PL(p) ⊆ PL;

PL(j) = Conjunto de plantas que pueden enviar producto terminado hacia el centro de distribución j ∈ CD, PL(j) ⊆ PL;

PT(i) = Conjunto de productos terminados que pueden ser fabricados en la planta i ∈ PL, PT(i) ⊆ PT.

Parámetros

CAPCD_j. = Capacidad del flujo a través del centro de distribución j ∈ CD para todos los productos p ∈ PT [kg/año];

CAPL_i. = Capacidad de producción en la planta i ∈ PL(p) para todos los productos p ∈ PT(i) [kg/año];

CF_j= Costo fijo de cierre del centro de distribución j ∈ CD [$/año];

CMANCD_jp = Costo de manipulación del producto terminado p ∈ PT en el centro de distribución j ∈ CD [$/kg. de p];

CP_ip = Costo variable de fabricación del producto terminado p ∈ PT en la planta . ∈ PL(p) [$/kg. de p];

CTRPTij = Costo de transporte de producto terminado p ∈ PT desde la planta i ∈ PL(p) hacia el centro de distribución j ∈ CD(i) [$/kg];

CTRPTCD_jj* = Costo de transporte de producto terminado p ∈ PT desde centro de distribución j ∈ CDRE(j*) hacia centro de distribución j* ∈ CDEN(j) [$/kg]; j ≠ j*

CTRPTC_jk = Costo de transporte de producto terminado p ∈ PT desde el centro de distribución j ∈ CD(k) hacia el cliente k ∈ C(j) [$/kg];

LTPjk= Tiempo promedio de reposición de productos desde centro de distribución j ∈ CD(k) hacia la zona de mercado k ∈ C(j) [días];

LTB_jj*= Tiempo promedio de reposición de productos desde centro de distribución j ∈ CDRE(j*) hacia centro de distribución j* ∈ CDEN(j), j ≠ j*[días];

CVLTP_jk = Coeficiente de variación del tiempo de reposición LTPjk;

CVLTP_jj* = Coeficiente de variación del tiempo de reposición LTBjj*;

DEM_kp = Demanda del producto terminado p ∈ PT en la zona de mercado k ∈ C [unidades de p/año];

FPT_jp = Factor de inventario de seguridad del producto p ∈ PT en centro de distribución j ∈ CD;

FPESO_p = Factor de peso por unidad de p ∈ PT [kg/unidad de p];

PCAP = Costo de penalización por cada unidad de la variable Capdip [$];

VIPT_jp = Costo de mantenimiento de inventario de producto terminado p ∈ PT en el centro de distribución j ∈ CD; [$/unidad de p].

Variables de decisión

x_ijp = Cantidad de producto terminado p ∈ PT(i) a fabricar en la planta i ∈ PL(p), i ∈ PL(j) enviado hacia el centro de distribución j e CD(i) [unidades de p/año];

u_jj*p = Cantidad de producto terminado p ∈ PT a enviar hacia centro de distribución j* ∈ CDEN(j) desde centro de distribución j ∈ CDRE(j*), j ≠ j*; [unidades de p/año];

y_jkp = Cantidad de producto terminado p ∈ PT enviado hacia las zonas de consumo k ∈ C(j) desde el centro de distribución j ∈CD(k)[unidades de p/ año];

w_j = Variable binaria asociada a cada centro de distribución j ∈ CD: "1" si se decide dejar abierto, "0" de lo contrario.

Capdip_i = Capacidad adicional de producción a considerar en la planta i ∈ PL para todos los productos terminados p PT(i)[kg/año];

Modelo matemático determinístico MILP

La función objetivo del modelo matemático MILP considera la minimización de la suma de costos de inventario de seguridad en CD's (1), costos de transporte de producto terminado de plantas hacia CD's (2), costos de transporte de producto terminado entre CD's y de CD's hacia clientes (3), costos variables de producción en plantas (4), costos de manejo de productos en CD's (5), costos fijos de cierre en CD's (6) y costos de penalización por excesos y capacidades adicionales en plantas (7):

Costos de inventario de seguridad en CD's

Costos de transporte de producto terminado de plantas hacia CD's

Costos de transporte de producto terminado entre CD's y de CD's hacia clientes

Costos variables de producción en plantas

Costos de manejo de productos en CD's

Costos fijos de cierre en CD's

Costos de penalización por excesos de capacidad

Las restricciones impuestas del modelo matemático son:

Note que el grupo de ecuaciones (8) se generan para cada planta del conjunto disponible, restringiendo la producción de una planta por encima de su capacidad. En el caso que fuera necesario utilizar capacidad extra, la variable Capdipl_i tomaría el valor de la capacidad extra. Esta variable es penalizada en la función objetivo con el parámetro PCAP. En el caso particular de la compañía bajo estudio se puede aumentar la capacidad de fabricación de las plantas en casos estrictamente necesarios ante el consumo del valor pre-establecido CAPL_i. La finalidad de dicha estrategia es presentar alternativas flexibles de aumento de capacidad de producción para abastecer plenamente la demanda y favorecer el flujo de inventarios a través del sistema para el total de productos p.

Las ecuaciones (9) y (10) restringen la capacidad en los centros de distribución (CD's). La modelación de la capacidadse establece de manera que el peso total despachado (y también recibido) de productos en un período determinado desde un centro de distribución no debe exceder su capacidad máxima. Las restricciones (11) determinan el balance de productos en los centros de distribución. El grupo de restricciones (12) asegura el cumplimiento de la demanda. Finalmente, el grupo de ecuaciones (13) están relacionadas con las restricciones de no negatividad.

Formulación del Modelo Matemático Estocástico (SILP)

Como se mencionó en la sección 2, varios trabajos de investigación y aplicación en la industria han abordado el problema de optimización de la cadena de abastecimiento mediante técnicas matemáticas de naturaleza determinística. Elementos comunes de estos trabajos, han sido el desarrollo de modelos de programación lineal entera mixta, modelos de programación

dinámicos y en algunas ocasiones modelos de programación no lineal que se pueden descomponer en un conjunto de modelos lineales. Todos los parámetros considerados en cada uno de estos trabajos se han asumido con condición de certeza.

Uno de los elementos que involucra sin lugar a dudas una gran fuente de incertidumbre a los problemas de diseño de redes es la demanda de productos para cada zona de consumo. Es evidente que se han desarrollado un "gran número" de teorías y metodologías alrededor de la solución a la problemática del pronóstico de demanda de una compañía. Sin embargo, a pesar del gran esfuerzo y los avances en la formulación de una considerable cantidad de modelos de toda índole, la variabilidad en la estimación de la demanda no ha sido eliminada por completo; esto debido a que existen factores de índole social, ambiental, cultural, etc., que no se pueden predecir con alguna exactitud matemática. El modelo estocástico propuesto considera dicha variabilidad.

Modelo matemático

Para efectos de la formulación del modelo estocástico SILP, es conveniente usar la siguiente notación compacta para el modelo (1) - (13):

Los vectores c, p, q, r,v, d y s corresponden a costos de cierre de centros de distribución, producción, transporte, inventario, penalización, demandas y capacidades de fabricación respectivamente.El grupo de restricciones (11) se puede reformular de la siguiente manera:

En la ecuación (16), las matrices N, Q y L representan la sumatoria de los lados izquierdos de la de la expresión (21). Por su parte la notación de las matrices D, F, R, y B corresponden a las sumatorias de los lados izquierdos de las expresiones (12),(8) y (9) o (10) respectivamente. Finalmente M representa el vector de capacidades de los centros de distribución.

Para convertir el modelo (14) - (20) en estocástico, se ha asumido que la demanda de cada ítem en cada zona de consumo es un parámetro aleatorio con distribución de probabilidad conocida. Para hallar cada una de las distribuciones de demanda, se hace necesario el uso de pruebas de bondad y ajuste estadístico. En el presente estudio se ha utilizado la técnica propuesta por Rici [26]. En este trabajo, el procedimiento consiste en encontrar una función matemática que represente en buena medida una variable aleatoria estadística. Se toma una serie de observaciones cuantitativas del parámetro x₁, x₂,... x_n y se le realizan diferentes pruebas a los datos para hallar la función de densidad de probabilidad f (x, q), donde q es un vector de los parámetros a estimar con los datos disponibles.

Para distinguir parámetros aleatorios de las realizaciones particulares de los modelos, se utilizará en particular el vector aleatorio ξ = (d). Los elementos de este vector son aleatorios y representan la componente probabilística del problema de optimización del diseño de la red de distribución. Considérese entonces la siguiente formalización del modelo estocástico de dos etapas SILP (22) - (29):

En el modelo estocástico de la primera etapa (22) - (23), c^Tes un vector de costos asociado a decisiones de primera etapa, que en este caso hace referencia al cierre o apertura de centros de distribución. En particular, el vector w es un vector de variables binarias relacionado con la toma de decisiones en esa primera etapa (repliegue de centros de distribución). En este caso el valor de |p| determina el valor del tamaño del problema en términos del número de variables de decisión abordadas en la primera etapa. Finalmente, Q(w,ξ) hace referencia al valor óptimo del problema de la siguiente etapa:

En el problema de la segunda etapa (24) - (29),x, y, u y z son vectores de variables de la segunda etapa, relacionadas con decisiones de cantidades de producto a fabricar y enviar entre nodos, y las capacidades adicionales a utilizar en planta. Los vectores p, q, r y v son los coeficientes de costo asociados a las variables representadas por los vectores x, u, y y z respectivamente. El conjunto de restricciones (25) - (27) están asociadas a las restricciones (16) - (18) una vez se han considerado las decisiones en la primera etapa. El grupo de restricciones (28) representa la relación entre las decisiones de la primera etapa (variables de decisión w) y las decisiones de la segunda etapa (variables de decisión x, y, u y z). En particular para nuestro caso, si la decisión de la primera etapa es determinar el cierre de un centro de distribución, el grupo de restricciones (28) impedirá que se asigne flujo de producto por ese centro de distribución.

De acuerdo con [19], existen dos fuentes potenciales de dificultad en la solución de modelos lineales estocásticos de dos etapas: i) la evaluación de la función objetivo f(w) para una configuración de centros de distribución determinada, supone el cálculo del valor esperado de la función lineal Q(w,ξ).En el caso de que los valores de demanda correspondan a distribuciones de probabilidad continuas, sería necesario calcular un número considerable de integrales para obtener el valor esperado de Q(w,ξ), lo que en la práctica es extremadamente difícil de resolver; ii) en el caso de que los valores de demanda correspondan a distribuciones de probabilidades discretas, para calcular el valor esperado de Q(w,ξ), sería necesario calcular un gran número de escenarios con las combinaciones de los parámetros. La cantidad de escenarios determinaría el número de subproblemas a resolver, que en nuestro caso en particular desde el punto de vista computacional sería complejo.

Descripción de la estrategia de solución SAA

El presente trabajo usa la metodología SAA para reducir la complejidad de solución del problema (22) - (29). En particular, se soluciona el problema SILP repetidas veces con un conjunto pequeño de escenarios. Para ello es necesario generar N muestras aleatorias del vector ξ, aproximando el valor de (22) mediante el promedio de los valores de la función objetivo obtenido en cada uno de los escenarios generados. De esta manera, el problema (22) - (29) es aproximado por el siguiente problema SAA:

La solución óptima del problema (30), Ŷ _N , y el valor óptimo, v _N , convergen con probabilidad de uno a la solución óptima del problema original (22) -(29) cuando el tamaño de la muestra se incrementa [3]. Asumiendo que se quiere resolver el problema SAA con una tolerancia (gap) de optimalidad δ> o, es posible estimar el tamaño de la muestra necesaria para garantizar una solución ε- optima al valor del problema original con una probabilidad de al menos 1 - α de acuerdo a la siguiente ecuación:

Donde ε > δy α Є (0,1). En (31), el término σ²_maxestá relacionado con la variabilidad de Q(w*ξⁿ) en la solución optima w* (para más detalles consultar [3]). En la práctica es usual tener en cuenta la relación entre la calidad de la solución óptima para el problema SAA y el esfuerzo computacional necesario para encontrarla. En efecto, solucionar el problema (30) con muestras independientes repetidamente puede ser más eficiente que incrementar el tamaño de muestras N. El procedimiento SAA descrito puede ser consultado en [2]

1) Generar M muestras independientes de tamaño N resolviendo el correspondiente problema SAA:

2) Calcular el promedio de todos los valores óptimos de la función objetivo

del problema SAA, y su varianza,

El promedio del valor objetivo estadísticamente provee un limite inferior para el valor objetivo del problema original (22) - (29). Para mayores detalles se puede observar Norkin et al [27] y Mak et al. [28].

3) Escoger una solución factible para el problema (22) - (29). Con dicha solución, estimar el valor objetivo del problema original usando una muestra aleatoria de tamaño N':

El estimador provee un límite superior en el valor objetivo del problema original. El valor de N' es generado de manera independiente a las muestras usadas en la solución de los problemas SAA. Debido a que la solución de la primera etapa es fija, se pueden seleccionar un gran número de escenarios para N'. La varianza de se puede estimar mediante la siguiente ecuación:

Calcular los estimadores de la brecha de la optimalidad y su varianza. Usando los valores calculados en los pasos 2 y 3, obtenemos dichos resultados mediante las siguientes ecuaciones:

Finalmente, en la metodología de solución propuesta, se han comparado las soluciones del modelo de programación estocástico (SILP) con el modelo de optimización de los valores promedios (MVP). Este modelo es obtenido reemplazando la distribución de la demanda en cada zona de consumo por sus valores promedio. La función objetivo del MVP está dada por la siguiente ecuación:

Note que la función objetivo del modelo MVP está compuesto por el vector de costos asociados a las decisiones de primera etapa y por el vector de costos asociadas a las decisiones de segunda etapa, considerando el valor promedio de la demanda para cada zona de consumo. Si los valores promedios de la variable randómica ξⁿ son fáciles de calcular, el problema de los valores medios es un problema deterministico que puede ser resuelto usando algoritmos de optimización deterministica logrando resolverlo en un tiempo computacional menor que un modelo estocástico.

El algoritmo SAA fue implementado en C++ y los experimentos fueron desarrollados en una Intel Core Duo CPU (2.00 GHz) bajo Linux Ubuntu 11.04 con 2 GB de memoria. Las distribuciones de demanda fueron obtenidas mediante el software Oracle Crystall Ball. Para resolver los diferentes subproblemas SAA se utilizó CPLEX 11.1.

En nuestro caso, se utilizaron cuatro diferentes tamaños de muestra para N =20 y 30,M = 20 y 30, yN'= 300. Los valores del tamaño de muestra han sido seleccionados verificando la convergencia de la función objetivo versus el tiempo computacional requerido para alcanzarlo.

En las Tablas 1 y 2 se resumen los resultados obtenidos al aplicar el procedimiento SAA sobre el caso de estudio. Nótese que se reportan experimentos computacionales para cada uno de los tamaños de muestra utilizados (N, M, N'). Para el primer caso (20, 20, 300), el algoritmo de optimización debe resolver, 400 subproblemas de programación lineal entera mixta en la primera etapa y 300 problemas de programación lineal en la segunda. La cantidad de subproblemas de la primera etapa se aumentan a 600 modelos para los casos (30, 20, 300) y (20, 30, 300), y de 900 para el caso (30, 30, 300).

La Tabla 1 muestra los valores obtenidos en los pasos 2 y 3 de la metodologia SAA:

La Tabla 2 muestra los valores de brecha de optimalidad de la metodología SAA:

La Tabla 3 muestra la comparación de la solución del MVP y el SILP en términos del valor objetivo promedio y su respectiva variabilidad. Los resultados son reportados para un tamaño de N = 30. Los resultados de la Tabla 3 muestran que el valor promedio y la variabilidad del modelo estocástico SILP es mucho menor a la presentada por el problema MVP; y en este caso la brecha de optimalidad, indica la superioridad de la solución del SAA. Note que al considerar tamaños pequeños de muestra, el modelo SILP obtiene soluciones cercanas al óptimo (gap < 1%).

Finalmente, es importante mencionar que los tiempos de ejecución de la metodología SAA no son reportados en este trabajo, debido a que son despreciables en relación al horizonte de tiempo de la decisión que se está abordando (orden estratégico-táctico).

Los resultados computacionales de la implementación del SAA reflejan la importancia y efectividad de la metodología propuesta como alternativa para el diseño de redes de suministro bajo condiciones de variabilidad de demanda. Los resultados obtenidos permitieron a la compañía: i) Reducir en aproximadamente 30% los costos de operación e infraestructura, ii) Disminuir la complejidad de la administración en la operación y la asignación de prioridades en el proceso de distribución, iii) Mejorar la ubicación del producto en inventario, iv) Mejorar el conflicto entre divisiones, regiones y gerencias. De igual forma cabe anotar que los resultados obtenidos por el modelo matemático difieren notoriamente de las decisiones de cierre y consolidación planeadas por la alta gerencia.

De manera general, se puede mencionar que la estrategia algorítmica SAA como metodología de solución de problemas reales de gran escala (problemas que incluyen miles de variables y restricciones) es bastante promisoria, debido a que se cuentan con diversos algoritmos, técnicas de descomposición, etc., que son eficientes para resolver problemas determinísticos.

Los resultados obtenidos en la investigación son importantes como precedente para el abordaje y solución de problemas reales de diseño de redes de abastecimiento de gran escala, en los cuales la consideración de la variabilidad y la incertidumbre son una preocupación cada vez mayor.

Queda abierto un interesante campo de investigación, en el cual se pueden realizar extensiones en cuanto a la optimización de redes de distribución con diversas medidas de desempeño como son el Valor Presente Neto (vpn) o el Flujo de Caja Descontado (fcd), consideraciones de modo de transporte y precios de transferencia para cadenas de suministro globales, formulación de modelos dinámicos estocásticos multiperiódicos para la solución de problemas reales y la revisión de las bondades de otras metodologías de solución para modelos lineales estocásticos.

Agradecimientos

El primer autor agradece al Ministerio de Educación e Investigación Italiano (MIUR), Italia y los profesores del Departamento de Electrónica, Informática y Sistema (DEIS) de la Universidad de Bologna.

^* Departamento de Ingeniería Civil e Industrial, Pontificia Universidad Javeriana, Cali (Colombia). johnwillmer.escobar2@unibo.it

^** Escuela de Ingeniería Industrial, Universidad del Valle, Cali (Colombia). juan.bravo@correounivalle.edu.co

^*** Escuela de Ingeniería Industrial, Universidad del Valle, Cali (Colombia). carlos.vidal@correounivalle.edu.co

[1] R.H. Ballou, Business logistics management: planning, organizing, and controlling the supply chain, 4a ed., New York: Prentice Hall Upper Saddle River, 2004, pp. 38-78.

[2] T.Santoso, et al., "A stochastic programming approach for supply chain network design under uncertainty". European Journal of Operational Research, vol. 167(1), pp.96-115, 2005.

[3] A.J.Kleywegt, et al., "The sample average approximation method for stochastic discrete optimization", SIAM Journal on Optimization, vol. 12(2), pp.479 -502, 2001.

[4] J.Current, et al., "Dynamic facility location when the total number of facilities is uncertain: A decision analysis approach". European Journal of Operational Research, vol. 110(3), pp.597-609, 1998.

[5] C.J.Vidal, and M.Goetschalckx,"A global supply chain model with transfer pricing and transportation cost allocation". European Journal of OperationalResearch, vol. 129(1), pp.134-158, 2001.

[6] O.Berman, and Z.Drezner,"A probabilistic one-centre location problem on a network". Journal of the Operational Research Society, vol. 54(8), pp.871-877, 2003.

[7] S.H.Owen, and M.S.Daskin,"Strategic facility location: A review". European Journal of Operational Research, vol. 111(3), pp.423- 447, 1998.

[8] B.M.Beamon, "Measuring supply chain performance". International Journal of Operations & Production Management, vol. 19(3), pp.275-292, 1999.

[9] C.Laval, et al.,"Hewlett-Packard combined OR and expert knowledge to design its supply chains". Interfaces, vol. 35(3), pp.238-247, 2005.

[10] M.Fleischmann, et al."Tntegrating closed-loop supply chains and spare-parts management at IBM". Interfaces, vol. 33(6), pp.44-56, 2003

[11] N.L. Ulstein, et al., "Elkem uses optimization in redesigning its supply chain". Interfaces, vol. 36(4), pp.314-325, 2006.

[12] A.Gupta, and C.D.Maranas, "Managing demand uncertainty in supply chain planning". Computers & Chemical Engineering, vol. 27(8-9), pp. 1219-1227, 2003.

[13] G.Chen, et al.,"A New Model for Stochastic Facility Location Modeling". Research Report, Department of Industrial and Systems Engineering, University of Florida, 2005.

[14] A.Dasci, and G.Laporte. "An analytical approach to the facility location and capacity acquisition problem under demand uncertainty". Journal of the Operational Research Society, vol. 56(4), pp. 397- 405, 2004.

[15] A.F.Gabor, and J.C.W.Van Ommeren. "An approximation algorithm for a facility location problem with stochastic demands and inventories". Operations research letters, vol. 34(3), pp. 257-263, 2006.

[16] C.S.Yu, and H.L.Li,"A robust optimization model for stochastic logistic problems". International Journal of Production Economics, vol. 64(1), pp.385-397, 2000.

[17] M.Goetschalckx, et al.,"Designing flexible and robust supply chains". Presented at IERC, Dallas Spring, 2001.

[18] C.J.Vidal, andM. Goetschalckx."Modeling the effect of uncertainties on global logistics systems". Journal of Business Logistics, vol. 21(1), pp. 95-120, 2000.

[19] F. Kiya, and H Davoudpour, "Stochastic programming approach to redesigning a warehouse network under uncertainty", Transportation Research Part E: Logistics and Transportation Review, vol. 48(5), pp. 919-936, 2012.

[20] L.V.Snyder, "Supply chain robustness and reliability: Models and algorithms", PhD dissertation, Evanston, Illinois, Northwestern University, 2003.

[21] A. Alonso-Ayuso, et al."An approach for strategic supply chain planning under uncertainty based on stochastic 0-1 programming". Journal of Global Optimization, vol. 26(1), pp.97-124, 2003.

[22] H.H.Toro."SAA-Sample Average Approximation Method, aplicado a la solución de modelos de Programación Lineal". Heurística, vol. 14, pp. 10-20,2007.

[23]T.Santoso, et al."Strategic design of robust global supply chains: two case studies from the paper industry". Presented at TAPPI Conference, Atlanta, 2004.

[24] J.W.Escobar, Modelación y optimización de redes de distribución de productos de consumo masivo con elementos estocásticos". Proceedings of XIV Latin American Summer Workshop on Operations Research (ELAVIO), El Fuerte, Mexico. 2009.

[25] J.W. Escobar, et al.,"Optimización de redes de distribución de productos de consumo masivo en condiciones de riesgo". Proceedings of XXXIII Congreso Nacional de Estadística e Investigación Operativa (SEIO), Madrid, Spain. May 2012.

[26] V.Rici (2010, March 12). Fitting Distributions with R [Online]. Available: http://www.fsf.org/.

[27] V.I.Norkin, et al."A branch and bound method for stochastic global optimization". Mathematical programming, vol. 83(3), pp.425-450, 1998.

[28] W.K.Mak, et al.,"Monte Carlo bounding techniques for determining solution quality in stochastic programs". Operations Research Letters, vol. 24(1), pp. 47-56, 1999.