En mi primera publicación de TDS, analicé un problema del mundo real en un programa lineal entero. En el segundo, hice que ese programa fuera sólido contra la incertidumbre. En el tercero, introduje la programación estocástica: cuatro formas basadas en principios de incluir incertidumbre en el modelo en lugar de descartarla con la mano.
El tercer post terminó con una promesa. Dije que los modelos de recursos de dos etapas son fáciles de escribir y ordenar en teoría, pero en la práctica tienen un problema embarazoso: a medida que se agregan más escenarios, el llamado “equivalente determinista” explota en tamaño y el solucionador comienza a quejarse. Mencioné la descomposición de Benders como la solución estándar y luego desaparecí sin mostrar cómo funciona realmente.
Esta es la publicación que cumple esa promesa.
Comenzaremos donde quedó la publicación anterior, mostraremos por qué el enfoque obvio se desmorona a escala, resolveremos la única observación que nos salva y luego recorreremos lentamente el algoritmo en sí. Hay un poco más de matemáticas que la última vez, pero nada exótico; La herramienta más avanzada que usaremos es la dualidad LP, e incluso eso viene empaquetado como "podemos escribir el mismo problema de una manera diferente".
Resumen rápido: el modelo que estamos intentando resolver
En la publicación anterior, nuestro héroe era una empresa de moda en Alemania que se abastecía de ropa de invierno de Bangladesh. La demanda ξ es una variable aleatoria; La producción x debe decidirse antes de conocer ξ. Una vez que se revela la demanda, se toma una decisión de recurso para cubrir cualquier déficit a un costo mayor (un lote de emergencia de Rumania, por ejemplo).
Si asumimos una distribución discreta con S escenarios (valores ξ¹, ξ²,…, ξS^S con probabilidades p₁,…, pS), el modelo de recurso de dos etapas colapsa en un solo gran LP, el equivalente determinista:
Una copia de la segunda etapa por escenario, todas pegadas al mismo xx de la primera etapa. Entrégaselo a HiGHS o Gurobi, espera y obtén una solución. Fin de la historia.
Excepto que a veces la espera se convierte en toda la historia.
Cuando el equivalente determinista choca contra una pared
Para SS pequeños todo está bien. Escribes el LP, lo resuelves y te vas a casa. Pero los escenarios tienden a multiplicarse más rápido que tu paciencia:
Una distribución de demanda continua que se aproxima con una aproximación promedio de muestra podría necesitar miles de muestras antes de que se estabilice el xx óptimo. Un modelo de múltiples etapas se despliega exponencialmente con el número de etapas. Tres etapas con diez ramas cada una ya te dan mil hojas; cinco etapas y estás en los cientos de miles. En aplicaciones reales, como la programación de la energía hidroeléctrica, la planificación de la cadena de suministro ante crisis de demanda o el despacho de energía, no es inusual que haya decenas de miles de SS. La gente que trabaja en energía hidroeléctrica modela regularmente árboles de escenarios que son demasiado grandes para siquiera enumerarlos, y mucho menos entregarlos a un solucionador.
El equivalente determinista crece diligentemente con SS. Tiene un bloque de restricciones por escenario, un conjunto de variables ysy_s por escenario y un TsxT_sx que acopla cada bloque a la primera etapa compartida xx. Si nos fijamos en la matriz de restricciones, tiene una forma muy particular:
Esto se llama estructura de bloques angulares. Cada escenario aporta su propio bloque diagonal WsW_s, más una columna de TsT_s a la izquierda que lo vincula a xx. La fila compartida en la parte inferior es la restricción de viabilidad de la primera etapa Ax≥bAx geq b.
La forma es informativa, pero para quien la resuelve es una mala noticia. El tiempo de ejecución de LP crece de forma superlineal en el tamaño del problema. Una búsqueda rápida en Google le dirá que es aproximadamente O(n2.5)O(n^{2.5}) para métodos de estilo simplex y O(n3.5)O(n^{3.5}) para métodos de punto interior en el peor de los casos, lo que significa que duplicar el número de escenarios no duplica el tiempo de ejecución, sino que lo sextuplica. En algún momento, el problema deja de caber en la memoria por completo y su solucionador amigable comienza a cambiar al disco, lo cual es una forma educada de decir "esto es demasiado".
Así que necesitamos algo más inteligente que "meter todo en un solucionador y esperar lo mejor".
¿Por qué los trucos ingenuos no ayudan?
Antes de llegar a la idea inteligente, descartemos las obvias.
Idea 1: utilice menos escenarios. A veces esto está bien; a veces la respuesta cambia significativamente a medida que agregas escenarios, y tu “óptimo” fue una coincidencia de la muestra que extrajiste. La literatura muestral sobre la aproximación promedio trata precisamente sobre cuándo puedes salirte con la tuya y cuándo no.
Idea 2: Resuelva cada escenario por separado y promedie los resultados. Tentador, pero equivocado. El objetivo de tener un xx es que tiene que funcionar para todos los escenarios simultáneamente. Resolver el escenario ss por sí solo le brinda el óptimo específico del escenario (la solución de “esperar y ver” de la última publicación), no una decisión de primera etapa con la que pueda comprometerse.
Idea 3: Resuelva el problema determinista con demanda media 𝔼[ξ] y luego ocúpese de las sobras. Esta es la solución para vehículos eléctricos. Puede estar muy equivocado, y ese error es exactamente lo que mide VSS.
Así que no podemos reducir el problema, no podemos desvincularlo y no podemos solucionarlo en promedio. Necesitamos algo que respete la estructura del equivalente determinista sin resolverlo como un LP gigante.
Ese algo existe desde 1962.
La única observación que nos salva
Mire la matriz de restricciones nuevamente:
Ahora imagina que xx es fijo y elige cualquier valor, sin importar cuál. El bloque de la primera etapa Ax≥bAx geq b se cumple o no. Y el resto de las restricciones se desacoplan: la restricción del escenario s se convierte en Wsys≥hs−TsxW_s y_s ≥ h_s − T_s x, lo que involucra solo yₛ. La conexión entre escenarios pasa por xx. Arreglar xx, y los escenarios se desmoronan en LP independientes de SS.
Ésta es la idea clave. La variable xx es una variable complicada: si se conociera, la vida sería mucho más fácil. Cada subproblema es pequeño (el tamaño de un escenario, no todos) y son independientes, por lo que incluso puedes resolverlos en paralelo.
La receta de descomposición de Benders se basa en esta observación. Apenas:
Adivina una decisión de la primera etapa xx. Resuelva los subproblemas del escenario S de forma independiente, dado xx. Utilice lo que le indiquen los subproblemas para actualizar su estimación de xx. Repita hasta que sus conjeturas dejen de mejorar.
La parte difícil (y la parte elegante) es el paso 3: ¿cómo nos dicen los subproblemas algo útil sobre xx? Ahí es donde la dualidad LP se gana la vida.
¡Dobladores!
Reduzcamos la velocidad y construyamos el algoritmo pieza por pieza.
La función de valor v(ξˢ, x)
Para una decisión fija de primera etapa xx y un escenario fijo ss, el subproblema de la segunda etapa es solo un pequeño PL:
Su dual es
Este dual tiene una característica en la que nos apoyaremos. Su región factible Λs={λ≥0:λ⊤Ws≤qs}Lambda_s = {lambda geq 0 : lambda^top W_s leq q_s} es un poliedro y no depende de xx en absoluto. El xx solo aparece en el objetivo. Combinado con el hecho de que siempre se alcanza un óptimo en un vértice de la región factible, obtenemos una reescritura útil:
donde los λkslambda^s_k son los (finitos) vértices de ΛsLambda_s.
Esa es una expresión muy agradable para mirar. Como función de xx, v(ξs,x)v(xi^s, x) es el máximo puntual de una colección finita de funciones afines. Eso lo hace lineal y convexo por partes. Imagínese el sobre superior de un abanico de líneas.
El costo de recurso esperado es solo una suma ponderada de estos:
Una suma ponderada de funciones convexas lineales por partes es en sí misma lineal y convexa por partes. Entonces Q(x)Q(x), lo que estamos tratando de minimizar en la primera etapa además de c⊤xc^top x, es una función convexa lineal por partes de xx. No tenemos una fórmula explícita para ello, porque eso requeriría enumerar cada vértice de cada ΛsLambda_s, pero conocemos su forma.
Esa forma es la palanca que tira Benders.
Aproximar Q con cortes
No conocemos QQ en forma cerrada, pero una función convexa lineal por partes tiene una propiedad muy indulgente: en cualquier punto x‾bar{x} donde conocemos el valor Q(x‾)Q(bar x) y un subgradiente g∈∂Q(x‾)g ∈ ∂Q(x̄), la función lineal
es un límite inferior global de QQ. Toca QQ en x‾x̄ y permanece debajo en cualquier otro lugar.
Específicamente para el costo de recurso, el subgradiente es algo que podemos calcular. Si resolvemos los subproblemas duales y obtenemos vértices duales óptimos λk∗slambda^s_{k*} para cada escenario, un subgradiente de QQ en x‾x̄ es
Entonces, cada vez que resolvemos los subproblemas en algún candidato xτx^tau, obtenemos un límite inferior afín ajustado en QQ de forma gratuita. Benders llama a cada uno de estos límites inferiores un corte de optimización.
El plan ahora es visible. Vamos a construir una secuencia de límites inferiores afines en QQ, uno por iteración, y los usaremos como sustituto de QQ en el problema de la primera etapa.
El problema maestro
Reemplazamos Q(x)Q(x) con una variable marcador de posición θtheta que tiene que estar por encima de cada límite inferior que hemos acumulado hasta ahora. La iteración tt tiene el problema maestro:
Aquí LL es cualquier límite inferior seguro de QQ (a menudo L=0L = 0 si los costos de recurso no son negativos, o algo más negativo si no lo son) y cada (αt,βt)(αₜ, βₜ) es un corte de optimización de una iteración anterior. El problema maestro es simplemente un PL normal, mucho más pequeño que el equivalente determinista porque sólo incluye las variables de la primera etapa y un montón creciente de recortes.
La solución candidata (xτ,θτ)(x^tau, θ^tau) de este maestro luego se envía a los subproblemas, que devuelven un corte nuevo, que se agrega y volvemos a recorrerlo.
el algoritmo
En conjunto, el algoritmo de Benders uni-cut es:
1. Resuelva el problema maestro para obtener (xτ,θτ)(x^tau, θ^tau).
2. Para cada escenario s=1,…,Ss = 1, …, S, resuelva el subproblema dual en xτx^tau para obtener un vértice óptimo λτsλˢ_τ y el valor óptimo v(ξs,xτ)v(ξˢ, x^tau).
3. Construya el corte de optimización θ≥αt+βt⊤xθ ≥ α_t + β_t^top x con
y agréguelo al maestro.
4. Deténgase si el θτθ^tau del maestro ya satisface el nuevo corte en xτx^tau, es decir, si θτ≥αt+βt⊤xτθ^tau ≥ α_t + β_t^top x^tau. De lo contrario, regrese al paso 1.
Algunas cosas son discretamente destacables acerca de este bucle.
Finitud. Cada corte proviene de algún vértice λksλ^s_k de algunos ΛsLambda_s. Cada ΛsLambda_s tiene un número finito de vértices. Entonces, en el peor de los casos, el algoritmo enumera todos los cortes y paradas posibles: no puede realizar un bucle indefinido. En la práctica, se detiene mucho antes: normalmente sólo un pequeño puñado de vértices acaban importando.
Descomposición bien hecha. El gran LP compartido ha sido reemplazado por un pequeño problema maestro más S pequeños subproblemas independientes. Los subproblemas se pueden resolver en paralelo. El maestro crece, pero sólo con una restricción por iteración.
Información, no enumeración. Nunca escribimos todos los escenarios en un LP. Dejamos que los subproblemas nos digan, a través de sus variables duales, lo que el maestro necesita saber sobre ellos. Los cortes son el idioma que hablan.
Unicorte versus multicorte
En la receta anterior, agregamos la información por escenario en un solo corte en QQ. La alternativa natural es mantener los marcadores de posición separados de SS θ1,…,θSθ₁,…, θ_S, uno por escenario, y agregar cortes separados de SS en cada iteración. Esta es la variante multicorte.
El corte múltiple utiliza más información por iteración, por lo que la aproximación maestra de QQ se ajusta más rápido. El costo es que estás agregando restricciones SS por iteración en lugar de una, y el maestro se vuelve difícil de manejar más rápido. No hay un ganador universal: You y Grossmann (2013) informan aceleraciones significativas a partir de modelos de corte múltiple en la cadena de suministro, pero esto se analiza caso por caso. Un valor predeterminado razonable es comenzar con unicorte y probar con un corte múltiple si la convergencia es lenta.
Una imagen de una iteración.
A veces, la forma más fácil de asimilar el algoritmo es repasar geométricamente lo que logra cada iteración.
Imagine Q(x)Q(x) como la función convexa lineal por partes verdadera (desconocida) trazada sobre xx. La aproximación actual del maestro Q^τ(x)hat Q^tau(x) es el máximo puntual de todos los cortes que hemos agregado hasta ahora: un límite inferior poliédrico que se encuentra debajo de Q en todas partes.
En la iteración τtau, el maestro minimiza c⊤x+Q^τ(x)c^top x + Q̂^tau (x) y aterriza en algún xτx^tau donde su aproximación cree que es óptimo. Luego preguntamos a los subproblemas: “¿cuál es el verdadero QQ en xτx^tau?” Responden tanto con un valor Q(xτ)Q(x^tau) como con un subgradiente (la solución dual). Ese subgradiente nos da una nueva función lineal que toca a QQ en xτx^tau: el corte. Lo agregamos; la aproximación del maestro aumenta, lo suficiente para ser exacta en xτx^tau; y xτ+1x^{tau+1} de la siguiente iteración será (normalmente) diferente.
El ciclo se detiene cuando la aproximación del maestro ya predice el valor verdadero en xτx^tau. Visualmente, esto significa que nuestra envolvente inferior poliédrica ha alcanzado a QQ exactamente en el punto que nos interesa.
No siempre lograrás una coincidencia exacta porque, por razones numéricas, la mayoría de las implementaciones se detienen en una tolerancia pequeña, pero la geometría es la misma.
No es una herramienta mágica
Hasta ahora esto se ha leído como un almuerzo gratis algorítmico: la misma respuesta, LP mucho más pequeños, paralelizables, finitos. Seamos honestos también acerca de las desventajas.
Convergencia lenta en la naturaleza. En teoría, Benders termina en un número finito de iteraciones; en la práctica, “un número finito” puede significar bastantes iteraciones y cada una resuelve un PL. Rahmaniani et al. (2017) es la revisión bibliográfica estándar sobre la docena de trucos que la gente usa para acelerar las cosas: selección de cortes, eliminación de cortes cuando los cortes antiguos se vuelven dominados, regiones de confianza para evitar xτx^tau que oscilan violentamente, etc.
El problema principal se convierte en el cuello de botella. No es inusual que más del 90% del tiempo de ejecución total se dedique a resolver el maestro (en continuo crecimiento). Las soluciones comunes incluyen no resolver el máster de forma óptima en las primeras iteraciones, eliminar cortes que ya no ayudan y explotar la estructura cuando el máster es en sí mismo un LP estructurado.
Los subproblemas también pueden dominar. Especialmente cuando SS es grande o cada subproblema es en sí mismo un PL difícil, la segunda etapa toma el control. La solución es principalmente de ingeniería: paralelizar entre escenarios, iniciar en caliente cada subproblema a partir de la iteración anterior, agruparlos cuando sea posible.
El supuesto oculto: recurso relativamente completo. Todo lo anterior supone silenciosamente que el subproblema de la segunda etapa es factible para cualquier xx que el maestro pueda proponer. Esa es la propiedad del recurso (relativamente) completo de la publicación anterior: para cada primera etapa x factible y cada escenario ξξ, existe un recurso y≥0y ≥ 0 que satisface Wy≥h(ξ)−T(ξ)xWy ≥ h(ξ) − T(ξ)x.
Cuando esa propiedad falla, algo de xτx^tau dejará el subproblema inviable: el LP de la segunda etapa no tiene y que funcione, por lo que el costo implícito es +∞+ infty, lo que significa que, para empezar, xτx^tau nunca fue una opción factible en la primera etapa. Benders soluciona esto con cortes de viabilidad, que se derivan resolviendo un LP auxiliar que mide qué tan inviable es el subproblema y usando su dual para producir un hiperplano que corta xτx^tau de la región factible del maestro. Sin pasar por el álgebra, la estructura es la misma que la de los cortes de optimización: resuelve un pequeño PL, mira el dual, construye una restricción lineal y devuélveselo al maestro. Entonces, en cada iteración de Benders, se agrega un corte de optimización (si cada subproblema era factible) o un corte de factibilidad (si algún subproblema no lo era). El modelo de recurso sin garantía de recurso total todavía funciona; solo necesitas ambos tipos de cortes.
Sólo se aplica plenamente a los modelos de dos etapas. Para programas estocásticos de múltiples etapas, la generalización natural no es simplemente "Benders con más etapas". Es programación dinámica dual estocástica (SDDP), que utiliza funciones de muestreo y valor aproximado para evitar construir el árbol de escenarios completo. El SDDP es el caballo de batalla en la programación de la energía hidroeléctrica y merece, y probablemente con el tiempo obtendrá, su propio puesto.
Las primeras etapas de números enteros se vuelven complicadas. Si xx tiene que ser un número entero (piense en decisiones de capacidad, ubicación de las instalaciones), entonces el maestro es un MILP, y la teoría convexa limpia de Benders pierde parte de su fuerza. Hay extensiones, como Benders generalizados, números enteros en forma de L, corte de rama y Benders, etc., pero están fuera del alcance de esta publicación.
Nada de esto significa que Benders sea malo. Significa que es un método, no un milagro. Usado cuidadosamente, con la división maestra/subproblema correcta, la variante de corte correcta y los trucos de aceleración correctos, es la diferencia entre resolver un problema o no.
Concluyendo
Comenzamos con un problema: los programas estocásticos son fáciles de escribir pero fáciles de hacer demasiado grandes para resolverlos. El equivalente determinista crece linealmente en el número de escenarios, y los solucionadores de LP crecen superlinealmente en el tamaño del problema, lo que es un curso de colisión cuadrático o peor.
La descomposición de Benders toma en serio esa estructura angular de bloques. Elimina las complicadas variables de la primera etapa, deja la segunda etapa desacoplada del escenario a un enjambre de pequeños subproblemas y utiliza la dualidad LP para transportar información entre las dos: la optimización se corta cuando los subproblemas son factibles, la factibilidad se corta cuando no lo son. El resultado es un esquema iterativo que converge en un número finito de pasos y muy a menudo mucho más rápido.
Para resumir lo que realmente está haciendo el trabajo:
La función de valor de segunda etapa v(ξ,x)v(ξ, x) es lineal por partes y convexa en xx. También lo es el coste de recurso esperado Q(x)Q(x). Una función convexa lineal por partes se puede aproximar arbitrariamente bien mediante una colección cada vez mayor de límites inferiores afines (cortes), cada uno de los cuales se obtiene de forma gratuita a partir de una solución dual. La primera etapa no necesita conocer QQ exactamente, solo limitarlo estrechamente en los puntos que le interesan. El bucle maestro/subproblema crea exactamente eso. Unicorte y multicorte son dos formas de agrupar la información por escenario. Los recortes de viabilidad abordan el caso en el que los subproblemas pueden fallar.
Si quitas una cosa, deja que sea el aspecto de bloque angular de la matriz de restricciones. Siempre que puedas reescribir un problema de optimización para que arreglar algunas variables haga que el resto sea separable, puedes probar Benders. Esta idea aparece en el diseño de redes, la programación y la planificación energética, en cualquier lugar donde se tenga un conjunto de decisiones “difíciles” que unen otras que de otro modo serían fáciles.
En la próxima publicación de esta serie, veremos qué hacer cuando incluso un Benders de dos etapas es el marco equivocado, cuando el mundo no se divide cortésmente en "decidir, observar, corregir" sino "decidir, observar, decidir, observar, decidir…" indefinidamente. Ese es el territorio del SDDP. Si alguna vez se ha preguntado cómo los operadores hidroeléctricos programan los embalses en un horizonte de varios años con miles de escenarios de afluencia conjunta, la respuesta está en la próxima publicación.
Créditos y referencias
Al igual que la publicación anterior, esta se basa en las conferencias del Dr. Ruben van Beesten (Universidad Noruega de Ciencia y Tecnología) de su curso de octubre de 2023 sobre Programación Estocástica, al que tuve el placer de asistir en Trondheim. La preparación pedagógica proviene directamente de sus diapositivas; cualquier torpeza en el recuento es mía. Además, todas las imágenes, excepto la miniatura, de este post están hechas por mí.
Cuatro referencias que vale la pena conocer:
Benders, JF (1962). Procedimientos de particionamiento para la resolución de problemas de programación de variables mixtas. Numerische Mathematik, 4(1), 238–252. El documento original. Sorprendentemente legible. Van Slyke, R. y Wets, R. (1969). Programas lineales en forma de L con aplicaciones para control óptimo y programación estocástica. Revista SIAM de Matemáticas Aplicadas, 17(4), 638–663. La versión con sabor a programación estocástica; De aquí proviene el nombre de “método en forma de L”. Fengqi You e Ignacio E Grossmann. (2013). Algoritmo de descomposición de dobladores multicorte para la planificación de la cadena de suministro de procesos en condiciones de incertidumbre. Anales de investigación de operaciones, 210:191–211. Rahmaniani, R., Crainic, TG, Gendreau, M. y Rei, W. (2017). El algoritmo de descomposición de Benders: una revisión de la literatura. Revista europea de investigación operativa, 259(3), 801–817. La referencia para aceleraciones, variantes y uso moderno.
Y las publicaciones anteriores de la serie: Cinco preguntas que le ayudarán a modelar mejor programas lineales enteros, Cuatro pasos para fortalecer su programa lineal y Una suave introducción a la programación estocástica.