Optimizacion Estocastica

OPTIMIZACIÓN ESTOCÁSTICA
Andrés Ramos
Santiago Cerisola
Septiembre 2010
[http://www.doi.icai.upcomillas.es/simio/intro.htm]
Alberto Aguilera 23 – E 28015 Madrid – Tel: 34 91 542 2800 – Fax: 34 91 541 1132 – www.doi.icai.upcomillas.es
ÍNDICE
I.1 INTRODUCCIÓN ........................................................................................... 3

I.1.1 Reseña histórica ....................................................................... 12
I.2 OPTIMIZACIÓN LINEAL BIETAPA Y MULTIETAPA ......................................... 15
I.3 TÉCNICAS DE DESCOMPOSICIÓN ................................................................. 17
I.4 DESCOMPOSICIÓN DE BENDERS ................................................................. 25
I.4.1 Descripción general ................................................................... 25
I.4.2 Cortes de infactibilidad .............................................................. 28
I.4.3 Algoritmo................................................................................. 31
I.4.4 Caso ejemplo ............................................................................ 34
I.4.5 Algoritmo escrito en GAMS ....................................................... 35
I.4.6 Algoritmo escrito en MATLAB ..................................................... 43
I.4.7 Problema de transporte con coste fijo .......................................... 44
I.4.8 Variables duales del problema lineal bietapa ................................. 49
I.5 DESCOMPOSICIÓN DE DANTZIG-WOLFE ..................................................... 51
I.5.2 Algoritmo................................................................................. 56
I.5.4 Algoritmo escrito en MATLAB ..................................................... 69
I.6 RELAJACIÓN LAGRANGIANA ...................................................................... 71
I.6.2 Cortes de acotamiento ............................................................... 72
I.6.3 Algoritmo................................................................................. 73
I.6.5 Problema de transporte con coste fijo .......................................... 91
I.7 OTRAS TÉCNICAS DE DESCOMPOSICIÓN ...................................................... 95
I.7.1 Descomposición Primal-Dual...................................................... 95
I.7.2 Descomposición anidada ............................................................ 96
I.7.3 Caso ejemplo .......................................................................... 100
I.7.4 Modelo hidrotérmico ............................................................... 102
I.8 OPTIMIZACIÓN LINEAL ESTOCÁSTICA BIETAPA Y MULTIETAPA ................. 109
I.8.1 Bietapa .................................................................................. 109
I.8.2 Multietapa.............................................................................. 112
I.9 MEJORAS EN LAS TÉCNICAS DE DESCOMPOSICIÓN .................................... 115
I.10 SIMULACIÓN EN OPTIMIZACIÓN LINEAL ESTOCÁSTICA BIETAPA ............... 119
I.10.1 Técnicas de reducción de varianza .......................................... 121
I.10.1.1 Números aleatorios comunes ........................................................................................ 123
I.10.1.2 Variables antitéticas ..................................................................................................... 124
I.10.1.3 Variable de control ....................................................................................................... 125
04/11/2012 i
I.10.1.4 Muestreo por importancia............................................................................................. 127
I.10.1.5 Muestreo estratificado .................................................................................................. 133
I.11 REFERENCIAS ........................................................................................ 137
I.12 BIBLIOTECA DE PROBLEMAS ................................................................... 143
ii 04/11/2012
I OPTIMIZACIÓN ESTOCÁSTICA
I.1 Introducción
La estocasticidad o incertidumbre aparece en todos los sistemas pero hasta
ahora no era posible la solución de problemas de optimización de grandes
sistemas considerando explícitamente ésta. La incertidumbre puede deberse a
carencia de datos fiables, errores de medida o tratarse de parámetros que
representan información sobre el futuro. Por ejemplo, en el caso de planificación
de sistemas de energía eléctrica la incertidumbre surge principalmente en: la
demanda y precios futuros de la electricidad o de los combustibles, las
aportaciones hidráulicas o la disponibilidad de los elementos de generación y
red. No toda la incertidumbre se encuentra en el mismo horizonte temporal.
En optimización determinista se supone que los parámetros del problema son
conocidos con certeza, aunque sea a su valor medio. En optimización estocástica
(stochastic optimization SP) se relaja esta condición. No se conocen sus valores,
sólo sus distribuciones y habitualmente se supone que éstas son discretas con un
número finito de estados posibles. La suposición de distribuciones discretas es
habitual en los optimizadores de optimización estocástica. Actualmente no
existen aplicaciones estándar o comerciales, potentes y fiables, para resolver
problemas estocásticos. Todavía no es un campo en desarrollo.
Los tipos de modelos que aparecen en programación lineal1 estocástica son
motivados principalmente por problemas con decisiones de tipo aquí y ahora
(here and now) [Dantzig:55], decisiones previas bajo futuro incierto2. Esto es,
decisiones que deben tomarse basándose en información a priori, existente o
supuesta, sobre situaciones futuras sin realizar observaciones adicionales.
Recurso3 es la capacidad de tomar una acción correctora después de que haya
1
La restricción a problemas lineales es por razones prácticas: primero, son los más
frecuentemente utilizados y segundo, aquéllos con métodos de solución más eficientes y
avanzados. Las técnicas que se presentan admiten también funciones no lineales mientras la
región factible y la función objetivo sean convexas.
2
Existen otros tipos de problemas de optimización estocástica, que no se tratan en el
documento, denominados con restricciones probabilistas (chance constraints o probabilistic
constraints) que incluyen restricciones lineales con parámetros aleatorios y se modela la
condición de cumplirlas restricciones con una cierta probabilidad. En el caso de distribuciones
discretas se necesita una variable binaria por cada suceso, definido como combinación de los
valores discretos de los parámetros.
3
El recurso se denomina completo cuando para cualquier decisión de la primera etapa
(independientemente de su factibilidad) existen decisiones factibles para la segunda etapa. Se
denomina relativamente completo si se cumple sólo para decisiones factibles de la primera etapa
y parcial cuando no siempre las decisiones de la segunda etapa son factibles para las decisiones
de la primera etapa.
04/11/2012 3
1 INTRODUCCIÓN
ocurrido un suceso aleatorio. Por ejemplo, se toman hoy un conjunto de

decisiones x 1 con valores de los parámetros conocidos (es decir, deterministas),
durante la noche se producen unos sucesos aleatorios (exógenos) y mañana se
toman un conjunto de acciones correctoras x 2ω que mitigan (corrigen) los efectos
de los sucesos aleatorios ω sobre las decisiones de hoy.
Las decisiones de la primera etapa no dependen del escenario que ocurra
realmente en el futuro. Las decisiones de la segunda etapa x 2ω son los recursos.
La función objetivo asociada a estas decisiones es la función de recursos,
depende de las decisiones de la primera etapa x 1 y de los valores ω de los
parámetros aleatorios. Esta circunstancia incrementa exponencialmente el
tamaño del problema. Este problema, denominado problema lineal estocástico
bietapa4 PLE-2, determina el valor óptimo de x 1 y x 2ω . Las decisiones relevantes
en optimización estocástica son exclusivamente las de la primera etapa x 1 ,
puesto que son éstas las que hay que tomar de manera inmediata.
Veamos un ejemplo sencillo de planificación lineal estocástica bietapa. Se
trata de un problema de planificación óptima de la expansión de la generación,
es decir, determinar la inversión óptima en generación. Se dispone de varios
tipos de generadores candidatos para satisfacer la demanda durante varios
periodos. El primero de los periodos tiene una demanda aleatoria. La decisión de
inversión debe ser única para todos los periodos. Además existe una restricción
presupuestaria a la inversión total y una potencia mínima a instalar.
Veamos la formulación matemática del problema. Sea fi el coste fijo de
inversión de cada generador candidato i en el alcance del modelo, v ji el coste
variable de producción del generador i en el periodo j , dsj la demanda de
potencia de cada periodo j para cada escenario s siendo ps la probabilidad del
escenario. Se define una potencia mínima a invertir P y un presupuesto
máximo disponible D . Las variables del problema serán X i la potencia a
instalar de cada generador i , Yjis potencia de operación de cada generador i en
cada periodo j para cada escenario s .
En primer lugar, veamos un escenario determinista cualquiera. La función
objetivo será minimizar la suma de costes fijos de inversión más los costes de
operación para dicho escenario.
4
Se denomina bietapa porque tiene dos etapas, una asociada a las decisiones anteriores a la
realización de la incertidumbre y otra a las decisiones correctoras una vez realizada la
incertidumbre.
4 04/11/2012
min ∑ fi X i + ∑ v jiYjis
XiYsji
i ji
∑X i ≥P
i
∑fX
i
i i ≤D
para cualquier escenario s
s
Y ≤ Xi
ji ∀ji
s
∑v Y ji ji ≥ d js ∀j
i
X i ,Yjis ≥ 0
Veamos ahora el problema estocástico de minimizar los costes totales
esperados para el conjunto de los tres escenarios
min ∑ fi Xi + ∑ ps v jiYjis
XiYsji
i sji
∑ Xi ≥ P
i
∑fX i i
≤D
i
Yjis ≤ X i ∀sji
s s
∑Y ji
≥d j
∀sj
i
X i ,Yjis ≥ 0
$title Planificación óptima de la expansión de la generación

sets
I generadores / gen-1 * gen-4 /
J periodos / per-1 * per-3 /
S escenarios de demanda / s-1 * s-3 /
parameters
F(i) coste fijo de inversión [€]
/ gen-1 10
gen-2 7
gen-3 16
gen-4 6 /
PROB(s) probabilidad de cada escenario [p.u.]
/ s-1 0.2
s-2 0.5
s-3 0.3 /
DEM(j) demanda para un escenario [MW]
table V(i,j) coste variable de operación [€ por MW]
per-1 per-2 per-3
gen-1 40 24 4
gen-2 45 27 4.5
gen-3 32 19.2 3.2
gen-4 55 33 5.5
table DEMS(s,j) demanda estocástica [MW]
per-1 per-2 per-3
s-1 3 3 2
s-2 5 3 2
s-3 7 3 2
04/11/2012 5
1 INTRODUCCIÓN
scalar
POTMIN potencia mínima a instalar [MW] / 12 /
PRSPTO límite presupuestario [€] / 120 /
variables
X(i) potencia a instalar [MW]
Y(j,i) potencia en operación [MW]
YS(s,j,i) potencia en operación estocástica [MW]
COSTE coste total
positive variables X, Y, YS
equations
COST coste total [€]
COSTS coste total estocástico [€]
PRESUP limitación presupuestaria [€]
INSMIN potencia mínima instalada [MW]
BALPOT potencia en operación menor que instalada [MW]
BALPOTS potencia en operación menor que instalada estocástica [MW]
BALDEM balance de demanda [MW]
BALDEMS balance de demanda estocástico [MW] ;
COST .. COSTE =E= sum(i, F(i) * X(i))

sum
+ sum((j,i),
sum V(i,j) * Y(j,i)) ;
COSTS .. COSTE =E= sum(i, F(i) * X(i))
sum
+ sum((s,j,i), PROB(s) * V(i,j) * YS(s,j,i)) ;
sum
PRESUP .. sum
sum(i, F(i) * X(i)) =L= PRSPTO ;
INSMIN .. sum
sum(i, X(i)) =G= POTMIN ;
BALPOT(j,i) .. Y(j,i) =L= X(i) ;
BALPOTS(s,j,i) .. YS(s,j,i) =L= X(i) ;
BALDEM(j) .. sum
sum(i, Y(j,i)) =G= DEM(j) ;
BALDEMS(s,j) .. sum
sum(i, YS(s,j,i)) =G= DEMS(s,j) ;
model DETERM / COST, INSMIN, PRESUP, BALPOT, BALDEM / ;
model ESTOCA / COSTS, INSMIN, PRESUP, BALPOTS, BALDEMS / ;
* cada escenario determinista por separado
loop (s,
DEM(j) = DEMS(s,j) ;
solve DETERM MINIMIZING COSTE USING LP ;
) ;
* escenario de demanda media
DEM(j) = sum
sum(s, PROB(s) * DEMS(s,j)) ;
* problema estocástico
solve ESTOCA MINIMIZING COSTE USING LP
El problema se resuelve en primer lugar de forma determinista para cada

escenario s de demanda por separado. A continuación se resuelve el problema
para el valor medio de la demanda y finalmente el problema estocástico y sus
resultados se presentan a continuación.
Determinista 1 Determinista 2 Determinista 3 Medio Estocástico

Generador 1 . 0.33 3.67 0.67 0.67
Generador 2 . . . . 2
Generador 3 3 4.67 3.33 4.53 4.33
Generador 4 9 7 5 6.8 5
Coste total 262 346.67 437.33 355.73 362.47
6 04/11/2012
Para la solución de problemas estocásticos se ha recurrido históricamente a

diferentes aproximaciones:
i) Una es solucionar el problema para cada posible escenario. Ésta se

denomina espera y observa (wait and see o scenario analysis o what-if
analysis), corresponde a tomar las decisiones una vez que se ha resuelto la
incertidumbre. Estas decisiones se denominan clarividentes, anticipativas, no
implantables. El problema se resuelve independientemente para cada
escenario. Un caso particular de problema determinista es el escenario medio,
correspondiente al valor medio de los parámetros y cuya función objetivo
denominamos (expected value EV). Potencialmente, se obtiene una solución
diferente para cada escenario (como sucede en el ejemplo anterior) e incluso
la solución de un escenario puede ser infactible en otro.
ii) A partir de las soluciones se observan sistemáticamente semejanzas y
diferencias entre ellas y mediante criterios heurísticos subjetivos elegir planes
robustos, es decir, planes cuyas decisiones forman parte de los planes óptimos
bajo numerosos escenarios o planes flexibles, cuando las decisiones que
componen un plan pueden ser modificadas conforme se resuelve la
incertidumbre. Evidentemente, este procedimiento no garantiza la obtención
de la solución óptima. De hecho, la solución óptima estocástica no tiene por
qué ser óptima para ningún escenario. Con un criterio heurístico el generador
2 en el problema de expansión previo no aparecería puesto que no aparece en
la solución para ninguno de los escenarios ni para el de demanda media. Sin
embargo, sí aparece en la solución óptima del problema estocástico.
El paso de estas aproximaciones a la optimización estocástica supone un

salto cualitativo en cuanto a complejidad de solución y tamaño del problema.
La optimización estocástica se debe utilizar cuando las decisiones pueden no ser
utilizables o realizables si no se tiene en cuenta expresamente la incertidumbre,
es una forma de protegerse frente a ella. En el ejemplo anterior se aprecia que
las decisiones de inversión en las que se está interesado cambian según se haga
la hipótesis sobre la ocurrencia de cada escenario futuro o del escenario de
demanda media.
Las decisiones en un momento dado dependen exclusivamente de la
información disponible en ese momento, no pueden utilizar información futura.
Las decisiones de la primera etapa (decisión única de inversión) son
independientes del escenario que ocurra, se toman con antelación a resolver la
incertidumbre. Esta característica se denomina propiedad de implantabilidad o
no anticipatividad de las decisiones.
04/11/2012 7
1 INTRODUCCIÓN
La optimización estocástica permite incorporar actitudes aversas al riesgo,

por ejemplo, mediante la penalización de los peores escenarios. La generalización
del concepto de protección frente al riesgo se hace mediante optimización
robusta. Se dice que una solución es robusta si es casi óptima para todos los
escenarios. Se dice que un modelo es robusto si es casi factible en todos los
escenarios. La optimización robusta trata de encontrar un equilibrio simultáneo
entre ambos objetivos de optimalidad y factibilidad.
Los diferentes estados que pueden tomar los parámetros aleatorios a lo largo
del tiempo se representan mediante un árbol de probabilidad o de escenarios,
figura 1.1. Se define un escenario como cualquier camino (trayectoria) que va
desde la raíz del árbol hasta las hojas de manera que los escenarios que
comparten una misma información hasta una cierta etapa también comparten
esa parte del árbol. El árbol de probabilidad es la forma natural y explícita de
representar la no anticipatividad de las decisiones estocásticas discretas. En el
caso ejemplo previo, el árbol es simplemente una raíz y tres hojas. Las
decisiones de la primera etapa se comparten (son las mismas para cualquiera de
los escenarios) y por eso la raíz es única. Las decisiones de la segunda etapa son
múltiples (dependen de cada escenario) y por eso tiene tres hojas. La
determinación del árbol de probabilidad debe considerar las dependencias
temporales y/o espaciales que pudieran existir entre los parámetros aleatorios.
Si el proceso aleatorio es markoviano la dependencia temporal sólo alcanzará un
periodo. En un caso general se puede extender más allá.
En la formulación que se utiliza en este ejemplo y en algunas de las técnicas
de descomposición (e.g., descomposición de Benders) que se van a explicar en el
documento se utiliza expresamente la estructura del árbol de probabilidad para
definir las variables y los parámetros aleatorios. Así en el ejemplo, las variables
de la primera etapa no dependen de los escenarios mientras que los de la
segunda sí. Existe otra técnica de resolución alternativa (scenario decomposition
with splitting variable formulation) que considera variables independientes para
cada escenario (como se muestra en la figura 1.2) y establece restricciones de no
anticipatividad que hacen que ciertas variables sean iguales entre sí en ciertos
periodos. En esta técnica se utiliza el método de relajación lagrangiana que se
verá más adelante, especialmente indicado cuando existe un número reducido de
restricciones de complicación. Es decir, el árbol de la figura 1.1 se convierte en el
de la figura 1.2 mediante restricciones. En el ejemplo previo sería equivalente a
haber definido una variable de inversión para cada escenario y haber impuesto
restricciones que las obligan a ser iguales.
8 04/11/2012
Y1
X Y2
Y3
Figura 1.1 Árbol de probabilidad.
X1 Y1
Y2
X2
X3 Y3
Figura 1.2 Árbol de probabilidad para agregación de escenarios.
La estructura de la matriz de restricciones del problema estocástico se

observa en la figura 1.3. En ella se muestran las restricciones que involucran
variables de cada etapa por separado y las que unen variables de ambas etapas.
El problema estocástico es de tamaño mucho mayor que el de un problema
determinista y tiene la estructura no anticipativa del árbol de probabilidad. Las
técnicas de descomposición aprovechan esta estructura repetitiva de la matriz de
restricciones.
Variables Variables
primera segunda
etapa etapa
X Y1, Y2, Y3
PRESUP Restricciones primera etapa
INSMIN
BALPOTS
BALDEMS Restricciones primer escenario
BALPOTS
BALDEMS Restricciones segundo escenario
BALPOTS
BALDEMS Restricciones tercer escenario
Figura 1.3 Estructura de la matriz de coeficientes de las restricciones.
04/11/2012 9
1 INTRODUCCIÓN
La función objetivo del problema estocástico (recourse problema RP) es el

valor de la parte correspondiente a las decisiones de la primera etapa más el
valor esperado de la correspondiente a las decisiones de la segunda etapa (362.47
en el ejemplo). Se denomina valor esperado con información perfecta (EVWPI)
a la suma ponderada para cada escenario de la función objetivo total sabiendo
que dicho escenario va a ocurrir con certeza (356.93 para el ejemplo). Es decir,
se revela la incertidumbre antes de tomar las decisiones de la primera etapa.
Este valor siempre será menor o igual, si se está minimizando, que la función
objetivo del problema estocástico. El valor de la solución estocástica (value of
stochastic solution VSS) es la diferencia entre la función objetivo del valor
esperado de la solución del valor medio de los parámetros estocásticos EEV
(280x0.2+347.73x0.5+454.73x0.3=366.27) y del problema estocástico (366.28–
362.47=3.81).
Para cada escenario, la solución del problema estocástico es siempre peor o
igual que la solución con información perfecta (280, 349.33 y 439.33
respectivamente). Se denomina valor esperado de la información perfecta
(expected value of perfect information EVPI) o arrepentimiento esperado a la
suma ponderada de la diferencia entre la solución del problema estocástico en
cada escenario y la solución con información perfecta en dicho escenario. El
arrepentimiento en cada escenario es 280–262=18, 349.33–346.67=2.66, 439.33–
437.33=2 y el arrepentimiento esperado EVPI=18x0.2+2.66x0.5+2x0.3=5.54.
Este valor siempre será positivo en un problema de minimización.
Cuando esta diferencia en algún escenario pueda ser elevada (de
consecuencias catastróficas) o existe gran no linealidad en la función objetivo
para los diferentes escenarios la optimización estocástica se puede formular como
minimización del máximo arrepentimiento (o criterio minimax). Éste es un
criterio conservador de protección frente el riesgo que se utiliza cuando no se
tienen o no se quieren utilizar las probabilidades de los escenarios porque no se
tiene confianza en dichos valores. Otro criterio semejante y complementario
sería el de maximización del mínimo beneficio.
Se cumplen las relaciones EVWPI + EVPI = RP , RP + VSS = EEV y, por
consiguiente, EVWPI ≤ RP ≤ EEV .
Veamos otro ejemplo de un problema de optimización de una compañía de
gas con demanda y precio estocásticos tomado de NEOS (http://www-
fp.mcs.anl.gov/otc/Guide/OptWeb/continuous/constrained/stochastic/index.ht
ml). Cuando la compañía compra gas suministra parte a sus consumidores y el
resto lo almacena. Cuando vende gas lo toma de sus compras o de sus depósitos.
Luego sus decisiones son cuánto gas comprar y vender directamente, cuánto
comprar y almacenar o cuánto tomar del depósito y vender. Las decisiones
óptimas dependerán del precio del gas ahora y en el futuro, del coste de
10 04/11/2012
almacenamiento, del tamaño del depósito y de la demanda en cada periodo. La

demanda anual y el precio del gas para cada escenario se presentan en la
siguiente tabla. El precio de almacenamiento de una unidad de gas es de 1
€/año. Se quieren tomar las decisiones óptimas para los dos primeros años
sabiendo que el año actual está siendo normal y el año próximo puede ser
normal, frío o muy frío con las probabilidades que aparecen en la tabla.
Escenario Coste del gas (€) Demanda Probabilidad

Normal 5.0 100 1/3
Frío 6.0 150 1/3
Muy frío 7.5 180 1/3
$title Planificación óptima de la compra de gas

sets
S escenarios de temperatura / nr, fr, mf /
O opciones de gestión / cyv, cya, ayv /
parameters
CVS(s) coste del gas en cada escenario [€]
/ nr 5.0
fr 6.0
mf 7.5 /
PROBS(s) probabilidad de cada escenario [p.u.]
/ nr 0.333333
fr 0.333333
mf 0.333334 /
DEMS(s) demanda en un escenario
/ nr 100
fr 150
mf 180 /
CV coste del gas
PROB probabilidad
DEM demanda en un escenario
scalar
CALM coste de almacenamiento [€ por año] / 1 /
variables
X(o) cantidad que gestiona el primer año
Y(o) cantidad que gestiona el segundo año
YS(s,o) cantidad que gestiona el segundo año estocástico
COSTE coste total
positive variables X, Y, YS
equations
COST coste total
COSTS coste total estocástico
BALDEM1 balance de demanda del primer año
BALDEM2 balance de demanda del segundo año
BALDEMS2 balance de demanda del segundo año
GSTDEP gestión del depósito ;
COST .. COSTE =E= CVS('nr') * (X('cyv')+X('cya'))
+ CALM * X('cya')
+ CV * (Y('cyv')+Y('cya')) ;
COSTS .. COSTE =E= CVS('nr') * (X('cyv')+X('cya'))

+ CALM * X('cya')
+ sum
sum(s, PROBS(s) * CVS(s) * (YS(s,'cyv')+YS(s,'cya'))) ;
BALDEM1 .. X('cyv') =E= DEMS('nr') ;
BALDEM2 .. X('ayv') + Y('cyv') =E= DEM ;

BALDEMS2(s) .. X('ayv') + YS(s,'cyv') =E= DEMS(s) ;
04/11/2012 11
1 INTRODUCCIÓN
GSTDEP .. X('cya') =E= X('ayv') ;

model DETERM / COST , BALDEM1, BALDEM2 , GSTDEP / ;
model PROBAB / COSTS, BALDEM1, BALDEMS2, GSTDEP / ;
* cada escenario determinista por separado
loop (s,
CV = CVS(s) ;
DEM = DEMS(s) ;
) ;
* escenario de demanda media

CV = sum(s,
sum PROBS(s) * CVS(s)) ;
DEM = sum(s,
sum PROBS(s) * DEMS(s)) ;
* problema estocástico
solve PROBAB MINIMIZING COSTE USING LP ;
Los resultados para cada escenario determinista, el escenario medio y el

problema estocástico se presentan a continuación.
Compra y Compra y Depósito a Compra y Compra y Coste

Vende venta vende
almacena almacena
Año 1 Año 2
Normal 100 . . 100 . 1000
Frío 100 . . 150 . 1400
Muy frío 100 180 180 . . 1580
Medio 100 143.33 143.33 . . 1360
. .
Estocástico 100 100 100 50 . 1400
80 .
Para este ejemplo también se observa que las soluciones en cada escenario
son claramente diferentes a las del escenario medio y a las decisiones del
problema estocástico.
I.1.1 Reseña histórica

La optimización estocástica apareció en 1955 como una extensión de la
programación lineal con énfasis en el gran número de variables y parámetros con
trabajos independientes de Dantzig [Dantzig:55] y Beale [Beale]. Por otra parte,
como otra extensión de la programación lineal para grandes sistemas con
estructuras especiales en la matriz de coeficientes de las restricciones aparecieron
las técnicas de descomposición [Benders, Dantzig:63, Dantzig:60, VanSlyke],
también denominadas de optimización matemática a gran escala. Los problemas
de optimización de gran tamaño con estructura especial pueden tener
subyacente una estructura espacial, temporal o debida a la aleatoriedad de los
12 04/11/2012
parámetros. Esta estructura es la que se rompe mediante el uso de las técnicas

de descompoción para la resolución del problema.
Aunque las primeras investigaciones aparecen muy temprano, sólo
recientemente el avance en la tecnología de los ordenadores ha permitido la
solución de problemas de muy gran tamaño y ha devuelto interés al tema de la
optimización estocástica produciendo además un avance en la teoría matemática
que lo sustenta. Los libros de Ermoliev y Wets (eds.) [Ermoliev], Kall y Walace
[Kall] (http://www.unizh.ch/ior/Pages/Deutsch/Mitglieder/Kall/bib/ka-wal-
94.pdf), Birge y Louveaux [Birge:97] y Prékopa [Prekopa] pueden servir de
compendio a la investigación reciente en este campo. En castellano como
refencia se pueden utilizar los libros de Ramos et al. [Ramos:08] y Alonso-Ayuso
et al. [Alonso:04]. Otros libros se pueden encontrar en
(http://stoprog.org/index.html?booksSP.html). En internet existe una dirección
que recoge las actividades de la Stochastic Programming Community
(http://stoprog.org/) y un repositorio de publicaciones electrónicas recientes
sobre programación estocástica Stochastic Programming E-Print Series
(http://www.speps.info/) así como una recopilación de bibliografía en
optimización estocástica (http://mally.eco.rug.nl/biblio/SPlist.html) y
optimización estocástica entera (http://mally.eco.rug.nl/biblio/SIP.HTML).
Existe específicamente una conferencia trianual dedicada a la optimización
estocástica The Internacional Conference on Stochastic Programming
(http://tucson.sie.arizona.edu/SPX) que ya va por su décima edición.
La resolución de este tipo de problemas abre nuevas posibilidades en
numerosos ámbitos. En el ámbito de la planificación, conceptualmente
asignación óptima de recursos, de sistemas de energía eléctrica se pueden citar
como ejemplos modelos de:
• Planificación a largo plazo de la expansión de la red [Perez, Stanford]

La expansión de la red se formula como un problema estocástico bietapa.
La primera etapa recoge las decisiones discretas de inversión en
corredores y la segunda las decisiones de explotación. Los distintos
escenarios se deben a variaciones en la demanda, la hidraulicidad y la
disponibilidad de elementos de generación y red. En el modelo PERLA la
expansión de la red se formula como un problema de planificación
estática que optimiza el coste total anual, suma del coste anualizado de
inversión, explotación e indisponibilidad. Este ha sido utilizado por REE
como herramienta de ayuda a la planificación táctica (con horizontes de
15 a 30 años). El modelo ha sido desarrollado en FORTRAN y utiliza
MINOS [Murtagh] y ZOOM para la solución de los problemas de
optimización lineal y entera respectivamente.
04/11/2012 13
1 INTRODUCCIÓN
• Coordinación hidrotérmica [Morton:93, Jacobs, Morton:96, Pereira]

Ésta consiste en la planificación de la explotación de un sistema
hidrotérmico en múltiples periodos considerando las aportaciones
hidráulicas como parámetro aleatorio. El objetivo es determinar la
explotación del subsistema hidráulico que minimiza los costes de
explotación esperados derivados de la utilización del subsistema térmico.
• Determinación del precio de la electricidad para diferentes tipos de
consumidores en sistemas hidráulicos [Mathiesen]
Este modelo analiza la formación de precios en el nuevo mercado de
electricidad de Noruega. Explícitamente considera la gran variación
interanual entre aportaciones hidráulicas (estocásticas) y el reparto anual
del agua teniendo en cuenta los diferentes estados. Considera varias
regiones y periodos en el año. El modelo propuesto se formula como
problema no lineal (cuadrático) estocástico multietapa. La principal
hipótesis es la estrecha relación entre el precio de la electricidad y otros
precios de la energía. Se supone que todos los agentes son neutrales frente
al riesgo y compran al precio estipulado. Se distingue entre dos tipos de
consumidores: aquellos que no tienen fuentes de energía alternativas
(demanda contractual cuyo precio está estipulado y garantizado durante
el año) y aquellos que sí (demanda puntual cuyo precio es incierto hasta
el comienzo del periodo).
• Planificación a largo plazo de la expansión de la generación [MIT,
Stanford] o de la red de distribución
• Planificación a largo plazo de la expansión de inversiones en elementos de
potencia reactiva
• Programación semanal de la explotación de un sistema eléctrico
• Optimización a corto plazo del margen de contribución de una compañía
eléctrica con incertidumbre en el comportamiento de los competidores
Otros ámbitos de aplicación pueden ser:
• Planificación de capacidad de una línea de producción [Ramos]

• Gestión de una cartera de valores [Dantzig:91c], [Mulvey]
• Planificación de la instalación de las subestaciones de alimentación al
AVE y los tipos de catenaria
14 04/11/2012
I.2 Optimización lineal bietapa y multietapa

El problema lineal bietapa PL-2 se representa matemáticamente de la forma
siguiente, considerando el vector x 1 las variables de la primera etapa y x 2 las
variables de la segunda:
min(c1T x 1 + cT2 x 2 )
x1 ,x 2
A1x 1 = b1
(2.1)
B1x 1 +A2x 2 = b2
x 1, x2 ≥0
donde A1 ∈ ℝ m1×n1 y A2 ∈ ℝ m2×n2 . Las dimensiones de los demás vectores y

matrices se derivan de éstas. Por ser optimización determinista se suponen
conocidas las matrices y vectores. El tamaño del problema completo es
(m1 + m2 ) × (n1 + n2 ) . La estructura de la matriz de restricciones del problema
(denominada triangular inferior por bloques) se presenta en la figura 2.1.
A1
B1 A2
Figura 2.1 Estructura de la matriz de coeficientes de las restricciones en PL-2.
Los dos ejemplos del capítulo anterior corresponden a problemas lineales

bietapa, deben proporcionar decisiones antes de resolver la incertidumbre,
siendo la segunda etapa estocástica. Un problema de planificación estática (para
un momento fijo en el tiempo) se formula frecuentemente como minimización de
una función objetivo suma de costes totales de inversión y explotación sujeta a
restricciones propias de inversión y de explotación. Las decisiones de inversión
corresponden a la primera etapa y las de explotación a la segunda.
Los problemas de planificación dinámica deben proporcionar decisiones
óptimas para momentos discretos del tiempo. Éstos se conocen genéricamente
como problemas lineales multietapa PL-P. Aparecen en el estudio de procesos
dependientes del tiempo donde las actividades de un periodo están conectadas
con las del periodo anterior pero no con otros [Dantzig:63]. Para P etapas el
problema se puede formular como:
04/11/2012 15
2 OPTIMIZACIÓN LINEAL BIETAPA Y MULTIETAPA
P
min ∑ cTp x p
xp
p =1
Bp−1x p−1 + Ap x p = bp p = 1, …, P (2.2)

xp ≥ 0
B0 ≡ 0
m ×n
donde las matrices Ap ∈ ℝ p p y Bp y los vectores bp y cp son conocidos por
ser optimización determinista y x p son las variables. El problema completo tiene
P P
m = ∑ p =1 m p filas y n = ∑ p =1 n p columnas.
Característica de estos problemas es una estructura de la matriz de
coeficientes de las restricciones en escalera (o diagonal por bloques), tal como se
presenta en la figura 2.2.
A1
B1 A2
B2 A3
⋱
BP −1 AP
Figura 2.2 Estructura de la matriz de coeficientes de las restricciones en problemas lineales multietapa.
El tamaño del problema crece linealmente con el número de etapas. Aunque

para problemas pequeños o medianos se pueden utilizar técnicas convencionales
de programación lineal, esto resulta imposible para problemas de gran tamaño.
Por ello, para su resolución se emplean técnicas de descomposición que
aprovechan su especial estructura.
Un problema de gestión óptima de recursos para un conjunto de periodos se
formula como minimización de una función objetivo para el conjunto sujeta a
restricciones que acoplan periodos consecutivos y otras específicas de cada uno
de ellos.
16 04/11/2012
I.3 Técnicas de descomposición

La resolución de problemas lineales estocásticos de gran tamaño requiere de
técnicas de optimización y/o simulación muy eficientes. En este documento se
presentan algunas de ellas:
• descomposición
Las técnicas de descomposición resuelven problemas de muy gran tamaño
con una estructura especial, que se aprovecha desde un punto de vista
teórico y computacional, mediante la solución iterativa de otros problemas
de menor tamaño con estructura similar. Tiene sentido aplicarlas a un
problema cuya estructura específica permite identificar partes del mismo que
son fácilmente resolubles de modo individual. Constituye una forma flexible
y elegante de resolver problemas.
Los tipos de estructuras especiales pueden ser:
• Multidivisional: conjunto de subsistemas interrelacionados
Por ejemplo, el problema de flujo de cargas multisistema (España,
Francia y Portugal o el sistema eléctrico centroamericano) o la
programación semanal de un sistema de energía eléctrica con múltiples
grupos.
• Dinámica: gran número de restricciones y variables replicadas para cada
periodo
Por ejemplo, la coordinación hidrotérmica a lo largo de múltiples
periodos.
• Estocástica o combinatorial: instanciaciones de variables aleatorias
discretas
Por ejemplo, la planificación de la generación para múltiples escenarios
futuros.
• simulación y reducción de varianza

Las técnicas de simulación y reducción de varianza son necesarias cuando la
resolución explícita del número de problemas que surgen debido a la
aleatoriedad en los parámetros resulta inviable.
Una técnica de descomposición para la solución de problemas lineales

deterministas bietapa fue propuesta por Dantzig-Wolfe (DW) [Dantzig:60] en su
forma primal y por Benders (Bd) [Benders] en su forma dual. En la primera se
optimiza en sucesivas subregiones de la región factible hasta detectar el óptimo.
04/11/2012 17
3 TÉCNICAS DE DESCOMPOSICIÓN
En la segunda se optimiza sobre una región factible mayor e incrementando las

restricciones se aproxima a la solución óptima en la verdadera región factible.
La relajación lagrangiana (RL) [Geoffrion:74] es equivalente a la descomposición
de DW excepto que el maestro se resuelve en su forma dual. Para la solución de
problemas lineales deterministas multietapa se emplea la descomposición
anidada, cuya primera aplicación se cita en [Dantzig:63]. Geoffrion [Geoffrion:70,
Geoffrion:71] unifica los dos conceptos fundamentales utilizados en optimización
matemática a gran escala: manipulación del problema y estrategia de solución.
El primero significa expresar el problema de forma alternativa esencialmente
equivalente pero más fácilmente solucionable. Como ejemplos: proyección,
linealización interior o exterior y dualización. La manipulación del problema
tiene tres objetivos: inducir separación entre problemas, inducir linealidad en un
problema parcialmente no lineal y aislar estructuras especiales para utilizar
algoritmos más eficientes. El segundo supone la reducción del problema de
optimización a una secuencia de problemas relacionados más sencillos. Como
ejemplos: linealización por tramos y relajación. En [Lasdon] se explican también
algunas de estas técnicas. Otras técnicas de descomposición incluidas en este
marco son: descomposición en cruz [Aardal, Hilger, Kim, VanRoy],
descomposición de Kornai-Liptàk [Aardal], aproximación externa [Duran].
La conjunción de técnicas de descomposición y simulación fue propuesta por
Dantzig [Dantzig:87, Dantzig:90] y ha dado resultados espectaculares por los
tamaños de los problemas de planificación resueltos [Dantzig:91a, Dantzig:91c,
Entriken:90, Infanger, Infanger:94, Stanford].
Estas técnicas hacen que el uso de procesamiento paralelo o distribuido5 sea
muy conveniente [Dantzig:87, Dantzig:90]. Con ellas la solución repetida de los
diferentes problemas se puede efectuar en paralelo entre los diferentes
procesadores disponibles [Entriken:88, Entriken:89, Ruszczynski].
Actualmente existe una herramienta para la solución de problemas lineales
estocásticos bietapa de gran tamaño denominada DECIS de Gerd Infanger que
se puede ejecutar en un ordenador paralelo [Dantzig:91b] o convencional
[Infanger]. Las reducciones de tiempo con respecto a su solución en ordenador
convencional logradas para este tipo de problema son del orden de la mitad del
número de procesadores utilizados. Por ejemplo, con un problema concreto
resuelto con 64 procesadores en paralelo tomando 600 muestras en la simulación
5
Las técnicas de descomposición son de uso matemático general y tienen sentido
independientemente de la tecnología del ordenador/es utilizados en su resolución. Sin embargo,
dada su naturaleza resulta muy conveniente la utilización de cálculo en paralelo (un ordenador
con múltiples CPUs débil o fuertemente acopladas) o cálculo distribuido (múltiples ordenadores
trabajando en colaboración) para reducir los tiempos de cálculo de manera sustancial.
18 04/11/2012
el tiempo de cálculo se redujo en 37.5 veces. Este decremento depende del

número de muestras y de procesadores.
Otra aplicación denominada MSLiP permite la solución de problemas
lineales estocásticos multietapa de tamaño medio para su uso en un ordenador
convencional [Gassmann:90]. Existe también una versión estocástica OSL-SE de
IBM del optimizador de programación lineal entera mixta OSL. Existe otro
optimizador denominado SLP-IOR de Peter Kall de la Universidad de Zurich
con capacidad de resolver problemas lineales estocásticos. Más recientemente
OptiRisk Systems ha desarrollado el código SPInE (http://www.optirisk-
systems.com/spine.asp) y en Chalmers University otro denominado BNBS
(http://www.math.chalmers.se/~alten/bnbs/bnbs.html). En [Ermoliev] se citan
además otros códigos de programación lineal estocástica.
La formulación de problemas estocásticos se está empezando a realizar
también con lenguajes algebraicos de modelado como GAMS de GAMS
Development Corporation (www.gams.com) [Ramos:98], AMPL
(www.ampl.com) [Gassmann:95, Gassmann:96] o MPL/SPInE de Maximal
Software (www.maximalsoftware.com). Tanto DECIS como OSL-SE están
disponibles conectados con GAMS.
Las técnicas de descomposición se pueden clasificar, ver [Flippo], según sean
las variables (Bd) o las restricciones (DW o RL) lo que complica la resolución
del problema. La estructura habitual de la matriz de restricciones es diagonal
por bloques con variables, restricciones o ambas que complican, tal como aparece
respectivamente en la figura 3.1. En principio, serían bloques independientes
excepto por la presencia de estos bloques adicionales de complicación que unen
los demás.
Figura 3.1 Matrices de restricciones diagonales por bloques.
Por ejemplo, un problema multiperiodo con esta estructura en su matriz de

restricciones tiene variables que complican y se pueden reorganizar para inducir
separabilidad en las restricciones específicas de los periodos, figura 3.2.
04/11/2012 19
⋮ ⋮
Figura 3.2 Matriz de restricciones de un problema multiperiodo.
En el caso de un problema multidivisional existen restricciones que

complican y pueden “eliminarse” para inducir separabilidad.
Figura 3.3 Matriz de restricciones de un problema multidivisional.
Obsérvese que el problema dual del segundo tiene la estructura del primero y
viceversa.
También pueden clasificarse según cuál sea la información que se envía desde
el maestro al subproblema. La descomposición de Bd envía variables del primal
y la de DW o RL variables del dual.
Matemáticamente las técnicas de descomposición se pueden utilizar siempre,
es decir, son formas exactas de resolver un problema de optimización. Sin
embargo, algorítmicamente sólo tienen sentido cuando existe una ventaja
computacional en su uso. En muchos casos, esta ventaja se deriva de la
separabilidad del problema completo en subproblemas similares de tamaño
mucho menor. Por ejemplo, en el caso de planificación multisistema subyace que
las relaciones (número de restricciones que acoplan) entre sistemas en el primer
caso o entre grupos en el segundo son débiles mientras que las relaciones dentro
del sistema son las más numerosas.
Veamos un modelo de coordinación hidrotérmica para el análisis de la
explotación de un sistema eléctrico a medio plazo. Se trata de un modelo de
explotación multiperiodo con gestión de las aportaciones hidráulicas para el
20 04/11/2012
conjunto de periodos. Se define como periodo a la unidad relevante de tiempo

para la gestión hidráulica (por ejemplo, una semana o un mes o intervalos con
duraciones intermedias crecientes según se alejan del presente). La estructura de
la matriz de restricciones presenta este aspecto.
Figura 3.4 Matriz de restricciones del problema de coordinación hidrotérmica.
En la figura 3.5 se observa claramente la estructura en escalera de la matriz

de restricciones de un problema real con 12 etapas. El tamaño del problema es
aproximadamente de 16500 restricciones y 25000 variables.
Figura 3.5 Estructura de la matriz de restricciones.
El planteamiento matemático para su resolución se puede hacer utilizando la

descomposición de Bd, de DW o descomposición anidada de Benders.
Alternativamente, este problema también puede resolverse por optimización
dinámica.
04/11/2012 21
Bajo el punto de vista de descomposición de Bd las variables de producción

total de cada hidráulica en cada periodo son las que complican. Por
consiguiente, el maestro elabora propuestas para estas variables y los
subproblemas son problemas de optimización intraperiodo exclusivamente. Cada
subproblema envía el valor (la derivada) de la producción total en cada periodo
de cada hidráulica (valor del agua) al maestro y recibe como dato la producción
hidráulica para cada periodo.
Figura 3.6 Matriz de restricciones del maestro y subproblema de Benders.
Bajo el prisma de descomposición de DW las restricciones de movimiento de

las reservas en los embalses para todos los periodos (pero expresadas entre
periodos consecutivos) complican la resolución. Cada subproblema hace una
propuesta de producción hidráulica para cada periodo y, por consiguiente, de
nivel final de reservas y el problema maestro, observando las restricciones de
gestión de los embalses, envía una señal económica de la restricción de gestión
hidráulica. Cada subproblema independientemente incluye esta señal en el coste
de explotación del periodo y envía al maestro un nuevo valor de la producción
hidráulica. La potencia hidráulica es una variable para el subproblema.
22 04/11/2012
Figura 3.7 Matriz de restricciones del maestro y subproblema de Dantzig-Wolfe.
La descomposición anidada de Benders es la aplicación recursiva de Benders

a los subproblemas que aparecen en el método de descomposición de Benders.
Cada subproblema padre envía las propuestas de niveles finales de los embalses
al subproblema hijo y cada subproblema hijo devuelve la variable dual asociada
a esa propuesta.
Figura 3.8 Matriz de restricciones del maestro y subproblema de Benders anidado.
La elección entre un algoritmo de resolución u otro se debe basar en dos

criterios fundamentales. Desde un punto de vista de la ingeniería influyen
• su adaptación al uso del modelo, por ejemplo en contextos

descentralizados,
• el significado físico o la manera más natural de manejar la información
que se intercambian maestro y subproblemas
Desde un punto de vista matemático se analizará:
• la idoneidad algorítmica en función de los resultados (tiempo de

resolución y número de iteraciones previsibles)
• tamaños de los respectivos maestros y subproblemas y adaptación a los
principios básicos de las descomposiciones
• potenciales problemas de infactibilidad
A continuación se explica cada uno de estos métodos, empezando por el

método de descomposición de Benders por ser el más utilizado actualmente y el
más intuitivo de entender.
04/11/2012 23
I.4 Descomposición de Benders

I.4.1 Descripción general
El método de descomposición de Benders (Bd) [Benders, Geoffrion:70,
Geoffrion:71, VanSlyke] recibe también el nombre de descomposición primal
porque el problema maestro fija variables del primal, descomposición en L
porque se aplica a problemas con matriz de restricciones con dicha forma,
descomposición por recursos porque el maestro asigna directamente las
decisiones sobre los recursos al subproblema.
Descompone el problema lineal bietapa PL-2 en un problema maestro y un
subproblema. El maestro representa la primera etapa más las condiciones
necesarias, denominadas cortes, derivadas de la segunda etapa. El subproblema
representa la segunda etapa para decisiones conocidas de la primera etapa. El
algoritmo es iterativo y alterna entre la solución del maestro y del subproblema.
Este método se utiliza cuando las variables x 1 de la primera etapa complican
la solución del problema. Es decir, ésta resulta mucho más sencilla cuando
aquéllas se fijan temporalmente. Implícitamente se supone que el número de
variables x 1 es reducido frente al número total de variables ( n1 ≪ n2 ). El
número de iteraciones de descomposición está implícitamente relacionado con el
número de variables de complicación. También tiene sentido su utilización
cuando maestro y subproblemas son de diferente naturaleza. Por ejemplo,
cuando el maestro es un problema MIP o cuando el subproblema es un
problema NLP. En su presentación se utilizarán problemas formulados
linealmente aunque este método ha sido generalizado y permite la minimización
de una función objetivo convexa con región factible formando un poliedro
convexo, ver [Geoffrion:72]. En [Paules] se encuentra una discusión de la
implantación del método de descomposición generalizada de Bd aplicado a un
problema de programación mixta no lineal.
El problema lineal bietapa PL-2 (2.1) se puede interpretar también de esta
manera:
min c1T x 1 + θ2 (x 1 )
x1 ,θ2 (x1 )
A1x 1 = b1 (4.1)
x1 ≥ 0
donde la función de recursos, θ2 (x 1 ) ∈ ℝ , es una función poligonal convexa
dependiente de x 1 . Representa la función objetivo de la segunda etapa como
función de las decisiones de la primera etapa y tiene la siguiente expresión:
04/11/2012 25
4 DESCOMPOSICIÓN DE BENDERS
θ2 (x 1 ) = min cT2 x 2
x2
A2x 2 = b2 − B1x 1 : π2 (4.2)

x2 ≥ 0
siendo π2 las variables duales (o precios sombra) de las restricciones.

Al problema (4.1) se le conoce como problema maestro y al (4.2) como
subproblema en la descomposición de Bd. Se supone que el subproblema es
factible para cualquier valor de x 1 (es decir, con recurso completo) o se
introducen penalizaciones por las infactibilidades. Esta condición es frecuente en
muchos problemas de planificación. Si no se cumple se verá más adelante cuál es
la modificación del algoritmo.
Expresando en su forma dual se obtiene:
θ2 (x 1 ) = max(b2 − B1x 1 )T π2
π2
(4.3)
T
A π2 ≤ c2
2
Sea Π = {π21, π22 , …, π2ν } el conjunto finito de vértices del poliedro convexo
(politopo) definido por la región factible A2T π2 ≤ c2 . Obsérvese que la región
factible del problema dual no depende del valor de x 1 . Se sabe que la solución
óptima de un problema lineal reside en un vértice, por lo tanto, el problema se
podría resolver por enumeración de todos ellos:
θ2 (x 1 ) = max {(b2 − B1x 1 )T π2l } l = 1, …, ν (4.4)
De la ecuación anterior se deriva que la función de recursos θ2 (x 1 ) es una

función poligonal de las variables x 1 . Si se expresa como problema lineal se
obtiene:
θ2 (x 1 ) = min θ2
θ2
θ2 ≥ (b2 − B1x 1 )T π21

(4.5)
⋮
θ2 ≥ (b2 − B1x 1 )T π2ν
donde θ2 ∈ ℝ y las restricciones se denominan cortes o planos de corte o

hiperplanos soporte o tangentes a la función objetivo en cada punto x 1 .
Constituyen una linealización (aproximación) exterior convexa de la función de
recursos.
Entonces, el problema original PL-2 se puede expresar como:
26 04/11/2012
min c1T x 1 + θ2
x1 ,θ
A1x 1 = b1
θ2 ≥ (b2 − B1x 1 )T π21
(4.6)
⋮
θ2 ≥ (b2 − B1x 1 )T π2ν
x1 ≥ 0
Esta formulación se denomina problema maestro completo, ya que contiene
todos los cortes posibles. Presenta todas las restricciones de la primera etapa
más todas las condiciones necesarias derivadas de la segunda etapa. Obsérvese
que la variable θ2 es libre.
Desde el punto de vista práctico, la resolución del problema maestro
completo implica disponer de forma explícita de todos los cortes de Benders, lo
cual es prácticamente imposible en problemas de tamaño realista. En lugar de
incluir todos los cortes (lo que implicaría disponer de forma explícita de la
función de recursos θ2 (x 1 ) ), el algoritmo introduce uno en cada iteración. De
esta forma, el problema maestro relajado para la iteración j se define como:
Maestro Benders
min cT1 x 1 + θ2
x1 ,θ2
A1x 1 = b1
(4.7)
lT lT
π B1x 1 + θ2 ≥ π b
2 2 2 l = 1, …, j
x1 ≥ 0
siendo θ2 ∈ ℝ , l el índice de iteraciones, x 1 y θ2 variarán en cada iteración6.
En cada iteración del algoritmo de Benders, la variable dual generada en el
subproblema es distinta del conjunto de variables duales generadas con
anterioridad por el algoritmo [VanSlyke]. Dado que el conjunto de posibles
valores duales vértices es finito, el número de cortes que se pueden obtener es
también finito luego el número de iteraciones está también limitado. El
algoritmo de descomposición de Bd converge en un número finito de iteraciones.
Se entiende por corte válido aquél que es externo a la función de recursos
aunque no necesariamente tangente. Esto implica que no es necesario resolver el
subproblema hasta optimalidad para garantizar que los cortes son válidos.
De cara a una implantación eficiente del algoritmo de descomposición, los
cortes de Benders aceptan la siguiente formulación como linealización de la
6
Se entiende por iteración un ciclo problema maestro - subproblema, en este orden.
04/11/2012 27
función de recurso en torno a los valores de las variables de la primera etapa y

de los valores de la variable θ2 obtenidos en cada iteración.
θ2 ≥ π2jT (b2 − B1x 1 ) = π2jT (b2 − B1x 1 + B1x 1j − B1x 1j ) =

(4.8)
= π2jT b2 − B1x 1j − B1 (x 1 − x 1j ) = π2jT b2 − B1x 1j  + π2jT −B1(x 1 − x 1j )
siendo x 1j y f2j = π2jT b2 − B1x 1j  los valores de las variables de la primera etapa
y el de la función objetivo de la segunda para la iteración j . Luego el conjunto
de cortes hasta la iteración j , l = 1, …, j , también se expresa como
θ2 − f2l ≥ π2lT B1 (x 1l − x 1 ) (4.9)
θ2 + π2lT B1x 1 ≥ f2l + π2lT B1x 1l (4.10)

Esta expresión indica que π2jT B1 es un subgradiente7 del valor de la función
de recursos θ2 (x 1 ) para la propuesta x 1j del maestro.
El problema maestro tiene ahora esta expresión
Maestro Benders
min cT1 x 1 + θ2
x1 ,θ2
A1x 1 = b1
(4.11)
π2lT B1x 1 + θ2 ≥ f2l + π2lT B1x 1l l = 1, …, j
x1 ≥ 0
El subproblema para cada iteración j se formula como:
Subproblema Benders
f2j = min cT2 x 2

x2
A2x 2 = b2 − B1x 1j : π2j (4.12)

x2 ≥ 0
I.4.2 Cortes de infactibilidad

La descripción anterior del algoritmo ha supuesto que el subproblema es
factible y acotado para cualquier propuesta del problema maestro. Esta
7
Se utiliza el subgradiente como una generalización del concepto de gradiente cuando la
función es convexa pero no es diferenciable en todos sus puntos. π2lT B1 es un subgradiente
porque la ecuación θ2 − f2l ≥ π2lT B1 (x 1l − x 1 ) es un plano de corte o hiperplano soporte del
epigrafo de la función.
28 04/11/2012
hipótesis, conocida en la literatura como recurso parcialmente completo8, no

suele satisfacerse en la práctica y el algoritmo de descomposición es modificado
cuando esto ocurre. La modificación del algoritmo consiste en la construcción de
otro tipo de corte, corte de infactibilidad, que elimina la solución propuesta en el
problema maestro. La construcción de este corte se comenta a continuación.
Si el subproblema es factible para un valor de x 1 los precios sombra de sus
restricciones son los obtenidos para su solución óptima y los cortes formados en
el maestro se denominan cortes de optimalidad. Si el subproblema es factible
pero no acotado el problema lineal bietapa es no acotado. Si el subproblema es
infactible entonces los precios sombra son los derivados de la fase I del simplex
[Dantzig:63, Dantzig:91d] y los cortes se denominan cortes de infactibilidad. En
esta fase la función objetivo del subproblema es la suma de infactibilidades:
min eT v + + eT v −
x 2 ,v + ,v −
A2x 2 + Iv + − Iv − = b2 − B1x 1 : π2 (4.13)

x 2, v +, v − ≥ 0
siendo v +, v − ∈ ℝm2 , eT = (1 ⋯ 1) , I la matriz identidad m2 × m2 y π2 las
variables duales de las restricciones para la solución óptima.
Recordamos que el lema de Farkas para el problema
min cT x
Ax = b (4.14)
x ≥0
dice que exactamente uno de los sistemas siguientes tiene solución.
• Sistema 1: AT x = b y x ≥ 0 para algún x ∈ ℝ n

• Sistema 2: AT y ≤ 0 y bT y > 0 para algún y ∈ ℝ m
Luego el problema (4.12) será factible si para todo vector de variables duales
π2 (o rayos extremos) que cumple A2T π2 ≤ 0 se cumple (b2 − B1x 1 )T π2 ≤ 0 . En
este caso, el subproblema se formula como
max(b2 − B1x 1 )T π2
π2
(4.15)
T
A π2 ≤ 0
2
8
En un problema bietapa por recurso completo se entiende que el subproblema de la
segunda etapa es siempre factible para cualquier valor de las variables de la primera etapa.
04/11/2012 29
Si la función objetivo del subproblema alcanza un valor estrictamente

positivo entonces se genera un corte de infactibilidad que excluye los valores de
x 1 que verifican
(b2 − B1x 1 )T π2j ≤ 0 (4.16)

Sin embargo, el problema (4.15) puede ser no acotado por ser la región
factible un cono luego para ello se introducen unas cotas a las variables (se
inscriben en el cubo unidad) formulando ahora este problema
max(b2 − B1x 1 )T π2
π2
T
A π2 ≤ 0
2 (4.17)
−1 ≤ π2 ≤ 1
Si obtenemos el problema dual de éste observamos que corresponde
precisamente al problema (4.13) de minimización de infactibilidades.
El problema maestro considerando ambos tipos de cortes, de optimalidad e
infactibilidad, se formula de la siguiente manera:
min c1T x 1 + θ2
x1 ,θ2
A1x 1 = b1
(4.18)
π2lT B1x 1 + δ1l θ2 ≥ π2lTb2 l = 1, …, j
x1 ≥ 0
o bien
min cT1 x 1 + θ2
x1 ,θ2
A1x 1 = b1
(4.19)
π2lT B1x 1 + δ1l θ2 ≥ f2l + π2lT B1x 1l l = 1, …, j
x1 ≥ 0
siendo δ1l = 1 para los cortes de optimalidad y δ1l = 0 para los de infactibilidad.
Los cortes de infactibilidad se pueden considerar como que fueran cortes de
optimalidad con pendiente infinita. Eliminan las propuestas del maestro que
hacen infactible el subproblema pero conservan las que no lo hacen infactible.
Normalmente es más conveniente, desde el punto de vista de cálculo,
formular el problema maestro con penalización por infactibilidad en el
subproblema que resolver éste y enviar los cortes de infactibilidad al maestro.
30 04/11/2012
I.4.3 Algoritmo
En cada iteración se modifica la región factible del maestro (aumenta en uno
el número de restricciones) y las cotas de las restricciones del subproblema. El
tamaño del problema lineal bietapa PL-2 era (m1 + m2 ) × (n1 + n2 ) . El tamaño
del problema maestro es (m1 + j ) × (n1 + 1) y el del subproblema m2 × n2 .
En una iteración del algoritmo se resuelve el problema maestro relajado y se
pasa el valor del vector x 1j al subproblema. Éste optimiza x 2 con los recursos
(b2 − B1x 1j ) y pasa las variables duales π2j (precios sombra de las restricciones)
de nuevo al maestro. Con éstas se forma un corte que se añade
consecutivamente a las restricciones del maestro. Esquemáticamente el
algoritmo se representa en la figura 4.1.
PROBLEMA MAESTRO
| ↑
j
propuesta x 1 precios sombra π2j
↓ |
SUBPROBLEMA
Figura 4.1 Método de descomposición de Benders.
Como la información relativa a la segunda etapa (cortes) aumenta con el

número de iteraciones, el maestro hace cada vez decisiones más próximas al
óptimo. Para cada iteración una cota superior z del valor óptimo de la función
objetivo del problema original PL-2 viene dada por (c1T x 1j + cT2 x 2j ) siendo x 1j y
x 2j soluciones factibles en maestro y subproblema en esa iteración. En cada
iteración, el valor obtenido por la función objetivo del problema maestro
relajado z , dada por (c1T x 1j + θ2j ) , es una cota inferior del problema completo
PL-2. La condición de convergencia para la terminación del algoritmo es la
coincidencia de ambas cotas con una tolerancia relativa ε (por ejemplo, 10−4 ).
z −z cT2 x 2j − θ2j
= T j ≤ε (4.20)
z c1 x 1 + cT2 x 2j
Se puede demostrar [Flippo] que si una solución (x 1j , θ2j ) del problema

maestro es ε1 -óptima, una solución x 2j del subproblema (primal) es ε2 -óptima y
una solución π2j del subproblema (dual) es ε3 -óptima, entonces (x 1j , x 2j ) es una
solución es ε ′ = (ε1 + ε2 + ε3 ) -óptima para el problema completo.
Cada corte tiene el efecto de eliminar la solución óptima de la iteración
anterior de la nueva región factible del maestro. Por esta razón, no es posible
generar más de una vez el mismo corte. Si en una iteración el algoritmo no ha
04/11/2012 31
convergido es porque cT2 x 2j − f2j > ε ′ . En la siguiente iteración el corte que se

introduce es éste
θ2 ≥ f2j + π2jT B1(x 1j − x 1 )
θ2 ≥ cT2 x 2j + π2jT B1(x 1j − x 1 ) (4.21)
θ2 ≥ θ2j + ε ′ + π2jT B1(x 1j − x 1 )
donde se observa que efectivamente el corte resulta infactible para la solución de
la iteración anterior. Por contradicción, se puede demostrar que si se pasa dos
veces la misma información desde el subproblema (o desde el problema maestro)
el algoritmo termina.
La sucesión de estas cotas inferiores es monótona creciente dado que en cada
iteración el problema maestro relajado contiene mayor número de restricciones,
mientras que la cota superior no es necesariamente decreciente. Por esta razón,
se toma como cota superior para evaluar convergencia el mínimo de todas las
cotas superiores previas. Por consiguiente, la solución óptima no corresponde
necesariamente a la de la última iteración.
Para obtener el valor del vector x 10 de las variables de la primera etapa a
enviar al subproblema en la primera iteración se pueden hacer estimaciones
razonables, si se conoce la naturaleza del problema, o solucionar el maestro sin
cortes, lo que equivale a fijar el valor de la variable de recurso a cero, θ2 = 0 .
En cada iteración del maestro se dispone de una base que es infactible sólo
por una variable básica, la variable de holgura del nuevo corte. Desde el punto
de vista del dual del maestro los cortes aparecen como columnas. Al añadir una
nueva columna la solución previa sigue siendo factible y la nueva solución
óptima puede obtenerse en pocas iteraciones del simplex. Por ello es
conveniente, teóricamente, resolver el maestro mediante el simplex dual.
Solamente los cortes activos en una iteración son necesarios en la siguiente, los
demás pueden ser eliminados, ya que la función objetivo del problema maestro
es monótona creciente al introducir nuevos cortes. Esta medida ayuda a reducir
el tiempo de cálculo [Birge:88].
Cada subproblema sólo cambia las cotas de las restricciones en cada
iteración. Por esta razón, suponiendo que ninguna solución del maestro ocasiona
infactibilidad en el subproblema, es conveniente resolverlo mediante el método
simplex primal (siempre que su tamaño lo aconseje).
El método de descomposición de Bd tiene la ventaja de disponer en cualquier
iteración de una solución factible, aunque no sea la óptima. Es decir, aunque no
se llegue a convergencia se puede disponer de una buena solución cuasióptima.
Esquemáticamente el algoritmo se formula a continuación:
32 04/11/2012
1. Inicialización: j = 0 , z = ∞ , z = −∞ , ε = 10−4
2. Resolución del problema maestro
min cT1 x 1 + θ2
x1 ,θ2
A1x 1 = b1
(4.22)
π2lT B1x 1 + δ1l θ2 ≥ f2l + π2lT B1x 1l l = 1, …, j
x1 ≥ 0
Obtener la solución (x 1j , θ2j ) y evaluar la cota inferior z = c1T x 1j + θ2j .
Mientras no se haya generado ningún corte de optimalidad, se fija el valor de
la variable de recurso θ2 a cero, pues en otro caso el problema maestro es no
acotado. Una vez obtenido algún corte de infactibilidad, esta variable pasa a
ser libre.
3. Resolución del subproblema de suma de infactibilidades
f2j = min eT v + + eT v −
x 2 ,v + ,v −
A2x 2 + Iv + − Iv − = b2 − B1x 1j : π2j (4.23)

x 2 , v +, v − ≥ 0
Si f2j ≥ 0 , obtener π2j , formar un corte de infactibilidad y añadirlo al

problema maestro, incrementar el número de iteraciones j = j + 1 e ir al
paso 2.
Si f2j = 0 , ir al paso 4.
4. Resolución del subproblema
f2j = min cT2 x 2
x2
A2x 2 = b2 − B1x 1j : π2j (4.24)

x2 ≥ 0
Obtener x 2j y actualizar cota superior z = c1T x 1j + cT2 x 2j .

5. Regla de parada
z −z cT2 x 2j − θ2j
Si = T j T j
≤ ε detener el algoritmo. En otro caso, obtener π2j ,
z c1 x 1 + c2 x 2
formar un corte de optimalidad y añadirlo al problema maestro, incrementar
el número de iteraciones j = j + 1 e ir al paso 2.
La siguiente sección presenta con un ejemplo los pasos del algoritmo de

descomposición de Bd.
04/11/2012 33
I.4.4 Caso ejemplo

Veamos gráficamente el funcionamiento algoritmo. Se desea resolver el
siguiente problema partiendo de una solución inicial de (x , y ) = (0, 4) .
min− 2x − y
x ,y
x ≤4
x +y ≤5
2x +3y ≤ 12
x, y ≥0
Suponemos que x es la variable de la primera etapa e y variable de la

segunda. El primer subproblema es el siguiente
θ(x ) = min− y
y
y ≤5−0
3y ≤ 12 − 2 ⋅ 0
y≥0
 ⋅ 
que da lugar a f = −4 y π = 
1 1
 . Luego la cota inferior del problema es
−1 3
z = −∞ y la cota superior es z = −4 .
El problema maestro corresponde a éste
min− 2x + θ
x ,θ
x ≤4
1
θ − (−4) ≥ ⋅ −1 3   (0 − x )
( )
2
x ≥0
manipulando el corte queda
min− 2x + θ
x ,θ
x ≤4
2
θ − x ≥ −4
3
x ≥0
que da como resultado x = 4 y θ = − 4 3 luego la cota inferior del problema es
z = − 28 3 . Formulando de nuevo el subproblema para x = 4
34 04/11/2012
θ(x ) = min− y
y
y ≤5−4 =1
3y ≤ 12 − 2 ⋅ 4 = 4
y ≥0
−1
que da lugar a f 2 = −1 y π 2 =   . Luego la cota superior es
 ⋅ 
z = min (−9, −4) = −9 .
El nuevo problema maestro es
min− 2x + θ
x ,θ
x ≤4
2
θ − x ≥ −4
3
θ − x ≥ −5
x ≥0
que da como resultado x = 4 y θ = −1 luego la cota inferior del problema es
z = −9 . Una vez formulado y resuelto de nuevo el subproblema se obtiene
z = min (−9, −9) = −9 y, por consiguiente, converge el algoritmo.
I.4.5 Algoritmo escrito en GAMS

Veamos con otro ejemplo los pasos del algoritmo previo. Se desea resolver el
siguiente problema lineal:
min− x − 2y − 2u − 3v
x ,y ,u ,v
x +y ≤ 600
x −2y ≤0
x +y +u +v ≤ 1000
x +u ≤ 500
−2u +v ≤0
x, y, u, v ≥0
04/11/2012 35
1000
−1 −2 600   1 1 

entonces, c1 =   , c2 =   , b1 =   , b2 =  500  , A1 =   ,
−2 −3  0  
 0 
 1 −2
 
1 1  1 1
   
 
B1 = 1 0 , A2 =  1 0
   
0 0 −2 1
   
Se supone una solución inicial factible para el problema maestro
3.000
x 1 =   .
3.000
Se calcula la cota inferior, siendo z = −9.000 . Al introducir en el
⌢
331.3
  y
subproblema estos valores se obtiene como solución x 2 = 
662.6⌢ 
 
 ⌢
−2.666
 
π2 =  0.000 

 
−0.333⌢ 
  ⌢
Se calcula la cota superior, siendo z = −2659.6 . Como la diferencia entre
ambas cotas es elevada se continúa iterando. Con los precios sombra se forma
un corte en el maestro y se resuelve éste de nuevo. La nueva solución obtenida
0.000
es x 1 =   .
0.000
⌢
Se calcula de nuevo la cota inferior, siendo z = −2666.6 . Al introducir en el
⌢
333.3
  y
subproblema estos valores se obtiene como solución x 2 = 
666.6⌢ 
 
 ⌢
−2.666
 
π2 =  0.000  .
 
−0.333⌢ 
  ⌢
Se calcula la cota superior en esta iteración, siendo z = −2666.6 . La
diferencia entre ambas cotas es ahora nula y, por lo tanto, acaba el algoritmo.
⌢ 0.000
La función objetivo es z * = −2666.6 y la solución óptima es x 1* =   y
0.000
⌢
333.3
  .
x 2* = 
666.6⌢ 
 
36 04/11/2012
A continuación se presenta el método de descomposición de Bd.
$title Descomposición de Benders (Bd)

sets
L máximo número de iteraciones / iter-1 * iter-5 /
J(l) iteración actual
N1 número variables primera etapa
M2 número restricciones segunda etapa
scalar
TOL tolerancia relativa / 1e-6 /
Z_INF cota inferior / -inf /
Z_SUP cota superior / inf /
parameters
PI2(l,m2) variables duales restricciones segunda etapa
DELTA1(l) coeficiente de cortes de infactibilidad (subproblema infactible)
X1_J(n1,l) valores variables primera etapa hasta iteración j
Z2_J(l) valores func obj segunda etapa hasta iteración j
* comienzo datos del problema
sets
M1 número restricciones primera etapa / r1-1 * r1-2 /
N1 número variables primera etapa / x1-1 * x1-2 /
M2 número restricciones segunda etapa / r2-1 * r2-3 /
N2 número variables segunda etapa / x2-1 * x2-2 /
parameters
C1(n1) coeficientes función objetivo primera etapa
/ x1-1 -1
x1-2 -2 /
C2(n2) coeficientes función objetivo segunda etapa
/ x2-1 -2
x2-2 -3 /
B1(m1) cotas restricciones primera etapa
/ r1-1 600
r1-2 0 /
B2(m2) cotas restricciones segunda etapa
/ r2-1 1000
r2-2 500
r2-3 0 /
table A1(m1,n1) matriz de restricciones primera etapa
x1-1 x1-2
r1-1 1 1
r1-2 1 -2
table BM1(m2,n1) matriz de restricciones segunda etapa

x1-1 x1-2
r2-1 1 1
r2-2 1 0
r2-3 0 0
table A2(m2,n2) matriz de restricciones segunda etapa

x2-1 x2-2
r2-1 1 1
r2-2 1 0
r2-3 -2 1
* fin datos del problema

positive variables
X1(n1) variables primera etapa
X2(n2) variables segunda etapa
variables
TT función de recursos
Z1 función objetivo primera etapa
Z2 función objetivo segunda etapa
equations
FOC función objetivo problema completo
FO1 función objetivo primera etapa
FO2 función objetivo segunda etapa
R1(m1) restricciones primera etapa
R2(m2) restricciones segunda etapa
04/11/2012 37
CORTESA(l) cortes de Benders

CORTESB(l) cortes de Benders ;
FOC .. Z1 =E= sum
sum(n1, C1(n1)*X1(n1)) + sum
sum(n2, C2(n2)*X2(n2)) ;
FO1 .. Z1 =E= sum

sum(n1, C1(n1)*X1(n1)) + TT ;
FO2 .. Z2 =E= sum(n2, C2(n2)*X2(n2)) ;
sum
R1(m1) .. sum
sum(n1, A1(m1,n1)*X1(n1)) =L= B1(m1) ;
R2(m2) .. sum
sum(n1, BM1(m2,n1)*X1(n1)) + sum
sum(n2, A2(m2,n2)*X2(n2)) =L= B2(m2) ;
CORTESA(j) .. sum
sum(m2, PI2(j,m2)* sum
sum(n1, BM1(m2,n1)*X1(n1))) + DELTA1(j) * TT =G=
sum(m2, PI2(j,m2)*B2(m2)) ;
sum
CORTESB(j) .. DELTA1(j) * TT =G= Z2_J(j) +
sum(m2, PI2(j,m2)*sum
sum sum(n1, BM1(m2,n1)*(X1_J(n1,j) - X1(n1)))) ;
sum
*model MAESTRO / FO1, R1, CORTESA /
model MAESTRO / FO1, R1, CORTESB /
model SUB / FO2, R2 /
model COMPLETO / FOC, R1, R2 /

file COPT / cplex.opt / ;
put COPT putclose 'scaind -1' / 'lpmethod 1' / 'preind 0' / 'epopt 1e-9' / 'eprhs 1e-9'
SUB.OptFile = 1 ;
* inicialización de parámetros del problema
J(l) = NO ;
TT.FX = 0 ;
DELTA1(l) = 0 ;
PI2(l,m2) = 0 ;
X1_J(n1,l) = 0 ;
Z2_J(l) = 0 ;
loop(l
loop $(ABS(1-Z_INF/Z_SUP) > TOL),
if (ord
ord(l)
ord = 1,
* alternativamente en la primera iteración se puede dar una solución
X1.L(n1) = 3 ;
Z1.L = sum(n1,
sum C1(n1)*X1.L(n1)) ;
* o resolver directamente el problema maestro
* solve MAESTRO USING LP MINIMIZING Z1 ;
else
solve MAESTRO USING LP MINIMIZING Z1 ;
) ;
X1_J(n1,l) = X1.L(n1) ;
* fijación de la variable de la primera etapa X1
X1.FX(n1) = X1.L(n1) ;
solve SUB USING LP MINIMIZING Z2 ;
* adquisición de los parámetros para formar un nuevo corte
if(SUB.MODELSTAT
if = 4,
* subproblema infactible
DELTA1(l) = 0 ;
Z2_J(l) = SUB.SUMINFES ;
else
* actualización de la cota inferior y superior
Z_INF = Z1.L ;
Z_SUP = MIN(Z_SUP, Z1.L - TT.L + Z2.L) ;
DELTA1(l) = 1 ;
Z2_J(l) = Z2.L;
TT.LO = -inf
inf ; TT.UP = inf ;
) ;
display Z_INF, Z_SUP ;
PI2(l,m2) = R2.M(m2) ;
X1.LO(n1) = 0 ; X1.UP(n1) = inf ;
* Incremento del conjunto de cortes
J(l) = yes ;
38 04/11/2012
) ;
abort $(ABS(1-Z_INF/Z_SUP) > TOL) 'Máximo número de iteraciones alcanzado' ;
* resolución del problema completo
solve COMPLETO USING LP MINIMIZING Z1
GAMS Rev 142 Intel /MS Window 04/14/05 07:13:26 Page 1

Descomposición de Benders (Bd)
C o m p i l a t i o n
COMPILATION TIME = 0.000 SECONDS 3.2 Mb WIN217-142 Apr 05, 2005
Equation Listing SOLVE SUB Using LP From line 135
---- FO2 =E= función objetivo segunda etapa
FO2.. 2*X2(x2-1) + 3*X2(x2-2) + Z2 =E= 0 ; (LHS = 0)
---- R2 =L= restricciones segunda etapa
R2(r2-1).. X1(x1-1) + X1(x1-2) + X2(x2-1) + X2(x2-2) =L= 1000 ; (LHS = 6)
R2(r2-2).. X1(x1-1) + X2(x2-1) =L= 500 ; (LHS = 3)
R2(r2-3).. - 2*X2(x2-1) + X2(x2-2) =L= 0 ; (LHS = 0)
Model Statistics SOLVE SUB Using LP From line 135
LOOPS L iter-1
MODEL STATISTICS
BLOCKS OF EQUATIONS 2 SINGLE EQUATIONS 4
BLOCKS OF VARIABLES 3 SINGLE VARIABLES 5
NON ZERO ELEMENTS 11
GENERATION TIME = 0.030 SECONDS 4.0 Mb WIN217-142 Apr 05, 2005
EXECUTION TIME = 0.050 SECONDS 4.0 Mb WIN217-142 Apr 05, 2005
Solution Report SOLVE SUB Using LP From line 135
L O O P S L iter-1
S O L V E S U M M A R Y
MODEL SUB OBJECTIVE Z2
TYPE LP DIRECTION MINIMIZE
SOLVER CPLEX FROM LINE 135
**** SOLVER STATUS 1 NORMAL COMPLETION
**** MODEL STATUS 1 OPTIMAL
**** OBJECTIVE VALUE -2650.6667
RESOURCE USAGE, LIMIT 0.010 1000.000
ITERATION COUNT, LIMIT 2 10000
GAMS/Cplex Apr 1, 2005 WIN.CP.CP 21.7 028.031.041.VIS For Cplex 9.0
Cplex 9.0.2, GAMS Link 28
User supplied options:
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
Optimal solution found.
Objective : -2650.666667
LOWER LEVEL UPPER MARGINAL
---- EQU FO2 . . . 1.0000
---- EQU R2 restricciones segunda etapa
r2-1 -INF 1000.0000 1000.0000 -2.6667
r2-2 -INF 334.3333 500.0000 .
r2-3 -INF . . -0.3333
---- VAR X1 variables primera etapa
x1-1 3.0000 3.0000 3.0000 2.6667
x1-2 3.0000 3.0000 3.0000 2.6667
---- VAR X2 variables segunda etapa
x2-1 . 331.3333 +INF .
x2-2 . 662.6667 +INF .
---- VAR Z2 -INF -2650.6667 +INF .
04/11/2012 39
**** REPORT SUMMARY : 0 NONOPT

0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
---- 153 PARAMETER Z_INF = -9.000 cota inferior
PARAMETER Z_SUP = -2659.667 cota superior
Equation Listing SOLVE MAESTRO Using LP From line 128
---- FO1 =E= función objetivo primera etapa
FO1.. X1(x1-1) + 2*X1(x1-2) - TT + Z1 =E= 0 ; (LHS = 0)
---- R1 =L= restricciones primera etapa
R1(r1-1).. X1(x1-1) + X1(x1-2) =L= 600 ; (LHS = 6)
R1(r1-2).. X1(x1-1) - 2*X1(x1-2) =L= 0 ; (LHS = -3)
---- CORTESB =G= cortes de Benders
CORTESB(iter-1).. - 2.66666666666667*X1(x1-1) - 2.66666666666667*X1(x1-2) + TT =G= -2666.66666666667 ; (LHS = -16)
Model Statistics SOLVE MAESTRO Using LP From line 128
LOOPS L iter-2
MODEL STATISTICS
Solution Report SOLVE MAESTRO Using LP From line 128
L O O P S L iter-2
MODEL MAESTRO OBJECTIVE Z1
Objective : -2666.666667
---- EQU FO1 . . . 1.0000
FO1 función objetivo primera etapa
---- EQU R1 restricciones primera etapa
r1-1 -INF . 600.0000 .
r1-2 -INF . . .
---- EQU CORTESB cortes de Benders
iter-1 -2666.6667 -2666.6667 +INF 1.0000
x1-1 . . +INF 1.6667
x1-2 . . +INF 0.6667
---- VAR TT -INF -2666.6667 +INF .
---- VAR Z1 -INF -2666.6667 +INF .
TT función de recursos
0 INFEASIBLE
0 UNBOUNDED
---- FO2 =E= función objetivo segunda etapa
FO2.. 2*X2(x2-1) + 3*X2(x2-2) + Z2 =E= 0 ; (LHS = 0)
40 04/11/2012
R2(r2-1).. X1(x1-1) + X1(x1-2) + X2(x2-1) + X2(x2-2) =L= 1000 ; (LHS = 994)

R2(r2-2).. X1(x1-1) + X2(x2-1) =L= 500 ; (LHS = 331.333333333333)
R2(r2-3).. - 2*X2(x2-1) + X2(x2-2) =L= 0 ; (LHS = 0)
LOOPS L iter-2
MODEL STATISTICS
L O O P S L iter-2
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
Objective : -2666.666667
---- EQU FO2 . . . 1.0000
r2-1 -INF 1000.0000 1000.0000 -2.6667
r2-2 -INF 333.3333 500.0000 .
r2-3 -INF . . -0.3333
x1-1 . . . 2.6667
x1-2 . . . 2.6667
x2-1 . 333.3333 +INF .
x2-2 . 666.6667 +INF .
---- VAR Z2 -INF -2666.6667 +INF .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
---- 153 PARAMETER Z_INF = -2666.667 cota inferior
PARAMETER Z_SUP = -2666.667 cota superior
Equation Listing SOLVE COMPLETO Using LP From line 166
---- FOC =E= función objetivo problema completo
FOC.. X1(x1-1) + 2*X1(x1-2) + 2*X2(x2-1) + 3*X2(x2-2) + Z1 =E= 0 ; (LHS = 0)
R1(r1-1).. X1(x1-1) + X1(x1-2) =L= 600 ; (LHS = 0)
R1(r1-2).. X1(x1-1) - 2*X1(x1-2) =L= 0 ; (LHS = 0)
R2(r2-1).. X1(x1-1) + X1(x1-2) + X2(x2-1) + X2(x2-2) =L= 1000 ; (LHS = 1000)
R2(r2-2).. X1(x1-1) + X2(x2-1) =L= 500 ; (LHS = 333.333333333333)
04/11/2012 41
R2(r2-3).. - 2*X2(x2-1) + X2(x2-2) =L= 0 ; (LHS = 0)

Model Statistics SOLVE COMPLETO Using LP From line 166
MODEL STATISTICS
Solution Report SOLVE COMPLETO Using LP From line 166
MODEL COMPLETO OBJECTIVE Z1
Objective : -2666.666667
---- EQU FOC . . . 1.0000
r1-1 -INF . 600.0000 .
r1-2 -INF . . .
r2-1 -INF 1000.0000 1000.0000 -2.6667
r2-2 -INF 333.3333 500.0000 .
r2-3 -INF . . -0.3333
x1-1 . . +INF 1.6667
x1-2 . . +INF 0.6667
x2-1 . 333.3333 +INF .
x2-2 . 666.6667 +INF .
---- VAR Z1 -INF -2666.6667 +INF .
0 INFEASIBLE
0 UNBOUNDED
**** REPORT FILE SUMMARY
COPT C:\cplex.opt
**** FILE SUMMARY
Input C:\BD.gms
Output C:\BD.lst
Los lenguajes de modelado algebraico como GAMS permiten el uso de

diferentes optimizadores para resolver los problemas. En el caso anterior, el
optimizador escogido ha sido CPLEX. Para que las variables duales devueltas
por el subproblema sean correctas para generar un corte de Benders tanto en el
caso de optimalidad como en el de infactibilidad, para este optimizador deben
ser desactivadas las opciones de preproceso (preind 0) y de escalado (scaind
42 04/11/2012
-1) y el subproblema debe ser resuelto mediante el algoritmo simplex primal

(lpmethod 1). Con el uso de otros optimizadores se debe comprobar la
corrección de las variables duales tanto para optimalidad como para
infactibilidad. Si las variables duales no son las adecuadas para generar un corte
de infactibilidad, el subproblema de minimización de infactibilidades debe de ser
planteado explícitamente en el código GAMS y resuelto.
I.4.6 Algoritmo escrito en MATLAB

Para mostrar la elegancia, potencia y sencillez del lenguaje MATLAB [Matlab,
OptTool] para el estudio de este tipo de algoritmos se presenta a continuación el
problema lineal bietapa y el método de descomposición de Bd. No obstante, no
es recomendable el uso de este lenguaje para la escritura de grandes problemas
de optimización.
%
% programa lineal bietapa PL-2
%
dt2
[m1,n1]=size(A1) ;
[m2,n2]=size(A2) ;
A = [[A1, zeros(m1,n2)];[B1,A2]] ;
b = [b1; b2] ;
c = [c1; c2]' ;
[x,lambda] = lp(c,A,b,zeros(n1+n2,1),[],[0 0 0 0]) ;
x1 = x(1:n1)
x2 = x(n1+1:n1+n2)
pi21 = lambda(1:m1)
pi22 = lambda(m1+1:m1+m2)
%
% descomposición de Benders
%
[m1,n1]=size(A1) ;
[m2,n2]=size(A2) ;
l = 0 ;
A = [A1, zeros(m1,1)] ;
b = b1 ;
zsup = Inf ;
zinf = -Inf ;
epsilon = 0.0001 ;
while 1-abs(zinf/zsup) > epsilon
l = l + 1 ;
if l == 1,
x1 = x11 ;
zinf = c1' * x1 ;
else
A = [A; [-pi2(:,l-1)'*B1, -1]] ;
b = [b; -pi2(:,l-1)'*b2] ;
c = [c1; 1]' ;
[x1] = lp(c,A,b,zeros(n1,1)) ;
zinf = c * x1 ;
x1(n1+1)=[] ;
end
[x2,pi2t] = lp(c2',A2,b2-B1*x1,zeros(n2,1)) ;
pi2(:,l) = -pi2t(1:m2) ;
zsup = c1' * x1 + c2' * x2 ;
end
04/11/2012 43
I.4.7 Problema de transporte con coste fijo

Veamos a continuación un caso ejemplo característico de la aplicación del
método de descomposición de Bd. Se trata del problema de transporte donde
algunos o todos los arcos tienen un coste fijo asociado a la decisión de su
instalación o a su uso. El problema consiste en la minimización de los costes
fijos y variables sujeto a las restricciones de respetar las ofertas máximas de los
orígenes y las demandas en los destinos. El problema se formula de la siguiente
manera
min ∑ (cij x ij + fij yij )

x ij ,yij
ij
∑x j
ij ≤ ai ∀i
∑x i
ij ≥ bj ∀j
x ij ≤ M ij yij ∀ij
x ij ≥ 0, yij ∈ {0,1}
siendo cij el coste variable unitario de transporte, fij el coste fijo asociado a la
decisión de inversión en el arco ij , ai la oferta máxima de producto en el origen
i , bj la demanda del destino j , x ij la variable que indica el flujo que recorre el
arco ij , yij la variable que representa la decisión de inversión en el arco ij y
M ij una cota superior de cualquier flujo en dicho arco ij (por ejemplo,
M ij = min {ai , bj } ).
Las variables yij son binarias. Una vez conocidas el problema anterior es un
problema clásico de transporte. Las variables yij son las variables que complican
la resolución y, por consiguiente, son asignadas al problema maestro en un
entorno de descomposición de Bd.
El subproblema se formula de la siguiente manera
min ∑ cij x ij
x ij
ij
∑xj
ij ≤ ai ∀i
∑x
i
ij ≥ bj ∀j
x ij ≤ M ij yijk ∀ij : πijk

x ij ≥ 0
y el problema maestro como
44 04/11/2012
min θ + ∑ fij yij

yij ,θ
ij
δθ + ∑ πijl M ij yij ≥ f l + ∑ πijl M ij yijl l = 1, …, k

ij ij
yij ∈ {0,1}
A continuación se expresa en GAMS este problema para un caso ejemplo en
el que se suponen cuatro orígenes del producto y tres destinos de demanda. El
problema debe decidir la combinación óptima de arcos de entre todos los
posibles, dados en la figura 4.2.
Figura 4.2 Arcos posibles.
$title Fixed-charge transportation problem (FCTP) solved by Benders' decomposition

* relative optimality tolerance in solving MIP problems
option OptcR = 0
sets
L iterations / l1 * l20 /
LL(l) iterations subset
I origins / i1 * i4 /
J destinations / j1 * j3 /
* Begin problem data
parameters
A(i) product offer
/ i1 10, i2 30, i3 40, i4 20 /
B(j) product demand
/ j1 20, j2 50, j3 30 /
table C(i,j) per unit variable transportation cost
j1 j2 j3
i1 1 2 3
i2 3 2 1
i3 2 3 4
i4 4 3 2
table F(i,j) fixed transportation cost

j1 j2 j3
i1 10 20 30
i2 20 30 40
i3 30 40 50
i4 40 50 60
* End problem data
abort $(sum
sum[i, A(i)] < sum
sum sum[j, B(j)]) 'Infeasible problem'
parameters
BdTol relative Benders' tolerance / 1e-6 /
Z_Lower lower bound / -inf /
Z_Upper upper bound / inf /
Y_L (l,i,j) first stage variables values in iteration l
PI_L (l,i,j) dual variables of second stage constraints in iteration l
04/11/2012 45
Delta(l) cut type (feasibility 0 optimality 1) in iteration l

Z2_L (l) subproblem objective function value in iteration l
positive variable
X(i,j) arc flow
binary variable
Y(i,j) arc investment decision
variables
Z1 first stage objective function
Z2 second stage objective function
Theta recourse function
equations
EQ_Z1 first stage objective function
EQ_Z2 second stage objective function
EQ_OBJ complete problem objective function
Offer (i ) offer at origin
Demand ( j) demand at destination
FlowLimit(i,j) arc flow limit
Bd_Cuts (l) Benders' cuts ;
EQ_Z1 .. Z1 =e= sum

sum[(i,j), F(i,j)*Y(i,j)] + Theta ;
EQ_Z2 .. Z2 =e= sum[(i,j), C(i,j)*X(i,j)] ;
sum
EQ_OBJ .. Z1 =e= sum
sum[(i,j), F(i,j)*Y(i,j)] + sum
sum[(i,j), C(i,j)*X(i,j)] ;
Offer (i ) .. sum
sum[j, X(i,j)] =l= A(i) ;
Demand ( j) .. sum
sum[i, X(i,j)] =g= B(j) ;
FlowLimit(i,j) .. X(i,j) =l= min[A(i),B(j)] * Y(i,j) ;
Bd_Cuts(ll) .. Delta(ll) * Theta =g= Z2_L(ll) -
sum[(i,j), PI_L(ll,i,j) * MIN[A(i),B(j)] * (Y_L(ll,i,j) - Y(i,j))] ;
sum
model Master_Bd / EQ_Z1 , Bd_Cuts /
model Subproblem_Bd / EQ_Z2 , Offer, Demand, FlowLimit /
model Complete / EQ_OBJ, Offer, Demand, FlowLimit / ;
X.up(i,j) = min[A(i),B(j)]
* to allow CPLEX correctly detect rays in an infeasible problem
* only simplex method can be used and no preprocessing neither scaling options
* optimality and feasibility tolerances are very small to avoid primal generation
file COPT / cplex.opt /
put COPT putclose 'ScaInd -1' / 'LPMethod 1' / 'PreInd 0' / 'EpOpt 1e-9' / 'EpRHS 1e-9'
/ ;
Subproblem_Bd.OptFile = 1 ;
* parameters initialization
LL (l) = no ;
Theta.fx = 0 ;
Delta (l) = 0 ;
Z2_L (l) = 0 ;
PI_L(l,i,j) = 0 ;
Y_L (l,i,j) = 0 ;
* Benders' algorithm iterations
loop (l $(abs(1-Z_Lower/Z_Upper) > BdTol),
* solving master problem
solve Master_Bd using MIP minimizing Z1 ;
* storing the master solution
Y_L(l,i,j) = Y.l(i,j) ;
* fixing first-stage variables and solving subproblem
Y.fx (i,j) = Y.l(i,j) ;
* solving subproblem
solve Subproblem_Bd using RMIP minimizing Z2 ;
* storing parameters to build a new Benders' cut
if (Subproblem_Bd.ModelStat = 4,
46 04/11/2012
Delta(l) = 0 ;
Z2_L (l) = Subproblem_Bd.SumInfes ;
else
* updating lower and upper bound
Z_Lower = Z1.l ;
Z_Upper = min(Z_Upper, Z1.l - Theta.l + Z2.l) ;
Delta(l) = 1 ;
Z2_L (l) = Z2.l ;
Theta.lo = -inf
inf ;
Theta.up = inf ;
) ;
PI_L(l,i,j) = FlowLimit.m(i,j) ;
Y.lo( i,j) = 0 ;
Y.up( i,j) = 1 ;
* increase the set of Benders' cuts

LL(l) = yes ;
) ;
solve Complete using MIP minimizing Z1
La solución óptima es y11 = y23 = y 31 = y 32 = y 42 = 1 y se alcanza en 12

teraciones9 con un coste total, fijo más variable, de 380. La evolución de las
combinaciones de arcos propuestas por el algoritmo de descomposición se
presenta en la figura 4.3. Debe destacarse que solamente las propuestas de las
iteraciones 7, 10, 11 y 12 son factibles para el subproblema, de modo que la
mayoría de cortes que el algoritmo genera son de infactibilidad y sólo dos cortes
son de optimalidad. Por último, la tabla 4.1 presenta la evolución de la cota
superior y la cota inferior del problema durante el algoritmo. En este ejemplo, el
algoritmo converge cuando la cota inferior es exactamente la misma que la cota
superior. La convergencia se alcanza también cuando aparece una propuesta
repetida para evitar entrar en un ciclo de soluciones. En este ejemplo, la
propuesta de arcos de la última iteración es, de hecho, la propuesta de la
iteración 12.
9
En este caso particular el subproblema tiene soluciones alternativas por lo que el algoritmo
puede converger en diferente número de iteraciones con otro optimizador u otra versión del
mismo optimizador.
04/11/2012 47
I I I I I
I I I
Figura 4.3 Evolución de arcos propuestos por el método de descomposición de Benders.
Iteration Lower Bound Upper Bound

1a6 −∞ ∞
7 140 390
8 140 390
9 140 390
10 360 390
11 370 390
12 380 380
Tabla 4.1 Evolución de cotas en el método de descomposición de Benders.
450
400
350
Objective function
300
250
200 Lower Bound
150 Upper Bound
100
50
0
7 8 9 10 11 12
Iterations
48 04/11/2012
Figura 4.4 Evolución de las cotas inferior y superior.
I.4.8 Variables duales del problema lineal bietapa

Al resolver el problema original bietapa mediante el método de
descomposición de Bd no se pueden calcular los valores de las variables duales
del problema completo con respecto a las ecuaciones de la segunda etapa. Esta
modificación del algoritmo [Entriken:87] permite obtenerlos. El problema
maestro para la iteración j se formula:
min c1T x 1 + θ2
x1 ,y1 ,θ2
A1x 1 = b1 : π1
B1x 1 − Iy1 = 0 : ρ1 (4.25)
π2lT y1 + θ2 ≥ π2lTb2 l = 1, …, j
x1 ≥ 0
donde θ2 ∈ ℝ , y1 ∈ ℝ m2 y π1 y ρ1 son las variables duales de las respectivas

restricciones.
El subproblema para cada iteración se formula como:
min cT2 x 2
x 2 ,z 2
A2x 2 + y1z 2 = b2 : π2
(4.26)
z2 = 1 : µ2
x2 ≥ 0
Con esta modificación, las variables duales del problema lineal bietapa PL-2
se pueden obtener a partir de las variables duales del maestro.
πˆ1*  π1* 
*    
πˆ =  *  =  *  (4.27)
πˆ2   ρ1 
   
donde π̂1* y π̂2* son las variables duales del primer y segundo conjunto de
restricciones respectivamente del problema lineal PL-2 para la solución óptima y
π1* y ρ1* son las variables duales calculadas en el maestro para la solución
óptima.
Como se observa tanto el problema maestro como el subproblema han
aumentado de tamaño. El problema maestro aumenta en m2 el número de
variables y de restricciones. El subproblema aumenta en una variable y una
restricción. Si el aumento de tamaño en el problema maestro no es razonable
una alternativa es resolver el problema lineal bietapa por el procedimiento
04/11/2012 49
original y una vez convergido resolver sólo el problema maestro con esta
modificación.
50 04/11/2012
I.5 Descomposición de Dantzig-Wolfe

El método de Dantzig-Wolfe (DW) [Dantzig:63, Dantzig:60, Geoffrion:70,
Geoffrion:71] también se denomina descomposición dual o descomposición por
precios porque el maestro envía precios o variables duales a los subproblemas,
generación de columnas (column generation) porque el maestro incrementa el
número de variables en cada iteración. Se utiliza cuando un conjunto de
restricciones complican la solución del problema10. Es decir, ésta resulta mucho
más sencilla cuando aquéllas no están presentes. El método de DW realiza una
relajación de estas restricciones. Implícitamente se supone que el número de
restricciones de este tipo es reducido frente al número total de restricciones
m1 ≪ m2 .
Un caso típico es la planificación de una compañía con múltiples divisiones,
que implica la coordinación de las decisiones que las divisiones realizan por
separado. Cada división opera con una autonomía considerable y optimiza
únicamente sus propias operaciones. Sin embargo, es necesaria una coordinación
corporativa para asignar recursos centrales entre las divisiones y conseguir una
operación óptima para toda la compañía. Otro caso en sistemas de energía
eléctrica es la programación semanal de grupos de generación donde las
restricciones que acoplan las decisiones individuales de los grupos son las
restricciones de cobertura de demanda en cada hora. Si esas restricciones no
existieran los grupos podrían hacer su programa semanal individualmente.
Descompone el problema lineal en un problema maestro y un subproblema.
El problema maestro representa una combinación lineal de las soluciones y las
restricciones que complican. El subproblema genera nuevas soluciones. El
algoritmo es iterativo y alterna entre la solución del maestro y del subproblema.
Supongamos que se desea resolver el siguiente problema de programación
lineal
min c1T x 1
x1
A1x 1 = b1
(5.1)
A2x 1 = b2
x1 ≥ 0
10
Aquí se identifican con las de la primera etapa aunque como tal el problema no es bietapa
por carecer de variables de la segunda etapa.
04/11/2012 51
5 DESCOMPOSICIÓN DE DANTZIG-WOLFE
donde A1x 1 = b1 son las restricciones que complican (en número reducido) y
A2x 1 = b2 son las restricciones con estructura diagonal por bloques y separable
(en número mucho mayor).
El problema lineal anterior se puede formular también como:
min cT1 x 1
x1
A1x 1 = b1 (5.2)
x1 ∈ K
siendo K la región convexa definida como:
K = {x 1 | A2x 1 = b2 , x 1 ≥ 0} (5.3)
Pero todo punto de un poliedro convexo (no vacío y acotado) se puede poner
como una combinación lineal convexa (es decir, una linealización interior) de
sus vértices x 1l , l = 1, …, ν .
 ν ν

K = ∑ x 1l λl | ∑ λl = 1, λl ≥ 0 (5.4)
 l =1 l =1

Entonces, el problema maestro completo se formula como:

ν
min ∑ (c1T x 1l )λl
λl
l =1
ν
l
∑ (A x )λ
l =1
1 1 l = b1 : π2
(5.5)
ν
∑λ
l =1
l =1 :µ
λl ≥ 0 l = 1, …, ν
Obsérvese que las variables de este problema son los pesos λl de los vértices
en lugar de las variables originales del problema x 1 . El número total de vértices
 n1 
de un politopo es m  siendo A2 ∈ ℝ m2×n1 . Algunos de ellos son factibles y otros
 2 
infactibles.
En lugar de enumerar todos los vértices, el algoritmo de descomposición
resuelve iterativamente el proceso, introduciendo un vértice nuevo en cada
iteración a medida que se necesitan. La región K viene definida por la
combinación lineal convexa del conjunto de vértices introducido hasta ese
momento. Luego, la resolución del maestro se puede interpretar como la del
problema completo original pero sobre una región factible K cada vez mayor,
correspondiente al segundo conjunto de restricciones (que no complican). Su
función objetivo disminuye en cada iteración al introducir una nueva variable y,
52 04/11/2012
por consiguiente, será monótona decreciente. Es una cota superior de la función

objetivo del problema completo.
El problema maestro restringido para la iteración j es:
Maestro Dantzig-Wolfe
j
λl
l =1
j
l
∑ (A x )λ
l =1
1 1 l = b1 : π2
(5.6)
j
∑λ
l =1
l =1 :µ
λl ≥ 0 l = 1, …, j
Supongamos que se dispone de una solución inicial factible x 11 ∈ K y el
problema maestro se resuelve por el método simplex [Dantzig:63, Dantzig:91d].
La condición de optimalidad de un problema lineal de minimización determina
que los costes reducidos deben ser no negativos:
 *
A1x 1 
T
c x − π
1
*
1 ( T
2 )
µ  ≥0
 1 
(5.7)
equivalente a:
θ2 = min (c1T − πT2 A1 ) x 1 − µ (5.8)

x1 ∈K
o expresando la minimización como problema lineal, el subproblema para la

iteración j es:
Subproblema Dantzig-Wolfe
θ2 = min (c1T − π2jT A1 ) x 1 − µ j

x1
A2x 1 = b2 (5.9)
x1 ≥ 0
Luego, el subproblema obtiene en cada iteración el vértice con el menor coste

reducido y éste se incorpora al maestro. Cuando el valor de la función objetivo
del subproblema es positivo o nulo se ha alcanzado el óptimo. Obsérvese que el
subproblema resulta separable si las restricciones A2x 1 = b2 tienen una
estructura diagonal por bloques.
04/11/2012 53
Veamos otra forma de deducir el algoritmo de descomposición de DW.

Formamos el lagrangiano del problema (5.2) relajando las restricciones de
complicación
L(x 1, π2′ ) = c1T x 1 + π2′T (A1x 1 − b1 ) (5.10)

siendo π2′ los multiplicadores de Lagrange o variable duales de las restricciones
que complican A1x 1 = b1 .
La función dual θ2 (π2′ ) tiene esta expresión
θ2 (π2′ ) = min L(x 1, π2′ ) = c1T x 1 + π2′T (A1x 1 − b1 )
x1 ∈K
(5.11)
= −π2′Tb1 + min(c1T + π2′T A1 )x 1
x1 ∈K
La función dual siempre es cóncava independientemente de que el problema

sea lineal o lineal entero mixto. El subproblema o función dual se puede
interpretar como el problema completo original pero en el que se han relajado
las restricciones de complicación. Por tanto, el algoritmo convergerá cuando la
función objetivo del subproblema, que es una cota inferior de la función objetivo
del problema, sea igual a la del maestro, que es una cota superior. Su cálculo se
puede expresar también como
Subproblema Dantzig-Wolfe
θ2 (π2 ) = min(c1T − πT2 A1 )x 1 + πT2 b1

x1
A2x 1 = b2 (5.12)
x1 ≥ 0
donde π2 = −π2′ .
Obsérvese que el subproblema formulado en (5.9) es diferente del formulado
en (5.12), aunque son equivalentes para el algoritmo de descomposición el valor
de la función objetivo cambia. En el primer caso, el problema (5.9) evalúa los
costes reducidos y, por consiguiente, tomará valor 0 al alcanzar el óptimo. En el
segundo caso, el problema (5.12) evalúa el lagrangiano y, por lo tanto, tomará el
valor óptimo del problema completo al alcanzar el óptimo el algoritmo de
descomposición.
Como el óptimo de un problema lineal se alcanza en un vértice podemos
tomar el mínimo de la función dual para todos los vértices x 1l , l = 1, …, ν
θ2 (π2 ) = πT2 b1 + min(

l
c1T − πT2 A1 )x 1l (5.13)
x1
que también se puede expresar como
54 04/11/2012
θ2 (π2 ) ≤ πT2 b1 + (cT1 − πT2 A1 )x 11

θ2 (π2 ) ≤ πT2 b1 + (cT1 − πT2 A1 )x 12
(5.14)
⋮
θ2 (π2 ) ≤ πT2 b1 + (cT1 − πT2 A1 )x 1ν
o expresado en forma de problema lineal y reagrupando términos
max θ2
θ2 ,π2
θ2 + (A1x 11 − b1 )T π2 ≤ c1T x 11 : λ1
θ2 + (A1x 12 − b1 )T π2 ≤ c1T x 12 : λ2 (5.15)
⋮
θ2 + (A1x 1ν − b1 )T π2 ≤ c1T x 1ν : λν
Si planteamos el problema dual de éste tenemos
ν
λl
l =1
ν
∑λ
l =1
l =1 : µ = θ2
(5.16)
ν
l
∑ (A x
l =1
1 1 − b1 )λl = 0 : π2
λl ≥ 0 l = 1, …, ν
Manipulando el segundo bloque de restricciones
ν ν ν ν
∑ (A1x1l − b1 )λl = ∑ (A1x1l )λl − ∑ b1λl = ∑ (A1x 1l )λl − b1 = 0

l =1 l =1 l =1 l =1
el maestro completo de DW se formula como en la ecuación (5.5).

Este método tiene una interesante interpretación económica. La función
objetivo del subproblema es su propio coste menos el valor π2 que tiene el
cumplimiento de las restricciones del maestro. Por eso este algoritmo se
denomina también descomposición por precios. Consigue indirectamente resolver
el problema de asignación de recursos escasos mediante la valoración por precios
de dichos recursos. Este valor se ajusta en cada iteración para cumplir las
restricciones del maestro. La solución del subproblema es la nueva propuesta
que se hace al maestro. El maestro utiliza todas las propuestas hechas hasta el
momento y las combina para conseguir una nueva solución óptima, de la que
obtiene los valores (variables duales) de las restricciones.
04/11/2012 55
I.5.2 Algoritmo
El tamaño del problema lineal completo es (m1 + m2 ) × n1 . El tamaño del
problema maestro es (m1 + 1) × j y el del subproblema m2 × n1 . El problema
maestro aumenta en una variable en cada iteración lo que ocasiona una pérdida
de optimalidad, luego conviene que sea resuelto por el método simplex primal.
El subproblema se ve afectado únicamente en su función objetivo luego también
se recomienda utilizar el método simplex primal.
En una iteración del algoritmo se resuelve el subproblema11. Mientras su
función objetivo (i.e., la condición de optimalidad del problema maestro cuando
es formulado según la ecuación (5.9)) sea negativa la solución x 1j obtenida en el
subproblema es enviada al maestro. Cuando es positiva o nula se alcanza el
óptimo y el algoritmo acaba. A continuación se resuelve el problema maestro y
se pasa el valor del vector de variables duales π2j y µ j al subproblema. En la
primera iteración se suponen unos valores iniciales razonables (o nulos) de la
variables duales π21 y µ1 y se resuelve el subproblema.
En general, el algoritmo progresa rápidamente al comienzo pero luego
converge muy lentamente. Por esta razón, a veces se termina prematuramente
con una solución subóptima. Un criterio de parada para detener previamente el
algoritmo, cuando el subproblema es formulado según la ecuación (5.12), puede
ser
z ≤ z 1* ≤ z (5.17)
siendo z y z los valores de las funciones objetivo de maestro y subproblema,
respectivamente, en una iteración cualquiera. Luego cuando la diferencia entre
ambas está por debajo de una cierta tolerancia se detiene el algoritmo. Esta
diferencia es teóricamente 0 en el óptimo cuando el problema completo es lineal.
Esquemáticamente el algoritmo se representa en la figura 4.4.
SUBPROBLEMA
| ↑
j j
precios sombra π , µ
2 propuesta x 1j
↓ |
PROBLEMA MAESTRO RESTRINGIDO
Figura 4.4 Método de descomposición de Dantzig-Wolfe.
11
En el método de descomposición de Dantzig-Wolfe la iteración empieza de forma natural
por el subproblema mientras que en Benders empezaba por el maestro.
56 04/11/2012
La solución óptima x 1* y el valor de la función objetivo z 1* , una vez

convergido el algoritmo, viene dada por la combinación lineal de las soluciones
de todas las iteraciones
j
x 1* = ∑ x 1l λl (5.18)
l =1
j
z 1* = ∑ (cT1 x 1l )λl (5.19)
l =1
El método puede ser aplicado a problemas de optimización no lineal pero se

requiere que el conjunto de restricciones que no complican definan un poliedro
convexo para que los puntos del interior se puedan definir como combinación
lineal convexa de sus vértices. Además para que tenga sentido utilizar las
variables duales de las restricciones que complican el problema maestro también
debe ser convexo.
En general en el método de DW, puede resultar difícil conseguir factibilidad
en las primeras iteraciones en las restricciones del maestro por lo que desde un
punto de vista de implantación es conveniente la introducción de variables de
holgura (en restricciones de tipo ≥ ) o exceso (en restricciones de tipo ≤ )
penalizadas en la función objetivo del maestro.

Veamos gráficamente con un ejemplo el funcionamiento del algoritmo. Se
desea resolver el problema de optimización siguiente
min− 4x − y − 6z
3x + 2y + 4z = 17
1≤x ≤2
1≤y ≤2
1≤z ≤2
La primera restricción se considera de complicación, luego irá al problema
maestro y las tres últimas formarán parte del subproblema. Empezamos
resolviendo el subproblema sin variables duales
min− 4x − y − 6z
1≤x ≤2
1≤y ≤2
1≤z ≤2
04/11/2012 57
El óptimo resulta ser (x , y, z ) = (2,2, 2) y se pasa este punto y el valor de la

función objetivo, que es la cota inferior z = −22 , al problema maestro
2
 
 
min −4 −1 −6 2 λ1
( )  
 
2
2
 
 
(3 2 4 2 λ1 = 17
)
 
2
 
λ1 = 1
λ1 ≥ 0
Este problema resulta infactible. Una posible solución es introducir unas
variables de exceso y holgura y penalizarlas con un valor elevado en la función
objetivo o bien dar un valor elevado a la variable dual de la restricción
π = −100 . La cota superior ahora es z = ∞ .
El subproblema se formula ahora como
min (−4 + 3 ⋅ 100) x + (−1 + 2 ⋅ 100) y + (−6 + 4 ⋅ 100) z − 17 ⋅ 100
1≤x ≤2
1≤y ≤2
1≤z ≤2
El óptimo es el punto (x , y, z ) = (1,1,1) , el valor de la función objetivo
z = −811 y un valor para el maestro de −11 . El nuevo problema maestro es
min− 22λ1 − 11λ2
18λ1 + 9λ2 = 17
λ1 + λ2 = 1
λ1, λ2 ≥ 0
de donde resulta π = −1.222 y cota superior z = −20.777 lo que da lugar al
subproblema que se formula ahora como
min− 0.333x + 1.444y − 1.111z − 20.777
1≤x ≤2
1≤y ≤2
1≤z ≤2
58 04/11/2012
El nuevo punto óptimo es (x , y, z ) = (2,1, 2) y el valor de la función objetivo

z = −22.222 y con un valor para el maestro de −21 . El nuevo problema
maestro es
min− 22λ1 − 11λ2 − 21λ3
18λ1 + 9λ2 + 16λ3 = 17
λ1 + λ2 + λ3 = 1
λ1, λ2 , λ3 ≥ 0
de donde resulta π = −0.5 , (λ1, λ2 , λ3 ) = (0.5, 0, 0.5) y cota superior z = −21.5 .
El subproblema es ahora
min− 2.5x − 4z − 8.5
1≤x ≤2
1≤y ≤2
1≤z ≤2
z = −21.5 . Al ser iguales la cota superior e inferior el algoritmo converge.
La solución óptima es la combinación lineal de las soluciones anteriores
    1 2  2 
x  2      
y  = 2 λ + 1 λ + 1 λ = 1.5
 
    1   2   3  
       
z  2 1 2  2 
Los pasos que sigue el algoritmo para resolver el siguiente problema lineal se
pueden ver en la implantación hecha en GAMS.
min− x − y − 2u − v
x +2y +2u +v ≤ 10
x +3y ≤ 40
2x +y ≤ 30
u ≤ 10
v ≤ 10
u +v ≤ 15
x, y, u, v ≥0
04/11/2012 59
40 1 3 
−1    
  30 2 1 
−1      
   10  
−2    (
entonces, c1 =   , b1 = (10) , b2 =   , A1 = 1 2 2 1 y A2 = ) 
1 

  10  1
  
−1  
15 

1 1
Obsérvese que en este problema, cuando se resuelve por descomposición, se

introduce una variable de exceso en la primera restricción que se penaliza en la
función objetivo para evitar que el maestro resulte infactible. Se necesitan tres
iteraciones para que converja por el método de descomposición de DW. El
algoritmo de DW escrito en GAMS se presenta a continuación.
$title Descomposición de Danztig-Wolfe (DW)

sets
L máximo numero de iteraciones / iter-1 * iter-5 /
N1 variables
scalar
Z_INF cota inferior de la solución / -inf /
Z_SUP cota superior de la solución / inf /
TOL tolerancia en optimalidad / 1e-6 /
PNLZ penalización restricciones primera etapa / 1000 /
parameters
X1_L(l,n1) soluciones del subproblema
sets
M1 número restricciones primera etapa / rdos-1 /
M2 número restricciones segunda etapa / runo-1 * runo-5 /
N1 variables / xuno-1 * xuno-4 /
parameters
C1(n1) coeficientes función objetivo primera etapa
/ xuno-1 -1
xuno-2 -1
xuno-3 -2
xuno-4 -1 /
/ rdos-1 10 /
/ runo-1 40
runo-2 30
runo-3 10
runo-4 10
runo-5 15 /
xuno-1 xuno-2 xuno-3 xuno-4
rdos-1 1 2 2 1
xuno-1 xuno-2 xuno-3 xuno-4
runo-1 1 3
runo-2 2 1
runo-3 1
runo-4 1
runo-5 1 1
positive variables
60 04/11/2012
LAMBDA(l) pesos de las soluciones

EXC(m1) variables de exceso en restricciones primera etapa <=
variables
Z2 función objetivo subproblema y completo
Z1 función objetivo maestro
equations
FOM función objetivo maestro
FOS función objetivo subproblema
R1C(m1) restricciones primera etapa
COMBLIN combinación lineal de los pesos de las soluciones ;
FOM .. Z1 =E= sum
sum((j,n1), C1(n1)*X1_L(j,n1)*LAMBDA(j)) + PNLZ * sum
sum(m1, EXC(m1)) ;
FOS .. Z2 =E= sum
sum(n1, (C1(n1) - sum
sum(m1, R1.M(m1)*A1(m1,n1)))*X1(n1)) - COMBLIN.M ;
FOC .. Z2 =E= sum

sum(n1, C1(n1)*X1(n1)) ;
R1C(m1) .. sum(n1,
sum A1(m1,n1)*X1(n1)) =L= B1(m1) ;
* se introducen variables de exceso en restricciones primera etapa <=
* para evitar infactibilidad
R1(m1) .. sum((j,n1),
sum A1(m1,n1)*X1_L(j,n1)*LAMBDA(j)) - EXC(m1) =L= B1(m1) ;
R2(m2) .. sum(n1,
sum A2(m2,n1)*X1(n1)) =L= B2(m2) ;
COMBLIN .. sum(j,
sum LAMBDA(j)) =E= 1 ;
model MAESTRO / FOM, R1, COMBLIN /
model SUB / FOS, R2 /
model COMPLETO / FOC, R1C, R2 /
file
put COPT putclose 'preind 0' / 'epopt 1e-9' / 'eprhs 1e-9'
SUB.OptFile = 1 ;
MAESTRO.OptFile = 1 ;
J(l) = NO ;
X1_L(l,n1) = 0 ;
LAMBDA.UP(l) = 1 ;
* inicialización de variables
Z2.L = - 10 * TOL ;
R1.M(m1) = 0 ;
COMBLIN.M = 0 ;
loop(l
loop $(abs(Z2.L) > TOL),
if (ord
ord(l)
ord > 1,
solve MAESTRO USING LP MINIMIZING Z1 ;
Z_SUP = Z1.L ;
) ;
solve SUB USING LP MINIMIZING Z2 ;

X1_L(l,n1) = X1.L(n1) ;
Z_INF = Z2.L + COMBLIN.M + sum
sum(m1, R1.M(m1)*B1(m1)) ;
display Z_INF, Z_SUP ;
J(l) = yes ;
) ;
abort $(abs(Z2.L) > TOL) 'Máximo numero de iteraciones alcanzado' ;
* resultados de la descomposición
X1.L(n1) = sum(j,
sum X1_L(j,n1)*LAMBDA.L(j)) ;
display X1.L ;
04/11/2012 61

Descomposición de Danztig-Wolfe (DW)
---- FOS =E= función objetivo subproblema
FOS.. X1(xuno-1) + X1(xuno-2) + 2*X1(xuno-3) + X1(xuno-4) + Z2 =E= 0 ; (LHS = -1E-5, INFES = 1E-5 ***)
R2(runo-1).. X1(xuno-1) + 3*X1(xuno-2) =L= 40 ; (LHS = 0)
R2(runo-2).. 2*X1(xuno-1) + X1(xuno-2) =L= 30 ; (LHS = 0)
R2(runo-3).. X1(xuno-3) =L= 10 ; (LHS = 0)
R2(runo-4).. X1(xuno-4) =L= 10 ; (LHS = 0)
R2(runo-5).. X1(xuno-3) + X1(xuno-4) =L= 15 ; (LHS = 0)
LOOPS L iter-1
MODEL STATISTICS
L O O P S L iter-1
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
Objective : -45.000000
---- EQU FOS . . . 1.0000
runo-1 -INF 40.0000 40.0000 -0.2000
runo-2 -INF 30.0000 30.0000 -0.4000
runo-3 -INF 10.0000 10.0000 -1.0000
runo-4 -INF 5.0000 10.0000 .
runo-5 -INF 15.0000 15.0000 -1.0000
xuno-1 . 10.0000 +INF .
xuno-2 . 10.0000 +INF .
xuno-3 . 10.0000 +INF .
xuno-4 . 5.0000 +INF .
---- VAR Z2 -INF -45.0000 +INF .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
---- 122 PARAMETER Z_INF = -45.000 cota inferior de la solución
PARAMETER Z_SUP = +INF cota superior de la solución
62 04/11/2012

---- FOM =E= función objetivo maestro
FOM.. 45*LAMBDA(iter-1) - 1000*EXC(rdos-1) + Z1 =E= 0 ; (LHS = 0)
R1(rdos-1).. 55*LAMBDA(iter-1) - EXC(rdos-1) =L= 10 ; (LHS = 0)
---- COMBLIN =E= combinación lineal de los pesos de las soluciones
COMBLIN.. LAMBDA(iter-1) =E= 1 ; (LHS = 0, INFES = 1 ***)
LOOPS L iter-2
MODEL STATISTICS
NON ZERO ELEMENTS 6
L O O P S L iter-2
**** OBJECTIVE VALUE 44955.0000
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
Objective : 44955.000000
---- EQU FOM . . . 1.0000
rdos-1 -INF 10.0000 10.0000 -1000.0000
---- EQU COMBLIN 1.0000 1.0000 1.0000 54955.0000
COMBLIN combinación lineal de los pesos de las soluciones
---- VAR LAMBDA pesos de las soluciones
iter-1 . 1.0000 1.0000 .
---- VAR EXC variables de exceso en restricciones primera etapa <=
rdos-1 . 45.0000 +INF .
---- VAR Z1 -INF 44955.0000 +INF .
0 INFEASIBLE
0 UNBOUNDED
FOS.. - 999*X1(xuno-1) - 1999*X1(xuno-2) - 1998*X1(xuno-3) - 999*X1(xuno-4) + Z2 =E= -54955 ;
(LHS = -55000, INFES = 45 ***)
04/11/2012 63

R2(runo-3).. X1(xuno-3) =L= 10 ; (LHS = 10)
R2(runo-4).. X1(xuno-4) =L= 10 ; (LHS = 5)
LOOPS L iter-2
MODEL STATISTICS
L O O P S L iter-2
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
Objective : -54955.000000
---- EQU FOS -54955.0000 -54955.0000 -54955.0000 1.0000
runo-1 -INF . 40.0000 .
runo-2 -INF . 30.0000 .
runo-3 -INF . 10.0000 .
runo-4 -INF . 10.0000 .
runo-5 -INF . 15.0000 .
xuno-1 . . +INF 999.0000
xuno-2 . . +INF 1999.0000
xuno-3 . . +INF 1998.0000
xuno-4 . . +INF 999.0000
---- VAR Z2 -INF -54955.0000 +INF .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
PARAMETER Z_SUP = 44955.000 cota superior de la solución
FOM.. 45*LAMBDA(iter-1) - 1000*EXC(rdos-1) + Z1 =E= 0 ; (LHS = 0)
R1(rdos-1).. 55*LAMBDA(iter-1) - EXC(rdos-1) =L= 10 ; (LHS = 10)
COMBLIN.. LAMBDA(iter-1) + LAMBDA(iter-2) =E= 1 ; (LHS = 1)
64 04/11/2012
LOOPS L iter-3
MODEL STATISTICS
NON ZERO ELEMENTS 7
L O O P S L iter-3
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOM . . . 1.0000
rdos-1 -INF 10.0000 10.0000 -0.8182
---- EQU COMBLIN 1.0000 1.0000 1.0000 EPS
iter-1 . 0.1818 1.0000 .
iter-2 . 0.8182 1.0000 .
rdos-1 . . +INF 999.1818
---- VAR Z1 -INF -8.1818 +INF .
0 INFEASIBLE
0 UNBOUNDED
FOS.. 0.181818181818182*X1(xuno-1) - 0.636363636363636*X1(xuno-2) + 0.363636363636364*X1(xuno-3)
+ 0.181818181818182*X1(xuno-4) + Z2 =E= EPS ; (LHS = -54955, INFES = 54955 ***)
R2(runo-3).. X1(xuno-3) =L= 10 ; (LHS = 0)
R2(runo-4).. X1(xuno-4) =L= 10 ; (LHS = 0)
LOOPS L iter-3
MODEL STATISTICS
04/11/2012 65

L O O P S L iter-3
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOS . . . 1.0000
runo-1 -INF 15.0000 40.0000 .
runo-2 -INF 30.0000 30.0000 -0.0909
runo-3 -INF 10.0000 10.0000 -0.1818
runo-4 -INF 5.0000 10.0000 .
runo-5 -INF 15.0000 15.0000 -0.1818
xuno-1 . 15.0000 +INF .
xuno-2 . . +INF 0.7273
xuno-3 . 10.0000 +INF .
xuno-4 . 5.0000 +INF .
---- VAR Z2 -INF -7.2727 +INF .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
PARAMETER Z_SUP = -8.182 cota superior de la solución
FOM.. 45*LAMBDA(iter-1) + 40*LAMBDA(iter-3) - 1000*EXC(rdos-1) + Z1 =E= 0 ; (LHS = 0)
R1(rdos-1).. 55*LAMBDA(iter-1) + 40*LAMBDA(iter-3) - EXC(rdos-1) =L= 10 ; (LHS = 10)
COMBLIN.. LAMBDA(iter-1) + LAMBDA(iter-2) + LAMBDA(iter-3) =E= 1 ; (LHS = 1)
LOOPS L iter-4
MODEL STATISTICS
L O O P S L iter-4
66 04/11/2012
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOM . . . 1.0000
rdos-1 -INF 10.0000 10.0000 -1.0000
---- EQU COMBLIN 1.0000 1.0000 1.0000 EPS
iter-1 . . 1.0000 10.0000
iter-2 . 0.7500 1.0000 .
iter-3 . 0.2500 1.0000 .
rdos-1 . . +INF 999.0000
---- VAR Z1 -INF -10.0000 +INF .
0 INFEASIBLE
0 UNBOUNDED
FOS.. - X1(xuno-2) + Z2 =E= EPS ; (LHS = -7.27272727272728, INFES = 7.27272727272728 ***)
R2(runo-3).. X1(xuno-3) =L= 10 ; (LHS = 10)
R2(runo-4).. X1(xuno-4) =L= 10 ; (LHS = 5)
LOOPS L iter-4
MODEL STATISTICS
L O O P S L iter-4
04/11/2012 67

scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
Objective : 0.000000
---- EQU FOS . . . 1.0000
runo-1 -INF . 40.0000 .
runo-2 -INF . 30.0000 .
runo-3 -INF . 10.0000 .
runo-4 -INF . 10.0000 .
runo-5 -INF . 15.0000 .
xuno-1 . . +INF EPS
xuno-2 . . +INF 1.0000
xuno-3 . . +INF EPS
xuno-4 . . +INF EPS
---- VAR Z2 -INF . +INF .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
PARAMETER Z_SUP = -10.000 cota superior de la solución
---- 132 VARIABLE X1.L variables primera etapa
xuno-1 3.750, xuno-3 2.500, xuno-4 1.250
FOC.. X1(xuno-1) + X1(xuno-2) + 2*X1(xuno-3) + X1(xuno-4) + Z2 =E= 0 ; (LHS = 10, INFES = 10 ***)
---- R1C =L= restricciones primera etapa
R1C(rdos-1).. X1(xuno-1) + 2*X1(xuno-2) + 2*X1(xuno-3) + X1(xuno-4) =L= 10 ; (LHS = 10)
R2(runo-1).. X1(xuno-1) + 3*X1(xuno-2) =L= 40 ; (LHS = 3.75)
R2(runo-2).. 2*X1(xuno-1) + X1(xuno-2) =L= 30 ; (LHS = 7.5)
R2(runo-3).. X1(xuno-3) =L= 10 ; (LHS = 2.5)
R2(runo-4).. X1(xuno-4) =L= 10 ; (LHS = 1.25)
R2(runo-5).. X1(xuno-3) + X1(xuno-4) =L= 15 ; (LHS = 3.75)
MODEL STATISTICS
68 04/11/2012

---- EQU FOC . . . 1.0000
---- EQU R1C restricciones primera etapa
rdos-1 -INF 10.0000 10.0000 -1.0000
runo-1 -INF . 40.0000 .
runo-2 -INF . 30.0000 .
runo-3 -INF . 10.0000 .
runo-4 -INF 10.0000 10.0000 EPS
runo-5 -INF 10.0000 15.0000 .
xuno-1 . . +INF EPS
xuno-2 . . +INF 1.0000
xuno-3 . . +INF .
xuno-4 . 10.0000 +INF .
---- VAR Z2 -INF -10.0000 +INF .
0 INFEASIBLE
0 UNBOUNDED
COPT C:\cplex.opt
**** FILE SUMMARY
Input C:\DW.gms
Output C:\DW.lst
I.5.4 Algoritmo escrito en MATLAB

El algoritmo escrito en MATLAB se presenta a continuación.
%
% descomposicion de Dantzig-Wolfe
%
[m1,n1]=size(A1) ;
[m2,n2]=size(A2) ;
l = 1 ;
optcnd = -1 ;
x1 = x11 ;
AA = B1*x1 ;
cc = c1'*x1 ;
while optcnd < 0
A = [[AA, A2]; [ones(1,l), zeros(1,n2)]] ;
b = [b2; 1] ;
c = [cc, c2'] ;
[x2,pi2t] = lp(c,A,b,zeros(l+n2,1)) ;
pi2(:,l) = -pi2t(1:m2) ;
[x1t] = lp([c1'-pi2(:,l)'*B1],A1,b1,zeros(n1,1)) ;
l = l + 1 ;
x1(:,l) = x1t ;
AA = [AA, B1*x1(:,l)] ;
cc = [cc, c1'*x1(:,l)] ;
optcnd = [c1'-pi2(:,l-1)'*B1]*x1(:,l) + pi2t(m2+1) ;
end
04/11/2012 69
I.6 Relajación lagrangiana

El método de relajación lagrangiana (RL) [Geoffrion:70, Geoffrion:74] es una
de las técnicas más extendidas en optimización discreta. Se utiliza en situación
análoga al método de Dantzig-Wolfe (DW), cuando existen restricciones que
complican la solución del problema. Esto es, la resolución del problema sin esas
ecuaciones tiene una estructura cuya resolución es más sencilla. Esta técnica se
basa en la dualización de esas ecuaciones y en la formulación de un problema
dual cuya resolución proporciona una aproximación del valor óptimo del
problema original mejor que la resolución de su relajación lineal12. En caso de
que el problema sea lineal y que el conjunto de ecuaciones de complicación sean
el conjunto de restricciones del problema, este problema dual es el problema
dual del problema original. Se formula el mismo subproblema que en DW y el
problema maestro es sustituido por su problema dual.
Partimos de la ecuación (5.15) que define el problema maestro completo
max θ2
θ2 ,π2
θ2 + (A1x 11 − b1 )T π2 ≤ c1T x 11 : λ1
θ2 + (A1x 12 − b1 )T π2 ≤ c1T x 12 : λ2 (5.20)
⋮
θ2 + (A1x 1ν − b1 )T π2 ≤ c1T x 1ν : λν
y definimos θ2 = b1T π2 + µ , entonces esta ecuación se convierte en
max b1T π2 + µ
π2 ,µ
(A1x 11 )T π2 + µ ≤ c1T x 11 : λ1
(A1x 12 )T π2 + µ ≤ c1T x 12 : λ2 (5.21)
⋮
(A1x 1ν )T π2 + µ ≤ c1T x 1ν : λν
luego el problema maestro restringido es
Maestro Relajación Lagrangiana
12
Por relajación lineal de un problema P se entiende el problema P en el que las variables
enteras son sustituidas por variables continuas.
04/11/2012 71
6 RELAJACIÓN LAGRANGIANA
max θ2
θ2 ,π2
(5.22)
θ2 + (A1x 1l − b1 )T π2 ≤ c1T x 1l : λl l = 1, …, j
o bien
max b1T π2 + µ
π2 ,µ
(5.23)
(A1x 1l )T π2 + µ ≤ cT1 x 1l : λl l = 1, …, j
donde las restricciones reciben el nombre de cortes duales o cortes de
optimalidad lagrangianos o corte de Lagrange. Esta formulación del maestro es
denominada método de planos de corte de Kelly.
I.6.2 Cortes de acotamiento

Dado que el maestro de RL es el dual del de DW, si el segundo es infactible
el primer resulta no acotado.
Consideremos un valor del multiplicador π2 y consideremos de nuevo el
subproblema de Lagrange
θ2 (π2 ) = min(c1T − πT2 A1 )x 1 + πT2 b1
x1
A2x 1 = b2 (5.24)
x1 ≥ 0
Supongamos que para este valor π2 del multiplicador existe v ≥ 0 y

A2v ≤ 0 13 tal que c1T − πT2 A1 < 0 . Entonces, si x 10 es una solución del
subproblema de Lagrange, observando que tv ≥ 0 y A2tv ≤ 0 para cualquier
valor positivo de t se tiene que
c1T (x 10 + tv ) + πT2 (A2 (x 10 + tv ) − b2 ) = cT1 x 10 + πT2 (A2x 10 − b2 ) + t(c1T + πT2 A2 )v (5.25)

con lo que el subproblema no está acotado dado que toma valores muy pequeños
cuando t → ∞ .
Del razonamiento anterior se deduce que el conjunto de multiplicadores
aceptables debe verificar que
(c1T + πT2 A2 )x 1 ≥ 0 ∀x 1 ∈ {x 1 ≥ 0, A2x 1 ≤ 0} (5.26)
La condición anterior, una vez propuesto un valor para el multiplicador,

puede ser comprobada resolviendo el siguiente problema, que denominamos
subproblema de acotamiento:
13
Es decir, v ∈ {x1 ≥ 0, A2x1 ≤ 0} que es el sistema homogéneo asociado a
v ∈ {x1 ≥ 0, A2x1 = b2 } .
72 04/11/2012
θ2* (π2 ) = min(c1T − πT2 A1 )x 1

x1
A2x 1 ≤ 0 (5.27)
0 ≤ x1 ≤ 1
y en caso de que tenga un valor negativo, se debe ser introducir un corte de

acotamiento (similar al corte de infactibilidad en la descomposición de Bd) en el
problema maestro de la RL de la forma:
(A1x 1l )T π2 ≤ c1T x 1l (5.28)

De este modo, el algoritmo de la RL itera entre un problema maestro,
formado por cortesduales y de acotamiento, y un subproblema que evalúa los
multiplicadores propuestos por el maestro. El problema maestro, considerando
ambos tipos de cortes, puede ser formulado de la siguiente forma
max θ2
θ2 ,π2
(5.29)
δ1l θ2 + (A1x 1l − δ1lb1 )T π2 ≤ c1T x 1l : λl l = 1, …, j
o bien
max b1T π2 + µ
π2 ,µ
(5.30)
(A1x 1l )T π2 + δ1l µ ≤ c1T x 1l : λl l = 1, …, j
siendo δ1l = 1 para los cortes duales y δ1l = 0 para los de acotamiento.
El subproblema es
Subproblema Relajación Lagrangiana
θ2 = min(c1T − π2jT A1 )x 1 + π2jTb1

x1
A2x 1 = b2 (5.31)
x1 ≥ 0
o bien
θ2 = min (c1T − π2jT A1 ) x 1 − µ j

x1
A2x 1 = b2 (5.32)
x1 ≥ 0
I.6.3 Algoritmo
Esquemáticamente el algoritmo se representa en la figura 4.5.
04/11/2012 73
SUBPROBLEMA
| ↑
j j
precios sombra π , µ 2 propuesta x 1j
↓ |
PROBLEMA MAESTRO RESTRINGIDO
Figura 4.5 Método de descomposición de Relajación Lagrangiana.
El algoritmo de la RL es descrito en los siguientes pasos.
1. Inicialización: j = 0 , z = ∞ , z = −∞ , ε = 10−4 .
2. Resolución del problema maestro de la RL
max θ2
θ2 ,π2
(5.33)
δ1l θ2 + (A1x 1l − δ1lb1 )T π2 ≤ c1T x 1l : λl l = 1, …, j
Actualizar la cota superior z = θ2 . Obtener el valor de π2 e ir al paso 3.
3. Resolución del subproblema de acotamiento
θ2* (π2 ) = min(c1T − π2jT A1 )x 1
x1
A2x 1 ≤ 0 (5.34)
0 ≤ x1 ≤ 1
Si θ2* (π2 ) ≥ 0 ir al paso 4. En otro caso obtener la solución x 1j y formar el

corte de acotamiento
(A1x 1j )T π2 ≤ c1T x 1j (5.35)

Ir al paso 2.
4. Resolución del subproblema
θ2 = min(c1T − π2jT A1 )x 1 + π2jTb1
x1
A2x 1 = b2 (5.36)
x1 ≥ 0
Actualizar la cota inferior z = θ2 . Obtener la solución x 1j y formar el corte

de Lagrange
θ2 + (A1x 1j − b1 )T π2 ≤ c1T x 1j (5.37)

Ir al paso 5.
5. Regla de parada
74 04/11/2012
Cálculo de la diferencia entre los multiplicadores propuestos en la iteración

j y la anterior j − 1 . Si d (π2j − π2j −1 ) < ε detener14. En otro caso, ir al paso 2.
La regla de parada se puede también establecer para la diferencia relativa
entre las cotas superior z e inferior z , establecidas por el maestro y el
subproblema respectivamente, o bien para el valor de la función objetivo del
subproblema en la segunda formulación.
La solución óptima es la combinación lineal de las soluciones obtenidas en

todas las iteraciones siendo los pesos las variables duales de los cortes de
Lagrange del maestro.
j
x 1* = ∑ x 1l λl (5.38)
l =1
j
z 1* = ∑ (cT1 x 1l )λl (5.39)
l =1
El tamaño del problema lineal completo es (m1 + m2 ) × n1 . El tamaño del

problema maestro es j × (m1 + 1) y el del subproblema m2 × n1 .
En el algoritmo de la RL, las variables π2 y µ son variables libres pero
puede ser conveniente introducir cotas algorítmicas suficientemente amplias a
sus valores que impiden que éste sea no acotado.
La función objetivo del maestro disminuye en cada iteración al introducir un
nuevo corte y, por consiguiente, será monótona decreciente. Por esta misma
razón se ocasiona una pérdida de factibilidad, luego conviene que sea resuelto
por el método simplex dual. Los cortes duales son una linealización exterior a la
función objetivo del subproblema.
El subproblema se ve afectado únicamente en su función objetivo luego se
recomienda utilizar el método simplex primal. El subproblema es no
diferenciable (es una poligonal).
Cuando la RL se aplica a problemas lineales, la sucesión de los valores
óptimos del problema maestro (que siempre es decreciente) converge al valor
óptimo del problema lineal. Por el contrario, esto no se satisface cuando el
problema es entero mixto. Este valor al que la sucesión decreciente converge es
simplemente una cota inferior del valor óptimo del problema. Esta diferencia
entre el valor óptimo del problema y el valor obtenido por la RL es lo que se
conoce en la literatura como intervalo de dualidad (duality gap). La aparición del
intervalo de dualidad origina que la solución primal propuesta por la RL pueda
ser infactible para las restricciones de complicación. Una alternativa para
obtener una solución factible para dichas restricciones es resolverlo mediante un
14
d (π2j − π2j −1 ) < ε representa distancia entre π2j y π2j −1 .
04/11/2012 75
método de penalización de programación no lineal como es el del lagrangiano

aumentado consiguiendo además diferenciabilidad en la función objetivo. Éste
consiste en introducir un término cuadrático que penaliza las infactibilidades
asociadas a las ecuaciones relajadas. Otra alternativa para evitar esta situación
consiste normalmente en realizar un postproceso de las soluciones obtenidas
para encontrar la solución factible. Este postproceso depende del problema
concreto que se resuelva y suele estar basado en el conocimiento específico de
dicho problema.
Dado que los problemas maestros de DW y RL son duales entre sí, ¿qué
ventaja puede tener uno frente a otro? La razón principal estriba en que para
DW las variables duales, que son cruciales en el algoritmo de descomposición,
son un subproducto de la optimización mientras que para RL son las variables
principales. De hecho pueden utilizarse formulaciones alternativas a la del
método de los planos de corte. Por ejemplo, una basada en técnicas de
programación no lineal puede ser el método del subgradiente
π2j +1 = π2j + αj p j (5.40)
donde p j = A1x 1j − b1 es el gradiente de la función dual y αj es un escalar que

indica la longitud de paso del método del subgradiente. αj disminuye al
aumentar las iteraciones y αj → 0 para j → ∞ . Un valor eficaz utilizado en la
práctica es
β j θ2 (π2j ) − c1T x 1* 

αj = 2 (5.41)
 
∑  j
∑ amn x n − bm 
m  n 
donde β j es un escalar entre 0 y 2, c1T x 1* es el valor de la función objetivo para
la mejor solución factible conocida del problema. Frecuentemente, se empieza la
secuencia por αj = 2 y se reduce su valor a la mitad cuando θ2 (π2j ) no ha
aumentado en un número prespecificado de iteraciones.
El bundle method es una variante del método de planos de corte donde la
función objetivo introduce un término cuadrático de penalización de la
desviación de los multiplicadores π2j con respecto a un punto πˆ2j que representa
el centro de gravedad de la región factible. En este caso el problema maestro
tiene función objetivo no lineal y además requiere el ajuste del parámetro de
penalización.
Incluso el problema maestro puede ser sustituido por algoritmos heurísticos
de actualización de los multiplicadores con tal de que se cumplan ciertas
condiciones. Por ejemplo, el método lambda de despacho económico en un
76 04/11/2012
sistema de energía eléctrica se puede interpretar como un algoritmo de RL

donde el multiplicador lambda se actualiza heurísticamente en cada iteración.

Volvamos a resolver el mismo problema que en el método de descomposición
de DW
min− 4x − y − 6z
3x + 2y + 4z = 17
1≤x ≤2
1≤y ≤2
1≤z ≤2
Empezamos resolviendo el subproblema con la variable dual π = 0 y el
óptimo resulta ser (x , y, z ) = (2,2, 2) y se pasa este punto y el valor de la función
objetivo, que es la cota inferior z = −22 , al problema maestro
max θ
θ + (18 − 17) π ≤ −22
Este problema resultaría no acotado. Para evitarlo se fija −1000 ≤ θ ≤ 1000
y −1000 ≤ π ≤ 1000 y se obtiene π = −1000 . La cota superior ahora es
z = 978 .
min (−4 + 3 ⋅ 1000) x + (−1 + 2 ⋅ 1000) y + (−6 + 4 ⋅ 1000) z − 17 ⋅ 1000
1≤x ≤2
1≤y ≤2
1≤z ≤2
El óptimo es el punto (x , y, z ) = (1,1,1) , el valor de la función objetivo
z = −8011 y un valor para el maestro de −11 . El nuevo problema maestro es
max θ
θ + (18 − 17) π ≤ −22
θ + (9 − 17) π ≤ −11
de donde resulta π = −1.222 y cota superior z = −20.777 lo que da lugar al
subproblema que se formula ahora como
04/11/2012 77
min− 0.333x + 1.444y − 1.111z − 20.777

1≤x ≤2
1≤y ≤2
1≤z ≤2
z = −22.222 y con un valor para el maestro de −21 . El nuevo problema
maestro es
max θ
θ + (18 − 17) π ≤ −22
θ + (9 − 17) π ≤ −11
θ + (16 − 17) π ≤ −21
de donde resulta π = −0.5 , variables duales de los cortes (λ1, λ2 , λ3 ) = (0.5, 0, 0.5)
y cota superior z = −21.5 . El subproblema es ahora
min− 2.5x − 4z − 8.5
1≤x ≤2
1≤y ≤2
1≤z ≤2
z = −21.5 . Al ser iguales la cota superior e inferior el algoritmo converge.
(x , y, z ) = (2,1.5, 2) .
Veamos de nuevo la resolución del problema con la formulación alternativa
de la RL. Igual que antes empezamos resolviendo el subproblema con las
variables duales π = 0 y µ = 0 y el óptimo resulta ser (x , y, z ) = (2,2, 2) y se
pasa este punto y el valor de la función objetivo, que en este caso es la cota
inferior z = −22 , al problema maestro
max 17 π + µ
18π + µ ≤ −22
Este problema resultaría no acotado. Para evitarlo se fija −1000 ≤ π ≤ 1000
y −1000 ≤ µ ≤ 1000 y se obtiene π = −56.777 y µ = 1000 . La cota superior
ahora es z = 34.777 .
78 04/11/2012
min (−4 + 3 ⋅ 56.777) x + (−1 + 2 ⋅ 56.777) y + (−6 + 4 ⋅ 56.777) z − 1000

1≤x ≤2
1≤y ≤2
1≤z ≤2
El óptimo es el punto (x , y, z ) = (1,1,1) , el valor de la función objetivo −500 ,
la cota inferior es z = −465.222 y un valor para el maestro de −11 . El nuevo
problema maestro es
max 17 π + µ
18π + µ ≤ −22
9π + µ ≤ −11
de donde resulta π = −1.222 y µ = 0 y cota superior z = −20.777 lo que da
lugar al subproblema que se formula ahora como
min− 0.333x + 1.444y − 1.111z
1≤x ≤2
1≤y ≤2
1≤z ≤2
El nuevo punto óptimo es (x , y, z ) = (2,1, 2) , el valor de la función objetivo
−1.444 , la cota inferior z = −22.222 y con un valor para el maestro de −21 .
El nuevo problema maestro es
max 17 π + µ
18π + µ ≤ −22
9π + µ ≤ −11
16π + µ ≤ −21
de donde resulta π = −0.5 y µ = −13 , variables duales de los cortes
(λ1, λ2 , λ3 ) = (0.5, 0, 0.5) y cota superior z = −21.5 . El subproblema es ahora
min− 2.5x − 4z + 13
1≤x ≤2
1≤y ≤2
1≤z ≤2
El nuevo punto óptimo es (x , y, z ) = (2,1, 2) , el valor de la función objetivo es
0 y la cota inferior z = −21.5 . El algoritmo converge en los dos criterios: por
ser iguales la cota superior e inferior o por ser 0 el valor del subproblema que
evalúa los costes reducidos.
04/11/2012 79

(x , y, z ) = (2,1.5, 2) .
La implantación en GAMS de un código de RL es similar a la implantación
de la descomposición de Bd. Se formula un modelo para el problema maestro y
otro para el subproblema y éstos son resueltos sucesivamente en un bucle de
iteraciones. Las ecuaciones que representan los cortes de Lagrange (una única
ecuación que engloba a los dos tipos de cortes) son declaradas sobre un conjunto
estático de índices y construidas sobre un subconjunto dinámico. Durante el
proceso algorítmico, este conjunto dinámico se va actualizado, aumentando por
tanto el número de cortes de Lagrange que contiene el problema maestro. El
algoritmo escrito en GAMS se presenta a continuación.
$title Descomposición de relajación lagrangiana (RL)

sets
L máximo número de iteraciones / iter-1 * iter-5 /
IT iteraciones del algoritmo / iter-1 * iter-9 /
N1 variables
scalar
Z_INF cota inferior de la solución / -inf /
Z_SUP cota superior de la solución / inf /
TOL tolerancia en optimalidad / 1e-6 /
GAMMA coeficiente cotas restricciones subproblema / 1 /
* fase 1 obtiene cortes hasta que el maestro resulte acotado
INFACT fase 1 (0) o fase 2 (1) de relajación lagrangiana / 0 /
N_CORTES número de cortes / 0 /
FORMUL tipo de formulación relajación lagrangiana / 2 /
parameters
DELTA1(l) coeficiente para cortes de acotamiento (subproblema no acotado)
X1_L(l,n1) soluciones del subproblema
sets
M1 restricciones primera etapa / runo-1 /
M2 restricciones segunda etapa / rdos-1 * rdos-6 /
N1 variables / xuno-1 * xuno-3 /
parameters
C1(n1) coeficientes función objetivo
/ xuno-1 -4
xuno-2 -1
xuno-3 -6 /
/ runo-1 17 /
/ rdos-1 -1
rdos-2 -1
rdos-3 -1
rdos-4 2
rdos-5 2
rdos-6 2 /
xuno-1 xuno-2 xuno-3
runo-1 3 2 4
xuno-1 xuno-2 xuno-3
rdos-1 -1
rdos-2 -1
rdos-3 -1
rdos-4 1
rdos-5 1
rdos-6 1
80 04/11/2012

variables
TT función objetivo maestro
PI2(m1) variables duales restricciones primera etapa
MU variable dual de la combinación lineal
equations
FOM1 función objetivo maestro
FOS1 función objetivo subproblema
CORTES1(l) cortes de relajación lagrangiana
CORTES2(l) cortes de relajación lagrangiana ;
FOM1 .. Z1 =E= sum
sum(m1, B1(m1)*PI2(m1)) + MU ;
FOM2 .. Z1 =E= TT ;
FOS1 .. Z2 =E= - sum

sum(m1, PI2.L(m1)*(sum
sum(n1, A1(m1,n1)*X1(n1)))) - MU.L +
sum
INFACT*sum
sum(n1, C1(n1)*X1(n1)) ;
sum
FOS2 .. Z2 =E= - sum
sum(m1, PI2.L(m1)*(sum
sum(n1, A1(m1,n1)*X1(n1))-B1(m1)*GAMMA)) +
sum
INFACT*sum
sum(n1, C1(n1)*X1(n1)) ;
sum
FOC .. Z2 =E= sum
sum(n1, C1(n1)*X1(n1)) ;
R1(m1) .. sum
sum(n1, A1(m1,n1) * X1(n1)) =L= B1(m1) ;
R2(m2) .. sum
sum(n1, A2(m2,n1) * X1(n1)) =L= B2(m2) * GAMMA ;
CORTES1(j) .. sum
sum(m1, PI2(m1)*(sum
sum(n1, A1(m1,n1)*X1_L(j,n1))))
sum +
DELTA1(j)*MU =L= INFACT*sum
sum(n1, C1(n1)*X1_L(j,n1)) ;
sum
CORTES2(j) .. sum
sum(m1, PI2(m1)*(sum
sum(n1, A1(m1,n1)*X1_L(j,n1))-DELTA1(j)*B1(m1))) +
sum
DELTA1(j)*TT =L= INFACT*sum
sum(n1, C1(n1)*X1_L(j,n1)) ;
sum
model MAESTRO1 / FOM1, CORTES1 /
model MAESTRO2 / FOM2, CORTES2 /
model SUB1 / FOS1, R2 /
model SUB2 / FOS2, R2 /
model COMPLETO / FOC , R1, R2 /

SUB1.OptFile = 1 ; SUB2.OptFile = 1 ;
MAESTRO1.OptFile = 1 ; MAESTRO2.OptFile = 1 ;
J(l) = NO ;
X1_L(l,n1) = 0 ;
DELTA1(l) = 1 ;
Z1.L = 1e3 ;
* inicialización de variables duales para subproblema

MU.L = 0 ;
PI2.L(m1) = 0 ;
* cotas de variables maestro en fase 1 (para acotar maestro)
TT.LO = - 1e3 ;
TT.UP = 1e3 ;
PI2.LO(m1) = - 1 ;
PI2.UP(m1) = 0 ;
MU.LO = - 1 ;
MU.UP = 1 ;
loop(it
loop $(INFACT = 0 or Z_INF < Z_SUP),
if (ord
ord(it)
ord > 1,
if (FORMUL = 1,
04/11/2012 81
solve MAESTRO1 USING LP MAXIMIZING Z1

else
solve MAESTRO2 USING LP MAXIMIZING Z1
) ;
Z_SUP = Z1.L
) ;
display Z_SUP
* se resuelve el subproblema salvo si en fase 1 el maestro es 0

if (INFACT = 1 or Z_SUP > TOL,
DELTA1(l) $(ord
ord(l) = N_CORTES + 1) = 1 ;
ord
if (FORMUL = 1,
solve SUB1 USING LP MINIMIZING Z2 ;
Z_INF $(SUB1.MODELSTAT = 1) = Z2.L + MU.L + sum
sum(m1, PI2.L(m1)*B1(m1)) ;
* cálculo del rayo de forma correcta
if (SUB1.MODELSTAT = 3 or SUB1.MODELSTAT = 19,
GAMMA = 0 ;
X1.LO(n1) = -1 ;
X1.UP(n1) = 1 ;
GAMMA = 1 ;
X1.LO(n1) = -inf ;
X1.UP(n1) = inf ;
Z_INF = -inf ;
DELTA1(l) $(ord
ord(l) = N_CORTES + 1) = 0 ;
ord
) ;
else
Z_INF $(SUB2.MODELSTAT = 1) = Z2.L ;
* cálculo del rayo de forma correcta
if (SUB2.MODELSTAT = 3 or SUB2.MODELSTAT = 19,
GAMMA = 0 ;
X1.LO(n1) = -1 ;
X1.UP(n1) = 1 ;
GAMMA = 1 ;
X1.LO(n1) = -inf ;
X1.UP(n1) = inf ;
Z_INF = -inf ;
DELTA1(l) $(ord
ord(l) = N_CORTES + 1) = 0 ;
ord
) ;
) ;
display Z_INF
* sólo se genera un corte si el valor del maestro es superior al del subproblema
if(Z_SUP
if > Z_INF,
X1_L(l,n1) $(ord
ord(l) = N_CORTES + 1) = X1.L(n1) ;
ord
J(l) $(ord
ord(l) = N_CORTES + 1) = yes ;
ord
N_CORTES = N_CORTES + 1 ;
) ;
) ;
* paso a la fase 2 de la RL cuando z1.l=0
if (INFACT = 0 and Z_SUP < TOL,
INFACT = 1 ;
PI2.LO(m1) = -inf ;
PI2.UP(m1) = 0 ;
MU.LO = -inf ;
MU.UP = inf ;
Z_SUP = 1e3 ;
Z_INF = -inf ;
);
) ;
abort $(ABS(Z_INF-Z_SUP) > TOL) 'Máximo número de iteraciones alcanzado' ;
* resultados descomposición
if (FORMUL = 1,
X1.L(n1) = sum
sum(j, X1_L(j,n1)*CORTES1.M(j))
else
X1.L(n1) = sum
sum(j, X1_L(j,n1)*CORTES2.M(j))
) ;

82 04/11/2012
Descomposición de relajación lagrangiana (RL)

E x e c u t i o n
---- 145 PARAMETER Z_SUP = +INF cota superior de la solución
Equation Listing SOLVE SUB2 Using LP From line 166
---- FOS2 =E= función objetivo subproblema
FOS2.. Z2 =E= 0 ; (LHS = 0)
R2(rdos-1).. - X1(xuno-1) =L= -1 ; (LHS = 0, INFES = 1 ***)
R2(rdos-4).. X1(xuno-1) =L= 2 ; (LHS = 0)
R2(rdos-5).. X1(xuno-2) =L= 2 ; (LHS = 0)
R2(rdos-6).. X1(xuno-3) =L= 2 ; (LHS = 0)
Model Statistics SOLVE SUB2 Using LP From line 166
LOOPS IT iter-1
MODEL STATISTICS
NON ZERO ELEMENTS 7
Solution Report SOLVE SUB2 Using LP From line 166
L O O P S IT iter-1
MODEL SUB2 OBJECTIVE Z2
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOS2 . . . 1.0000
rdos-1 -INF -1.0000 -1.0000 EPS
rdos-2 -INF -1.0000 -1.0000 EPS
rdos-3 -INF -1.0000 -1.0000 EPS
rdos-4 -INF 1.0000 2.0000 .
rdos-5 -INF 1.0000 2.0000 .
rdos-6 -INF 1.0000 2.0000 .
xuno-1 -INF 1.0000 +INF .
0 INFEASIBLE
0 UNBOUNDED
04/11/2012 83
E x e c u t i o n
---- 181 PARAMETER Z_INF = 0.000 cota inferior de la solución
Equation Listing SOLVE MAESTRO2 Using LP From line 141
---- FOM2 =E= función objetivo maestro
FOM2.. Z1 - TT =E= 0 ; (LHS = 1000, INFES = 1000 ***)
---- CORTES2 =L= cortes de relajación lagrangiana
CORTES2(iter-1).. TT - 8*PI2(runo-1) =L= 0 ; (LHS = 0)
Model Statistics SOLVE MAESTRO2 Using LP From line 141
LOOPS IT iter-2
MODEL STATISTICS
NON ZERO ELEMENTS 4
Solution Report SOLVE MAESTRO2 Using LP From line 141
L O O P S IT iter-2
MODEL MAESTRO2 OBJECTIVE Z1
TYPE LP DIRECTION MAXIMIZE
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOM2 . . . 1.0000
---- EQU CORTES2 cortes de relajación lagrangiana
iter-1 -INF . . 1.0000
---- VAR TT -1000.0000 . 1000.0000 .
---- VAR PI2 variables duales restricciones primera etapa
runo-1 -1.0000 . . 8.0000
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
---- 145 PARAMETER Z_SUP = 0.000 cota superior de la solución
FOM2.. Z1 - TT =E= 0 ; (LHS = 0)
CORTES2(iter-1).. TT - 8*PI2(runo-1) =L= -11 ; (LHS = 0, INFES = 11 ***)
LOOPS IT iter-3
84 04/11/2012
MODEL STATISTICS
NON ZERO ELEMENTS 4
L O O P S IT iter-3
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOM2 . . . 1.0000
iter-1 -INF -11.0000 -11.0000 1.0000
---- VAR Z1 -INF -11.0000 +INF .
---- VAR TT -1000.0000 -11.0000 1000.0000 .
runo-1 -INF . . 8.0000
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
---- 145 PARAMETER Z_SUP = -11.000 cota superior de la solución
FOS2.. 4*X1(xuno-1) + X1(xuno-2) + 6*X1(xuno-3) + Z2 =E= 0 ; (LHS = 11, INFES = 11 ***)
R2(rdos-1).. - X1(xuno-1) =L= -1 ; (LHS = -1)
R2(rdos-2).. - X1(xuno-2) =L= -1 ; (LHS = -1)
R2(rdos-3).. - X1(xuno-3) =L= -1 ; (LHS = -1)
R2(rdos-4).. X1(xuno-1) =L= 2 ; (LHS = 1)
R2(rdos-5).. X1(xuno-2) =L= 2 ; (LHS = 1)
R2(rdos-6).. X1(xuno-3) =L= 2 ; (LHS = 1)
LOOPS IT iter-3
MODEL STATISTICS
04/11/2012 85

L O O P S IT iter-3
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOS2 . . . 1.0000
rdos-1 -INF -2.0000 -1.0000 .
rdos-2 -INF -2.0000 -1.0000 .
rdos-3 -INF -2.0000 -1.0000 .
rdos-4 -INF 2.0000 2.0000 -4.0000
rdos-5 -INF 2.0000 2.0000 -1.0000
rdos-6 -INF 2.0000 2.0000 -6.0000
---- VAR Z2 -INF -22.0000 +INF .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
FOM2.. Z1 - TT =E= 0 ; (LHS = 0)
CORTES2(iter-1).. TT - 8*PI2(runo-1) =L= -11 ; (LHS = -11)
CORTES2(iter-2).. TT + PI2(runo-1) =L= -22 ; (LHS = -11, INFES = 11 ***)
LOOPS IT iter-4
MODEL STATISTICS
NON ZERO ELEMENTS 6
L O O P S IT iter-4
86 04/11/2012

scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOM2 . . . 1.0000
iter-1 -INF -11.0000 -11.0000 0.1111
iter-2 -INF -22.0000 -22.0000 0.8889
---- VAR Z1 -INF -20.7778 +INF .
---- VAR TT -1000.0000 -20.7778 1000.0000 .
runo-1 -INF -1.2222 . .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
FOS2.. 0.333333333333334*X1(xuno-1) - 1.44444444444444*X1(xuno-2) + 1.11111111111111*X1(xuno-3) + Z2 =E=
-20.7777777777778 ; (LHS = -22, INFES = 1.22222222222222 ***)
R2(rdos-1).. - X1(xuno-1) =L= -1 ; (LHS = -2)
R2(rdos-2).. - X1(xuno-2) =L= -1 ; (LHS = -2)
R2(rdos-3).. - X1(xuno-3) =L= -1 ; (LHS = -2)
R2(rdos-4).. X1(xuno-1) =L= 2 ; (LHS = 2)
R2(rdos-5).. X1(xuno-2) =L= 2 ; (LHS = 2)
R2(rdos-6).. X1(xuno-3) =L= 2 ; (LHS = 2)
LOOPS IT iter-4
MODEL STATISTICS
L O O P S IT iter-4
scaind -1
lpmethod 1
preind 0
epopt 1e-9
04/11/2012 87
eprhs 1e-9
---- EQU FOS2 -20.7778 -20.7778 -20.7778 1.0000
rdos-1 -INF -2.0000 -1.0000 .
rdos-2 -INF -1.0000 -1.0000 -1.4444
rdos-3 -INF -2.0000 -1.0000 .
rdos-4 -INF 2.0000 2.0000 -0.3333
rdos-5 -INF 1.0000 2.0000 .
rdos-6 -INF 2.0000 2.0000 -1.1111
---- VAR Z2 -INF -22.2222 +INF .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
FOM2.. Z1 - TT =E= 0 ; (LHS = 0)
CORTES2(iter-1).. TT - 8*PI2(runo-1) =L= -11 ; (LHS = -11)
CORTES2(iter-2).. TT + PI2(runo-1) =L= -22 ; (LHS = -22)
CORTES2(iter-3).. TT - PI2(runo-1) =L= -21 ; (LHS = -19.5555555555556, INFES = 1.44444444444444 ***)
LOOPS IT iter-5
MODEL STATISTICS
NON ZERO ELEMENTS 8
L O O P S IT iter-5
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOM2 . . . 1.0000
88 04/11/2012

iter-1 -INF -17.5000 -11.0000 .
iter-2 -INF -22.0000 -22.0000 0.5000
iter-3 -INF -21.0000 -21.0000 0.5000
---- VAR Z1 -INF -21.5000 +INF .
---- VAR TT -1000.0000 -21.5000 1000.0000 .
runo-1 -INF -0.5000 . .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
FOS2.. 2.5*X1(xuno-1) + 4*X1(xuno-3) + Z2 =E= -8.5 ; (LHS = -9.22222222222222, INFES = 0.722222222222221 ***)
R2(rdos-1).. - X1(xuno-1) =L= -1 ; (LHS = -2)
R2(rdos-2).. - X1(xuno-2) =L= -1 ; (LHS = -1)
R2(rdos-3).. - X1(xuno-3) =L= -1 ; (LHS = -2)
R2(rdos-4).. X1(xuno-1) =L= 2 ; (LHS = 2)
R2(rdos-5).. X1(xuno-2) =L= 2 ; (LHS = 1)
R2(rdos-6).. X1(xuno-3) =L= 2 ; (LHS = 2)
LOOPS IT iter-5
MODEL STATISTICS
NON ZERO ELEMENTS 9
L O O P S IT iter-5
scaind -1
lpmethod 1
preind 0
epopt 1e-9
eprhs 1e-9
---- EQU FOS2 -8.5000 -8.5000 -8.5000 1.0000
rdos-1 -INF -2.0000 -1.0000 .
rdos-2 -INF -1.0000 -1.0000 EPS
rdos-3 -INF -2.0000 -1.0000 .
rdos-4 -INF 2.0000 2.0000 -2.5000
rdos-5 -INF 1.0000 2.0000 .
rdos-6 -INF 2.0000 2.0000 -4.0000
04/11/2012 89

---- VAR Z2 -INF -21.5000 +INF .
0 INFEASIBLE
0 UNBOUNDED
E x e c u t i o n
FOC.. 4*X1(xuno-1) + X1(xuno-2) + 6*X1(xuno-3) + Z2 =E= 0 ; (LHS = 0)
R1(runo-1).. 3*X1(xuno-1) + 2*X1(xuno-2) + 4*X1(xuno-3) =L= 17 ; (LHS = 17)
R2(rdos-1).. - X1(xuno-1) =L= -1 ; (LHS = -2)
R2(rdos-2).. - X1(xuno-2) =L= -1 ; (LHS = -1.5)
R2(rdos-3).. - X1(xuno-3) =L= -1 ; (LHS = -2)
R2(rdos-4).. X1(xuno-1) =L= 2 ; (LHS = 2)
R2(rdos-5).. X1(xuno-2) =L= 2 ; (LHS = 1.5)
R2(rdos-6).. X1(xuno-3) =L= 2 ; (LHS = 2)
MODEL STATISTICS
---- EQU FOC . . . 1.0000
runo-1 -INF 17.0000 17.0000 -0.5000
rdos-1 -INF -2.0000 -1.0000 .
rdos-2 -INF -1.5000 -1.0000 .
rdos-3 -INF -2.0000 -1.0000 .
rdos-4 -INF 2.0000 2.0000 -2.5000
rdos-5 -INF 1.5000 2.0000 .
rdos-6 -INF 2.0000 2.0000 -4.0000
90 04/11/2012
---- VAR Z2 -INF -21.5000 +INF .

0 INFEASIBLE
0 UNBOUNDED
COPT C:\cplex.opt
**** FILE SUMMARY
Input C:\RL.gms
Output C:\RL.lst
I.6.5 Problema de transporte con coste fijo

Consideremos el problema de transporte con coste fijo anterior y
supongamos que son las restricciones de capacidad límite las que son relajadas e
introducidas en la función objetivo a través del lagrangiano. Para un valor λij el
subproblema de Lagrange presenta la forma:
min ∑ (cij x ij + fij yij ) + λij (x ij − M ij yij )

x ij ,yij
ij
∑xj
ij ≤ ai ∀i
∑xi
ij ≥ bj ∀j
x ij ≥ 0, yij ∈ {0,1}
que reformulado queda como
min ∑ (cij + λij ) x ij + ( fij − λij M ij ) yij

x ij ,yij
ij
∑x j
ij ≤ ai ∀i
∑x i
ij ≥ bj ∀j
x ij ≥ 0, yij ∈ {0,1}
Esta reformulación permite observar que el subproblema es separable en dos
problemas. Un problema de transporte en el que el coste variable es ligeramente
modificado por el multiplicador y un segundo problema que es entero puro y
cuya solución puede obtenerse de modo inmediato.
Problema de transporte
04/11/2012 91
min ∑ (cij + λij ) x ij

x ij
ij
∑x j
ij ≤ ai ∀i
∑x ij ≥ bj ∀j
i
x ij ≥ 0
Problema entero puro
min ∑ ( fij − λij M ij ) yij

yij
ij
yij ∈ {0,1}
Esta formulación de la RL para el problema del transporte con coste fijo
puede implantarse en GAMS con el siguiente código.
$title Problema de transporte con coste fijo (RL)

OPTION OPTCR = 0
sets
L índice de iteraciones / l1 * l100 /
LL(l) subconjunto de índices
I orígenes / i1 * i4 /
J destinos / j1 * j3 /
* Datos del problema
parameters
A(i) ofertas de producto
/ i1 10, i2 30, i3 40, i4 20 /
B(j) demandas de producto
/ j1 20, j2 50, j3 30 /
table C(i,j) coste variable unitario de transporte

j1 j2 j3
i1 1 2 3
i2 3 2 1
i3 2 3 4
i4 4 3 2
table F(i,j) coste fijo de transporte
j1 j2 j3
i1 10 20 30
i2 20 30 40
i3 30 40 50
i4 40 50 60
* Fin datos del problema
abort $(sum
sum[i, A(i)] < sum
sum sum[j, B(j)]) 'Problema infactible'
scalar
TOL tolerancia relativa / 1e-6 /
Z_INF cota inferior / -INF /
Z_SUP cota superior / inf /
DIF diferencia en iteraciones
parameters
X_J(i,j,l) valores de las variables de flujo en la iteración j
Y_J(i,j,l) valores de las variables de inversión en la iteración j
LAMBDA_J(i,j,l) multiplicadores en la iteración j
positive variables
X(i,j) flujo por los arcos
LAMBDA(i,j) multiplicador
92 04/11/2012
binary variable
Y(i,j) decisiones de inversion en los arcos
variables
Z variable objetivo primal (subproblema)
W variable dual
* Declaración de las ecuaciones

equations
EQ_OBJ función objetivo problema completo
OFERTA(i) ofertas de los orígenes
DEMANDA(j) demanda de los destinos
LIMITE(i,j) límite de uso del arco
CORTES_RL(l) cortes de Relajación Lagrangiana ;
EQ_OBJ .. Z =E= sum
sum[(i,j), F(i,j)*Y(i,j)] + sum
sum[(i,j), C(i,j)*X(i,j)]
+ sum
sum[(i,j), LAMBDA.L(i,j)*(X(i,j)-MIN[A(i),B(j)] * Y(i,j))] ;
OFERTA(i) .. sum
sum[j, X(i,j)] =L= A(i) ;
DEMANDA(j) .. sum
sum[i, X(i,j)] =G= B(j) ;
LIMITE(i,j) .. X(i,j) =L= MIN[A(i),B(j)] * Y(i,j) ;

CORTES_RL(ll).. W=L=sum
sum[(i,j), F(i,j)*Y_J(i,j,ll)] + sum
sum sum[(i,j), C(i,j)*X_J(i,j,ll)] +
sum[(i,j), LAMBDA(i,j)*(X_J(i,j,ll)-MIN[A(i),B(j)] * Y_J(i,j,ll))] ;
sum
model MAESTRO_RL / CORTES_RL /
model SUB_RL / EQ_OBJ, OFERTA, DEMANDA /
model COMPLETO / EQ_OBJ, OFERTA, DEMANDA, LIMITE / ;
X.UP(i,j) = MIN[A(i),B(j)]
* para que los resultados de un problema infactible sean los correctos con CPLEX sólo
* se puede utilizar el método simplex y sin preproceso
SUB_RL.OptFile = 1 ; COMPLETO.OptFile = 1 ;
* inicialización de parámetros del problema
ll(l) = NO ;
Y_J(i,j,l) = 0 ;
X_J(i,j,l) = 0 ;
LAMBDA.UP(i,j) = 10 ;
DIF = inf ;
* Iteraciones del algoritmo
loop(l
loop $(DIF > TOL),
* Resolución del problema maestro de Relajación Lagrangiana
if(ord
if ord(l)
ord > 1,
solve MAESTRO_RL USING LP MAXIMIZING W
else
LAMBDA.L(i,j) = 0
) ;
LAMBDA_J(i,j,l) = LAMBDA.L(i,j) ;
* Resolución del subproblema de Relajación Lagrangiana
solve SUB_RL USING MIP MINIMIZING Z ;
X_J(i,j,l) = X.L(i,j) ;
Y_J(i,j,l) = Y.L(i,j) ;
* Actualización del conjunto de cortes de Lagrange

ll(l) = yes ;
* Actualización de la diferencia de multiplicadores

DIF $(ord
ord(l)>1) = sum
ord sum[(i,j), ABS(LAMBDA_J(i,j,l)-LAMBDA_J(i,j,l-1))]
) ;
LAMBDA.FX(i,j) = 0 ;
solve COMPLETO USING MIP MINIMIZING Z
04/11/2012 93
En este ejemplo, la sucesión de valores del problema maestro y la sucesión de

óptimos del subproblema en la RL convergen en la iteración 69 a la solución
dada por la relajación lineal del problema15 con valor de la función objetivo de
355 u.m., ver figura 4.6. El valor de la función objetivo del problema completo
es de 380 u.m.
Cota inferior Cota superior
1500
1000
500
0
1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 52 55 58 61 64 67
-500
-1000
Iteración
Figura 4.6 Iteraciones del método de descomposición de Relajación Lagrangiana.
Por otra parte, los valores obtenidos por las variables de decisión no son
factibles en ninguna de las iteraciones del método. Se necesita un postproceso de
estas soluciones para obtener una solución. En este problema, una posibilidad es
considerar la solución dada para las variables continuas x ij por el problema de
transporte (lineal) y adecuar las variables binarias a esta solución, es decir,
si x ij > 0 entonces yij = 1

Con este postproceso de las soluciones se puede comprobar que la solución
obtenida es, en este caso, la solución óptima del problema entero mixto:
y1 = y23 = y 31 = y 32 = y 42 = 1
15
Cuando esto ocurre se dice que el problema satisface la propiedad de integralidad
(integrality property) y tiene poco sentido el uso de la RL como técnica de aproximación del
valor óptimo del problema.
94 04/11/2012
I.7 Otras técnicas de descomposición

I.7.1 Descomposición Primal-Dual
El método de descomposición Primal-Dual [Krishna] nace con el objetivo de
evitar el crecimiento en la cota superior que suele aparecer en las últimas
iteraciones del proceso de descomposición de Bd. Para suprimirlo en lugar de
utilizar en el subproblema el último valor de x 1 obtenido por el maestro se
formula la mejor combinación lineal de los valores obtenidos en un subconjunto
de iteraciones anteriores. Este método tiene un comportamiento más robusto
que la descomposición de Bd frente a errores en el cálculo de sensibilidades por
tener en cuenta todas las propuestas previas que ha recibido del maestro. Esto
incrementa ligeramente el tamaño del subproblema pero acelera el proceso de
convergencia del algoritmo y disminuye sustancialmente el número de
iteraciones. Es una generalización de los métodos de descomposición previos,
DW y Bd.
El problema maestro restringido para la iteración j se plantea como en el
método de Bd:
Maestro Primal-Dual
min cT1 x 1 + θ2
x1 ,θ2
A1x 1 = b1
(6.1)
π2lT B1x 1 + θ2 ≥ π2lTb2 l = 1, …, j
x1 ≥ 0
El subproblema para cada iteración j ahora se formula de manera similar al

maestro en DW:
Subproblema Primal-Dual
t
min ∑ (c1T x 1l )λl + cT2 x 2
λl ,x 2
l =1
t
l
∑ (B x )λ
l =1
1 1 l + A2x 2 = b2 : π2
(6.2)
t
∑λ
l =1
l =1 :µ
λl , x 2 ≥ 0 l = 1, …, j
donde t ≤ j asociado al tamaño del problema, siendo 10 un valor máximo
razonable. Se pueden intentar diferentes estrategias de combinación lineal de
04/11/2012 95
7 OTRAS TÉCNICAS DE DESCOMPOSICIÓN
soluciones anteriores x 1l y actualización durante el proceso iterativo cuando

t < j . La estrategia actualmente utilizada es substituir aquella solución con el
menor valor de λl por la nueva solución x 1l y la de mayor valor de λl por una
combinación de todas las t soluciones.
Si λl = 0 para cierta l ello indica que la solución x 1l no forma parte de la
combinación lineal óptima y puede ser eliminada. Eliminando solamente estas
soluciones previas se garantiza que la cota superior es monótona decreciente. Si
λl = 1 se detiene la ejecución del proceso iterativo.
El tamaño del problema lineal bietapa PL-2 es (m1 + m2 ) × (n1 + n2 ) . El
tamaño del problema maestro es (m1 + j ) × (n1 + 1) y el del subproblema
(m2 + 1) × (t + n2 ) . El número de restricciones del maestro aumenta en uno en
cada iteración. El tamaño del subproblema aumenta con respecto al del método
de Bd en t variables y 1 restricción. El subproblema del método Primal-Dual se
convierte en el de Bd si t = 1 y en el maestro de DW si t = j .
Nótese que las variables duales del subproblema del método Primal-Dual
serán diferentes, en general, de las variables duales del subproblema del método
de Bd.
I.7.2 Descomposición anidada

La descomposición anidada [Birge:85] consiste en aplicar de forma recursiva
el principio de descomposición [Dantzig:63]. Si se aplica al problema lineal
multietapa PL-P definido como
P
min ∑ cTp x p
xp
p =1
Bp−1x p−1 + Ap x p = bp p = 1, …, P (6.3)

xp ≥ 0
B0 ≡ 0
la etapa 1 es el problema maestro y las etapas 2, …, P el subproblema. Para
solucionar el subproblema se aplica de nuevo el principio de descomposición
siendo la etapa 2 el problema maestro y las etapas 3, …, P el subproblema. De
nuevo se aplica descomposición hasta llegar a la etapa P − 1 como problema
maestro y la etapa P como subproblema. En cada momento sólo se necesita
resolver un problema correspondiente a una etapa, como problema maestro para
las etapas p = 1, …, P − 1 o como subproblema para la etapa P . Este método
resuelve repetidamente una secuencia de problemas lineales de menor tamaño
para solucionar el problema lineal multietapa.
La descomposición anidada se puede utilizar con cualquiera de los métodos
de descomposición anteriores. Cuando se utiliza el método de descomposición de
96 04/11/2012
Bd una etapa cualquiera p en su papel de maestro recibe la información para

generar los cortes correspondientes a los subproblemas p + 1, …, P y pasa la
solución a dichos subproblemas. Debido a la estructura en escalera, la
información para generar los cortes procede sólo del subproblema p + 1 y sólo la
cota de dicho subproblema se ve afectada por la solución de la etapa p . De
forma similar, la etapa p en su función de primer subproblema de la serie
p, …, P genera cortes para su maestro, la etapa p − 1 , y ve modificada su cota
por la información procedente de su maestro p − 1 .
Veamos un ejemplo de problema lineal de 4 etapas.
min c1T x 1 + cT2 x 2 + cT3 x 3 + cT4 x 4
x1 ,x 2 ,x 3 ,x 4
A1x 1 = b1
B1x 1 +A2x 2 = b2
(6.4)
B2x 2 +A3x 3 = b3
B 3x 3 +A4x 4 = b4
x 1, x 2, x 3, x4 ≥0
Se resuelve el problema de la etapa 4 para un cierto valor de x 3l

obteniéndose el valor de la función objetivo θ4 y las variables duales
correspondientes π4 . Por claridad a partir de ahora se suprime el superíndice l
en las variables duales que forman parte de una restricción.
min cT4 x 4
x4
A4x 4 = b4 − B3x l3 : π4 (6.5)

x4 ≥ 0
El problema maestro restringido de la etapa 3 resulta ser:
min cT3 x 3 + θ4
x 3 ,θ4
A3x 3 = b3 − B2x 2l : π3
(6.6)
θ4 + πT4 B3x 3 ≥ πT4 b4 : η3
x3 ≥ 0
siendo la función objetivo del problema dual max πT3 (b3 − B2x 2l ) + ηT3 (πT4 b4 )
π3 ,η3
Si en lugar de resolver este maestro se hubiera resuelto el problema de las
dos últimas etapas tendríamos:
04/11/2012 97
min cT3 x 3 + cT4 x 4

x 3 ,x 4
A3x 3 = b3 − B2x 2l : π3′

(6.7)
B3x 3 + A4x 4 = b4 : µ3
x 3, x 4 ≥ 0
siendo la función objetivo del dual max π3′T (b3 − B2x 2l ) + µT3 b4
π3′ ,µ3
El problema maestro restringido de la etapa 3 es una cota inferior de la
función objetivo del problema de las dos últimas etapas. Luego los óptimos
correspondientes a los problemas duales verifican la misma relación, esto es:
π3′T (b3 − B2x 2l ) + µT3 b4 ≥ πT3 (b3 − B2x 2l ) + ηT3 (πT4 b4 ) (6.8)
siendo esta relación válida para cualquier valor de x 2
π3′T (b3 − B2x 2 ) + µT3 b4 ≥ πT3 (b3 − B2x 2 ) + ηT3 (πT4 b4 ) (6.9)
El corte que se introduciría en la etapa anterior si el problema se hubiera

considerado bietapa (las dos primeras por una parte y las dos últimas por otra)
sería
θ3 ≥ π3′T (b3 − B2x 2 ) + µT3 b4 ≥ πT3 (b3 − B2x 2 ) + ηT3 (πT4 b4 ) (6.10)
Esta ecuación indica que el corte pasado a la etapa anterior es válido aunque
desde la etapa actual hasta la última no se haya resuelto el problema hasta
optimalidad [Morton:93]. Ahora la expresión del corte de Benders queda
θ3 + πT3 B2x 2 ≥ πT3 b3 + ηT3 (πT4 b4 ) (6.11)
El problema maestro restringido de la etapa 2 resulta ser:

min cT2 x 2 + θ3
x 2 ,θ3
A2x 2 = b2 − B1x 1l : π2
(6.12)
θ3 + πT3 B2x 2 ≥ πT3 b3 + ηT3 (πT4 b4 ) : η2
x2 ≥ 0
siendo la función objetivo del problema dual
max π2 (b2 − B1x 1 ) + η2  π3 b3 + η3 (π4 b4 )
T l T  T T T
π2 ,η2
Si en lugar de resolver este maestro se hubiera resuelto el problema de las
tres últimas etapas tendríamos:
98 04/11/2012
min cT2 x 2 + cT3 x 3 + cT4 x 4

x 2 ,x 3 ,x 4
A2x 2 = b2 − B1x 1l : π2′

B2x 2 + A3x 3 = b3 : µ2 (6.13)
B3x 3 + A4x 4 = b4 : µ3
x 2, x 3 , x 4 ≥ 0
siendo la función objetivo del dual max π2′T (b2 − B1x 1l ) + µT2 b3 + µT3 b4
π2′ ,µ2 ,υ2
El corte a introducir en la etapa 1 será:
θ2 ≥ π2′T (b2 − B1x 1 ) + µT2 b3 + µT3 b4 ≥ πT2 (b2 − B1x 1 ) + ηT2  πT3 b3 + ηT3 (πT4 b4 ) (6.14)
 
θ2 + πT2 B1x 1 ≥ πT2 b2 + ηT2  πT3 b3 + ηT3 (πT4 b4 ) (6.15)
 
Viendo simultáneamente los cortes para las diferentes etapas se puede
deducir la expresión general del corte en cualquier etapa:
θ4 + πT4 B3x 3 ≥ πT4 b4
θ3 + πT3 B2x 2 ≥ πT3 b3 + ηT3 (πT4 b4 ) (6.16)
θ2 + πT2 B1x 1 ≥ πT2 b2 + ηT2  πT3 b3 + ηT3 (πT4 b4 )
θp +1 + πTp +1Bp x p ≥ q p = πTp +1bp +1 + ηTp +1q p +1 (6.17)

Si en lugar de estos cortes se utilizan los de la linealización alrededor de un
punto se tiene:
θ4 + πT4 B3x 3 ≥ f4l + πT4 B3x l3
(6.18)
θ4 ≥ f4l + πT4 B3 (x 3l − x 3 )
θ3 ≥ πT3 (b3 − B2x 2 ) + ηT3 (πT4 b4 ) =

= πT3 (b3 − B2x 2 + B2x 2l − B2x 2l ) + ηT3 (πT4 b4 ) =
(6.19)
= πT3 (b3 − B2x 2l ) + πT3 B2 (x 2l − x 2 ) + ηT3 (πT4 b4 ) =
= f3l + πT3 B2 (x 2l − x 2 )
y la expresión genérica para cualquier etapa es:
θp +1 ≥ fpl +1 + πTp +1Bp (x lp − x p )

(6.20)
θp +1 + πTp +1Bpx p ≥ fpl +1 + πTp +1Bp x lp
La formulación del problema a resolver para cualquier etapa p , p = 1, …, P
es la siguiente:
04/11/2012 99
Subproblema descomposición anidada
min cTp x p + θp +1
x p ,θp +1
Ap x p = bp − Bp−1x lp−1 : πp
πlT lT lT
p +1Bp x p + θp +1 ≥ q p = πp +1bp +1 + ηp +1q p +1 : ηp l = 1, …, j
xp ≥ 0
(6.21)
θP +1 ≡ 0
B0 ≡ 0
πPl +1 ≡ 0
ηPl +1 ≡ 0
o bien con la otra formulación de los cortes:
min cTp x p + θp +1
x p ,θp +1
Ap x p = bp − Bp−1x lp−1 : πp
πlT l lT l
p +1Bp x p + θp +1 ≥ fp +1 + πp +1Bp x p : ηp l = 1, …, j
xp ≥ 0
(6.22)
θP +1 ≡ 0
B0 ≡ 0
πPl +1 ≡ 0
ηPl +1 ≡ 0
donde el superíndice l de los vectores significa que son valores de variables
conocidos en alguna iteración previa y θp +1 ∈ ℝ .
En términos generales, el algoritmo consiste en resolver cierta etapa p y
entonces o bien proceder con la siguiente etapa p + 1 modificando las cotas o
proceder con la etapa previa p − 1 formando un corte. Existen varias estrategias
para recorrer las etapas del problema que se explican más adelante. El
procedimiento iterativo finaliza cuando se produce la convergencia en la primera
etapa.
I.7.3 Caso ejemplo

Volvamos a resolver el mismo problema que en el método de descomposición
de DW
100 04/11/2012
min− 4x − y − 6z
3x + 2y + 4z = 17
1≤x ≤2
1≤y ≤2
1≤z ≤2
Para ello lo primero es poner las ecuaciones de manera que sean diagonal por
bloques introduciendo una variable auxiliar v = 2y + 4z
min− 4x − y − 6z
3x + v = 17
v − 2y − 4z = 0
1≤x ≤2
1≤y ≤2
1≤z ≤2
Empezamos resolviendo el problema maestro
min− 4x
1≤x ≤2
que da lugar a x = 2 . Esa solución se le pasa al siguiente subproblema
min− y
v = 17 − 3 ⋅ 2
1≤y ≤2
que da lugar a y = 2 y v = 11 . Esa solución se le pasa al siguiente subproblema
min− 6z
−4z = −11 + 2 ⋅ 2
1≤z ≤2
que da lugar a z = 1.75 , función objetivo −10.5 y la variable dual asociada a la
primera ecuación es 1.5 . Esa solución se le pasa al subproblema anterior para
formar un corte
min− y + θ3
v = 17 − 3 ⋅ 2
θ3 + 10.5 ≥ 1.5(11 − 2 ⋅ 2 − v + 2y )
1≤y ≤2
04/11/2012 101
que da lugar a y = 1 , función objetivo −14.5 y la variable dual asociada a la

primera ecuación es −1.5 . Esa solución se le pasa al subproblema anterior para
formar un corte
min− 4x + θ2
θ2 + 14.5 ≥ −1.5(2 − x )
1≤x ≤2
que da lugar a x = 2 , función objetivo z = −22.5 , que es una cota inferior del
problema. La cota superior será z = −20.5 . Como hay diferencia entre cotas se
debe proseguir con las iteraciones.
I.7.4 Modelo hidrotérmico

A continuación se presenta un caso ejemplo escrito en GAMS de operación
de un sistema hidrotérmico resuelto por descomposición anidada. Las
restricciones que acoplan las etapas corresponden a las de gestión del nivel de
reserva de cada embalse. Es interesante observar en la programación las
siguientes características:
• El modelo se ha formulado de manera que tanto el problema completo como

cada subproblema sean resueltos con exactamente las mismas ecuaciones.
Por una parte, las ecuaciones se definen sobre un conjunto dinámico en
etapas. Por otra, hay un escalar que controla las diferencias entre una y otra
formulación.
• La ecuación de gestión hidráulica interetapa es única, incluyendo las
condiciones iniciales y finales de situación de los embalses que son prefijadas
y la formulación simultánea para el problema completo y descompuesto.
• El algoritmo de descomposición recorre las etapas en sentido descendente de
2 hasta P (excepto en la primera etapa donde también se resuelve en
sentido descendente la etapa 1) y en sentido ascendente desde P − 1 a 1.
$title Modelo de planificación de la expansión y explotación de generación y red

* Andrés Ramos
* Departamento de Organización Industrial
* Escuela Técnica Superior de Ingeniería
* Universidad Pontificia Comillas
* Junio 2000
* especificación de la división temporal del modelo
* dimensionamiento del sistema generación red
* atributos del sistema
sets
T años / ann-1 /
K(t) subconjunto de años / ann-1 /
P periodos / per-1 * per-3 /
QQ(p) subconjunto de periodos
B bloques / blq-1 * blq-3 /
ND nudos / nudo-1 * nudo-9 /
102 04/11/2012
TR generadores térmicos / gent-1 * gent-7 /

HD grupos hidráulicos / hidr-1 * hidr-4 /
GD gestión de demanda / gstd-1 * gstd-2 /
ATG atributos de los generadores o acciones de gestión de demanda
/ pmx, pmn, cvr, efor, q1, cfj, tcf, rmx, rmn, rvai, rvaf /
ATL atributos de las líneas / r, x, flmx /
NDTR(nd,tr) localización de generadores térmicos en nudos
/
nudo-1 . gent-4
nudo-1 . gent-5
nudo-2 . gent-3
nudo-2 . gent-6
nudo-4 . gent-2
nudo-7 . gent-1
nudo-9 . gent-7
/
NDHD(nd,hd) localización de generadores hidráulicos en nudos
/
nudo-1 . hidr-3
nudo-5 . hidr-4
nudo-7 . hidr-1
nudo-9 . hidr-2
/
LE(nd,nd) conjunto de líneas existentes
L iteraciones de Benders / 1 * 5 /
SNT sentido en la pasada de la descomposición / 1aP, Pa1 /
ALIAS (nd,ni,nf), (p,pp)
* parámetros generales
scalar
SBASE potencia base (GW) / 0.1 /
MR margen de reserva (p.u.) / 0.1 /
CENS coste ENS (€ por kWh) / 300 /
PNNR penalización por energía (€ por kWh) / 1e3 /
PNPT penalización por potencia (€ por kW) / 1e3 /
CSL coeficiente seguridad líneas / 0.7 /
OPCPRD opción de pérdidas no lineales / 0 /
OPCDA opción de descomposición anidada / 0 /
TMSB tiempo de ejecución del subproblema / 0 /
ZINF cota inferior / -inf /
ZSUP cota superior / inf /
TOL tolerancia en la descomposición / 1e-9 /
CSTT coste completo / 0 /
parameters
THETA_L(l,t,p) función de recursos de cada periodo e iteración
RVA_L(l,t,p,hd) reserva de cada hidráulica en cada periodo e iteración
PIRVA_L(l,t,p,hd) variable dual de la reserva hidráulica
ORDEN(p) orden de cada periodo
DSP(p) desplazamiento del índice
INCDEM(t) incremento acumulado de la demanda (p.u.)
INC(t) incremento anual de la demanda (p.u.) / ann-1 0.00 /
PCTND(nd) reparto de la demanda total por nudos (p.u.)
/ nudo-1 0.10
nudo-2 0.10
nudo-3 0.10
nudo-4 0.05
nudo-5 0.20
nudo-6 0.05
nudo-7 0.05
nudo-8 0.20
nudo-9 0.15 / ;
ORDEN(p) = ord(p)
ord ;
INCDEM(t) = PROD(k
PROD $(ord
ord(k) <= ord(t)),
ord ord 1+INC(k)) ;
table DTGT(tr,atg) datos generadores térmicos

pmx pmn cvr efor q1 cfj tcf
* MW MW € p.u. p.u. € p.u.
* kWh W
gent-1 800 800 1.0 0.00 0.00 300 0.10
gent-2 800 800 1.0 0.00 0.00 300 0.10
gent-3 250 100 3.0 0.05 0.05 180 0.12
gent-4 400 250 3.5 0.05 0.05 192 0.12
gent-5 300 100 4.0 0.05 0.05 146 0.13
04/11/2012 103
gent-6 300 100 4.0 0.05 0.05 146 0.13

gent-7 300 100 4.0 0.05 0.05 146 0.13
table DTHD(hd,atg) datos generadores hidráulicos
pmx rmx rmn rvai rvaf cfj tcf
* MW GWh GWh GWh GWh € p.u.
* W
hidr-1 600 2500 100 1250 1250 240 0.12
hidr-2 600 2500 200 1250 1250 230 0.12
hidr-3 500 2500 300 1250 1250 240 0.12
hidr-4 500 2500 400 1250 1250 220 0.12
table DTGD(gd,atg) datos acciones de gestión de demanda
cvr cfj tcf
* € € p.u.
* kWh W
gstd-1 5.0 100 0.13
gstd-2 3.5 150 0.14
table APORT(hd,p) aportaciones en grupos hidráulicos (GWh)
per-1 per-2 per-3
hidr-1 525.0 700.0 175.0
hidr-2 450.0 600.0 150.0
hidr-3 412.5 550.0 137.5
hidr-4 412.5 550.0 137.5
table DEM(b,p) demanda del bloque b en el periodo p (MW)

per-1 per-2 per-3
blq-1 5000 3500 3000
blq-2 4000 2800 2400
blq-3 2400 1680 1440
table DUR(b,p) duración del bloque b en el periodo p (h)
per-1 per-2 per-3
blq-1 166.66 333.33 166.66
blq-2 500.00 1000.0 500.00
blq-3 333.33 666.66 333.33
table DTLE(ni,nf,atl) datos líneas existentes
r x flmx
* p.u. p.u. MW
nudo-1 . nudo-6 0.002 0.02 800
nudo-2 . nudo-3 0.004 0.03 800
nudo-2 . nudo-6 0.006 0.04 800
nudo-3 . nudo-4 0.008 0.05 800
nudo-3 . nudo-6 0.007 0.06 800
nudo-4 . nudo-5 0.005 0.07 800
nudo-4 . nudo-6 0.003 0.08 800
nudo-4 . nudo-9 0.001 0.07 800
nudo-6 . nudo-7 0.003 0.06 800
nudo-6 . nudo-8 0.005 0.05 800
nudo-7 . nudo-8 0.007 0.04 800
nudo-8 . nudo-9 0.006 0.03 800 ;
LE(ni,nf) $DTLE(ni,nf,'x') = yes ;
* Escalación de los datos de generadores y líneas de MW a GW y GWh a TWh

APORT(hd,p) = APORT(hd,p) / 1e3 ;
DTHD(hd,'rmn') = DTHD(hd,'rmn') / 1e3 ;
DTHD(hd,'rmx') = DTHD(hd,'rmx') / 1e3 ;
DTHD(hd,'rvai') = DTHD(hd,'rvai') / 1e3 ;
DTHD(hd,'rvaf') = DTHD(hd,'rvaf') / 1e3 ;
DTGT(tr,'pmx') = DTGT(tr,'pmx') / 1e3 ;
DTGT(tr,'pmn') = DTGT(tr,'pmn') / 1e3 ;
DTHD(hd,'pmx') = DTHD(hd,'pmx') / 1e3 ;
DEM(b,p) = DEM(b,p) / 1e3 ;
DTLE(ni,nf,'flmx') = DTLE(ni,nf,'flmx') / 1e3 ;
variables
ACP(t,p,b,tr) decisión de acoplamiento de los grupos térmicos (0 o 1)
PTP(t,p,b,tr) potencia térmica producida (GW)
PHP(t,p,b,hd) potencia hidráulica producida (GW)
RVA(t,p,hd) reserva hidráulica (GWh)
MOVHD(p,hd) movimiento del embalse (GWh)
PDM(t,p,b,nd,gd) potencia de gestión de demanda (GW)
PNS(t,p,b,nd) potencia no suministrada (GW)
EXC(t,p,b,nd) exceso de potencia (GW)
DRR(t,p) defecto de reserva rodante (GW)
ERR(t,p) exceso de reserva rodante (GW)
104 04/11/2012
TT(t,p,b,nd) ángulo del nudo (rad)

FLE(t,p,b,ni,nf) flujo en líneas existentes (GW)
PTI(t,tr) potencia térmica instalada (GW)
PHI(t,hd) potencia hidráulica instalada (GW)
PDI(t,nd,gd) potencia de gestión demanda instalada (GW)
PRDAS(t,p,b,nd) pérdidas en el nudo (GW)
CSTTL coste total inversión (gen y red) y explotación (G€)
THETA(t,p) función de recursos ;
positive variables ACP, PTP, PHP, RVA, PDM, PNS, EXC, DRR, ERR, PTI, PHI, PDI
equations
FO Costes totales de inversión y explotación
KR1(t,p,b,nd) Primera ley de Kirchhoff para cada nudo
EPRDAS(t,p,b,nd) Pérdidas en cada nudo
FLJE(t,p,b,ni,nf) Flujo líneas existentes
ACOPLA(t,p,b,tr) Acoplamiento de los grupos térmicos
PRDTRM(t,p,b,tr) Potencia térmica producida inferior a la acoplada
PTPMIN(t,p,tr) Potencia térmica producida mínima
PRDHID(t,p,b,hd) Potencia hidráulica producida inferior a la instalada
RSRVAS(t,p,hd) Reservas hidráulicas
MOVEMHD(t,p,hd) Movimiento de embalses hidráulicos
PRDDMN(t,p,b,nd,gd) Gestión de demanda utilizada inferior a la instalada
MARGEN(t,p) Margen de reserva
CORTES(l,t,p) Cortes de Benders ;
FO .. CSTTL =E= sum

sum((t,qq), THETA(t,qq)) $(OPCDA = 1) +
sum((t,tr),
sum PTI(t,tr) * DTGT(tr,'cfj')*DTGT(tr,'tcf') * (card
card(t)-ord
card ord(t)+1))
ord
+ sum((t,hd),
sum PHI(t,hd) * DTHD(hd,'cfj')*DTHD(hd,'tcf') * (card
card(t)-ord
card ord(t)+1))
ord
+ sum((t,nd,gd),
sum PDI(t,nd,gd) * DTGD(gd,'cfj')*DTGD(gd,'tcf') * (card
card(t)-ord
card ord(t)+1))
ord
+ sum((t,qq,b,tr),
sum DTGT(tr,'cvr') * PTP(t,qq,b,tr) * DUR(b,qq))/1e3
+ sum((t,qq,b,nd,gd),
sum DTGD(gd,'cvr') * PDM(t,qq,b,nd,gd) * DUR(b,qq))/1e3
+ CENS * sum((t,qq,b,nd),
sum PNS(t,qq,b,nd) * DUR(b,qq))/1e3
+ PNNR * sum((t,qq,b,nd),
sum EXC(t,qq,b,nd) * DUR(b,qq))/1e3
+ PNPT * sum((t,qq),
sum DRR(t,qq) + ERR(t,qq))/1e3 ;
KR1(t,qq,b,nd) ..
sum(NDTR(nd,tr),
sum PTP(t,qq,b,tr)) + sum
sum(NDHD(nd,hd), PHP(t,qq,b,hd))
+ sum(gd,
sum PDM(t,qq,b,nd,gd))
+ sum(LE(ni,nd),
sum FLE(t,qq,b,ni,nd)) - sum
sum(LE(nd,nf), FLE(t,qq,b,nd,nf))
+ PNS(t,qq,b,nd) - EXC(t,qq,b,nd)
=E= PCTND(nd) * DEM(b,qq) * INCDEM(t) + PRDAS(t,qq,b,nd) $OPCPRD ;
EPRDAS(t,qq,b,nd) .. PRDAS(t,qq,b,nd) =E=
+ SBASE * sum(LE(ni,nd),
sum (1-cos(TT(t,qq,b,ni) - TT(t,qq,b,nd)))
* DTLE(ni,nd,'r') / (DTLE(ni,nd,'r')**2+DTLE(ni,nd,'x')**2))
+ SBASE * sum(LE(nd,nf),
sum (1-cos(TT(t,qq,b,nf) - TT(t,qq,b,nf)))
* DTLE(nd,nf,'r') / (DTLE(nd,nf,'r')**2+DTLE(nd,nf,'x')**2)) ;
FLJE(t,qq,b,LE(ni,nf)) ..
FLE(t,qq,b,ni,nf)/SBASE =E= (TT(t,qq,b,ni)-TT(t,qq,b,nf))/DTLE(ni,nf,'x') ;
PRDTRM(t,qq,b,tr) ..
PTP(t,qq,b,tr) =L= (DTGT(tr,'pmx') +
sum(k
sum $(ord
ord(k) <= ord(t)),
ord ord PTI.L(k,tr))) * (1-DTGT(tr,'efor')) * ACP(t,qq,b,tr) ;
ACOPLA(t,qq,b+1,tr) .. ACP(t,qq,b+1,tr) =L= ACP(t,qq,b,tr) ;

PTPMIN(t,qq,tr) ..
PTP(t,qq,'blq-3',tr) =G= DTGT(tr,'pmn')*(1 - DTGT(tr,'q1'))*ACP(t,qq,'blq-3',tr) ;
PRDHID(t,qq,b,hd) ..
PHP(t,qq,b,hd) =L= DTHD(hd,'pmx') + sum
sum(k $(ord
ord(k) <= ord(t)),
ord ord PHI.L(k,hd)) ;
* restricción para la descomposición
RSRVAS(t,qq,hd) ..
sum(pp
sum $SAMEAS(pp+1,qq), RVA(t,pp,hd)) $(OPCDA = 0)
+ sum(pp
sum $SAMEAS(pp+1,qq), RVA.L(t,pp,hd)) $(OPCDA = 1)
+ DTHD(hd,'rvai') $(ORDEN(qq) = 1) + MOVHD(qq,hd) =E=
RVA(t,qq,hd) $(ORDEN(qq) <> card
card(p)) + DTHD(hd,'rvaf') $(ORDEN(qq) = card
card(p)) ;
MOVEMHD(t,qq,hd) ..
MOVHD(qq,hd) =E= APORT(hd,qq) - sum
sum(b, PHP(t,qq,b,hd) * DUR(b,qq)) / 1e3 ;
PRDDMN(t,qq,b,nd,gd) .. PDM(t,qq,b,nd,gd) =L= sum
sum(k $(ord
ord(k) <= ord(t)),
ord ord PDI(k,nd,gd)) ;
MARGEN(t,qq) ..
sum(tr,
sum (DTGT(tr,'pmx') + sum
sum(k $(ord
ord(k) <= ord
ord ord(t)), PTI.L(k,tr)))
* (1-DTGT(tr,'efor')) * ACP(t,qq,'blq-1',tr))
04/11/2012 105
+ sum(hd,
sum DTHD(hd,'pmx') + sum
sum(k $(ord
ord(k) <= ord(t)),
ord ord PHI(k,hd)))
+ sum((nd,gd),
sum sum(k $(ord
sum ord(k) <= ord(t)),
ord ord PDI(k,nd,gd)))
+ DRR(t,qq) - ERR(t,qq)
=E= DEM('blq-1',qq) * (1+MR) * INCDEM(t) ;
CORTES(j,t,qq) $(OPCDA = 1) ..
THETA(t,qq) - THETA_L(j,t,qq) =G=
sum(hd,
sum PIRVA_L(j,t,qq,hd) * (RVA(t,qq,hd) - RVA_L(j,t,qq,hd))) ;
* Se anula la expansión de la generación para que el modelo sea lineal
PTI.FX(t,tr) = 0 ;
PHI.FX(t,hd) = 0 ;
PDI.FX(t,nd,gd) = 0 ;
* cotas a las variables
RVA.LO(t,p,hd) = DTHD(hd,'rmn') ;
RVA.UP(t,p,hd) = DTHD(hd,'rmx') ;
ACP.UP(t,p,b,tr) = 1 ;
PTP.UP(t,p,b,tr) = DEM(b,p) * INCDEM(t) ;
PHP.UP(t,p,b,hd) = DEM(b,p) * INCDEM(t) ;
PDM.UP(t,p,b,nd,gd) = PCTND(nd) * DEM(b,p) * INCDEM(t) ;
PNS.UP(t,p,b,nd) = PCTND(nd) * DEM(b,p) * INCDEM(t) ;
EXC.UP(t,p,b,nd) = PCTND(nd) * DEM(b,p) * INCDEM(t) ;
DRR.UP(t,p) = (1+MR) * DEM('blq-1',p) * INCDEM(t) ;
ERR.UP(t,p) = (1+MR) * DEM('blq-1',p) * INCDEM(t) ;
PTI.UP(t,tr) = smax((b,p), DEM(b,p))
smax * INCDEM(t) ;
PHI.UP(t,hd) = smax((b,p), DEM(b,p))
smax * INCDEM(t) ;
PDI.UP(t,nd,gd) = PCTND(nd) * smax
smax((b,p), DEM(b,p)) * INCDEM(t) ;
TT.LO(t,p,b,nd) = -1 ;
TT.UP(t,p,b,nd) = 1 ;
TT.FX(t,p,b,'nudo-1') = 0 ;
FLE.LO(t,p,b,ni,nf) = - DTLE(ni,nf,'flmx') * CSL ;
FLE.UP(t,p,b,ni,nf) = DTLE(ni,nf,'flmx') * CSL ;
model PLGNRD / FO, KR1, FLJE, ACOPLA, PRDTRM, PTPMIN, PRDHID, RSRVAS, PRDDMN, MARGEN,
CORTES / ;
model PLGNRDNL / FO, KR1, FLJE, ACOPLA, PRDTRM, PTPMIN, PRDHID, RSRVAS, PRDDMN, MARGEN,
CORTES, EPRDAS / ;
file SAL / sal.txt /
put SAL ;
J(l) = NO ;
THETA_L(l,t,p) = 0 ;
RVA_L(l,t,p,hd) = 0 ;
PIRVA_L(l,t,p,hd) = 0 ;
* proceso iterativo de descomposición anidada de Benders
OPCDA = 1 ;
QQ(p) = NO ;
loop (l $(ABS(1-ZINF/ZSUP) > TOL and ZINF < ZSUP),
* iteraciones hacia adelante snt=1 y hacia atrás snt=2

loop (snt,
if (ord
ord(snt) = 1,
ord
DSP(p) = 0 ;
else
DSP(p) = card
card(p) - 2*ord
ord(p) + 1 ;
ord
) ;
loop (p,
QQ(p+DSP(p)) = yes ;
* solución del problema en 2 a P y P-1 a 1
if ((ord
ord(l) = 1 and ord(snt)
ord ord = 1) or (ord
ord(snt)=1 and ord(p)+DSP(p)
ord ord > 1) or
(ord
ord(snt)=2 and ord(p)+DSP(p)
ord ord < card
card(p)),
if (ord
ord(snt) = 1 and ord
ord ord(p) = card
card(p),
THETA.FX(t,qq) = 0 ;
elseif ord
ord(l) = 1 and ord
ord(snt) = 1,
THETA.LO(t,qq) = 0 ;
else
THETA.LO(t,qq) = -inf ;
106 04/11/2012
) ;
solve PLGNRD USING LP MINIMIZING CSTTL ;
put 'Suproblema (M€)'
sum(qq, ORDEN(qq)):3:0 CSTTL.L:12:6 '
sum tiempo (s)' PLGNRD.resusd:7:1
PLGNRD.numequ:7:0 PLGNRD.numvar:7:0 PLGNRD.numnz:7:0 PLGNRD.iterusd:7:0/;
TMSB = TMSB + PLGNRD.resusd ;
abort $(PLGNRD.modelstat > 2) 'Terminación anormal' ;
) ;
* incremento de iteraciones al llegar al ultimo hacia adelante
J(l) $(ord
ord ord(p) = card
card(p)) = yes ;
* se guardan las variables interperiodo
if (ord
ord ord(p)+DSP(p) < card
card(p),
RVA_L(l,t,qq,hd) = RVA.L(t,qq,hd) ;
) ;
* las variables duales se guardan en sentido ascendente
if ((ord
ord(snt) = 1 and ord(p)+DSP(p)
ord ord = card
card(p)) or
(ord
ord(snt) = 2 and ord(p)+DSP(p)
ord ord < card
card(p)),
PIRVA_L(l,t,qq,hd) = - RSRVAS.M(t,qq,hd) ;
THETA_L(l,t,pp) $SAMEAS(pp+1,p+DSP(p)) = CSTTL.L ;
) ;
* cálculo de los costes totales
if (ord
ord(snt) = 1 or (ord
ord ord(snt) = 2 and ord
ord ord(p)+DSP(p) = 1),
CSTT = CSTT + CSTTL.L - sum
sum((t,qq), THETA.L(t,qq)) ;
) ;
QQ(p+DSP(p)) = NO ;
) ;
* calculo de las cotas inferior y superior
ZSUP $(ord
ord(snt)
ord = 1) = MIN(ZSUP,CSTT) ;
ZINF $(ord
ord(snt)
ord = 2) = MAX(ZINF,CSTTL.L) ;
put $(ord
ord(snt) = 1) /
ord
put $(ord
ord(snt) = 2) L.TL 'Cota Inf' ZINF:12:6 '
ord Cota Sup' ZSUP:12:6
' Converg' ABS(1-ZINF/ZSUP):10:6 ' tiempo (s)' TMSB:7:1 // ;
CSTT $(ord
ord(snt)
ord = 2) = 0 ;
) ;
) ;
* modelo completo
OPCDA = 0 ;
J(l) = NO ;
QQ(p) = yes ;
solve PLGNRD USING LP MINIMIZING CSTTL ;

put 'Completo (M€) ' CSTTL.L:12:6 ' tiempo (s)' PLGNRD.resusd:7:1
PLGNRD.numequ:7:0 PLGNRD.numvar:7:0 PLGNRD.numnz:7:0 PLGNRD.iterusd:7:0 // ;
abort $(PLGNRD.modelstat > 2) 'Terminación anormal' ;
04/11/2012 107
I.8 Optimización lineal estocástica bietapa y

multietapa
I.8.1 Bietapa
El problema lineal estocástico bietapa PLE-2, donde se considera
aleatoriedad en la segunda etapa, se formula matemáticamente como16:
minω cT1 x 1 + ∑ p ωc2ωT x 2ω

x1 ,x 2
ω ∈Ω
A1x 1 = b1
(7.1)
B xω
1 1 +A x ω ω
2 2 = b2ω
x 1, x 2ω ≥0
donde la función objetivo corresponde a las decisiones de la primera etapa más

el valor esperado correspondiente a las decisiones de la segunda etapa y donde
ω ∈ Ω toma un conjunto finito Ω , cardinal de Ω , de valores cada uno con
probabilidad p ω . En un modelo general la estocasticidad puede aparecer en
cualquier término del problema lineal. En la práctica sólo un número reducido
de coeficientes suele ser estocástico. Si A2ω es independiente de ω este problema
se denomina con recurso fijo. La estructura de la matriz de restricciones de un
problema PLE-2 para tres escenarios en la segunda etapa se representa en la
siguiente figura.
A1
B11 A21
B12 A22
B13 A23
Figura 8.1 Estructura de la matriz de coeficientes de las restricciones en PLE-2.
El problema de minimización del máximo arrepentimiento tiene una

estructura similar al anterior y se formula como:
16
El superíndice ω indica la dependencia de las matrices y vectores con respecto al valor del
escenario ω .
04/11/2012 109
8 OPTIMIZACIÓN LINEAL ESTOCÁSTICA
min α
α,x1 ,x 2ω
α −c1T x 1 −c2ωT x 2ω ≥ −f ω ∀ω ∈ Ω
A1x 1 = b1 (7.2)
B1ω x 1 +A2ω x 2ω = b2ω
x 1, x 2ω ≥0
siendo f ω el valor de la función objetivo con información perfecta para cada

escenario ω (i.e., la función objetivo para la decisión óptima sabiendo que va a
ocurrir un determinado escenario ω ).
Otro problema diferente pero similar a éste en estructura sería minimizar el
máximo coste, que puede tener sentido si no se confía en las probabilidades
asociadas a los escenarios.
min α
α,x1 ,x 2ω
α −c1T x 1 −c2ωT x 2ω ≥0 ∀ω ∈ Ω
A1x 1 = b1 (7.3)
B1ω x 1 +A2ω x 2ω = b2ω
x 1, x 2ω ≥0
Si se supone que ω puede tomar tres posibles valores Ω = {ω1, ω2 , ω3 } ,

entonces el problema determinista equivalente [Dantzig:55] del (7.1), donde se
enumeran los valores de ω , es:
min c1T x 1 + p ω1c2ω1T x 2ω1 + p ω2 c2ω2T x 2ω2 + p ω3c2ω3T x 2ω3
x1 ,x 2 ,x 2 ,x 2ω3
ω1 ω2
A1x 1 = b1
B1ω1 x 1 +A2ω1 x 2ω1 = b2ω1
(7.4)
B1ω2 x 1 +A2ω2 x 2ω2 = b2ω2
B1ω3 x 1 +A2ω3 x 2ω3 = b2ω3
x 1, x 2ω1 , x 2ω2 , x 2ω3 , ≥0
La enumeración es posible siempre que el conjunto de posibles estados de ω ,

Ω , sea suficientemente pequeño. El tamaño del bloque de la segunda etapa
depende linealmente del número de escenarios. Por otra parte, si el número de
parámetros fuera 25 cada uno con 2 posibles estados, el número total de estados
sería de 225 ≈ 3.3 107 . El espacio de estados crece exponencialmente con el
número de parámetros aleatorios. En este caso, es necesario recurrir a técnicas
de simulación para obtener muestras de los mismos. En la sección siguiente a la
110 04/11/2012
planificación estocástica multietapa se explica el uso de la simulación en el

contexto de optimización estocástica.
El uso de técnicas de descomposición se justifica por la estructura especial
que presenta la matriz de restricciones del problema completo. Los
subproblemas son separables para cada escenario y tienen la misma estructura
en las restricciones. Por primera vez este problema equivalente se resolvió
aplicando el método de descomposición DW a su problema dual [Dantzig:61]
donde son las restricciones las que complican su resolución. Posteriormente se
utilizó la descomposición de Bd [VanSlyke] al problema primal original, donde
son las variables las que complican. Este último método tiene ventaja sobre el
anterior ya que el primal tiene más variables que restricciones (o dicho de otra
forma, el dual generalmente tiene más filas que el primal) y, por lo tanto,
requiere más cálculo por iteración [Birge:88].
Para su resolución utilizando la descomposición de Bd se considera
 p ω1c ω1  B ω1  Aω1  b ω1  x ω1 
 2   1   2   2   2 
 ω2 ω2   ω2     ω2   ω2 
c2 =  p c2  , B1 = B1  , A2 =  A2ω2  , b2 = b2  y x 2 = x 2 
      ω   ω 
p ω3 c ω3  B ω3   A2ω3 

b2 3  x 2 
3
 2   1  
Claramente, el subproblema es separable en tantos subproblemas como
estados de ω , siendo cada uno de ellos:
min
ω
c2ωT x 2ω
x2
A2ω x 2ω = b2ω − B1ω x 1l : π2ω (7.5)

x 2ω ≥ 0
La formulación del problema maestro restringido se puede hacer con una
única variable θ2 ∈ ℝ para todos los subproblemas, denominado método
monocorte:
min cT1 x 1 + θ2
x1 ,θ2
A1x 1 = b1
(7.6)
∑p ω
π2ωlT B1ω x 1 + θ2 ≥ ∑p ω
π2ωlTb2ω l = 1, …, j
ω ∈Ω ω ∈Ω
x1 ≥ 0
o múltiples variables θ2ω ∈ ℝ , una para cada subproblema, denominado método

multicorte:
04/11/2012 111
min
ω
cT1 x 1 + ∑ p ω θ2ω
x1 ,θ2
ω ∈Ω
A1x 1 = b1
(7.7)
π2ωlT B1ω x 1 + θ2ω ≥ π2ωlTb2ω ω ∈Ω l = 1, …, j
x1 ≥ 0
Este último se basa en la idea intuitiva de que utilizando explícitamente

todos los cortes se envía más información sobre los subproblemas que con uno y,
por lo tanto, se necesitarán menos iteraciones del algoritmo de descomposición.
Pero, por otra parte, el problema maestro del método multicorte se incrementa
en ( Ω − 1) variables y en cada iteración en ( Ω − 1) restricciones con respecto al
monocorte. El método multicorte será especialmente eficiente cuando m2 sea
grande y Ω no sea significativamente mayor que n1 [Birge:88]. Cuando Ω es
mucho mayor que n1 es conveniente utilizar un método híbrido consistente en
agrupar subconjuntos de cortes hasta obtener un número reducido de éstos. No
obstante, la ventaja de un método sobre otro depende del problema y debe
obtenerse experimentalmente.
En el caso de que sólo las cotas de las restricciones sean parámetros
estocásticos cualquier corte obtenido en cualquier nodo de un cierto periodo es
un corte válido para cualquier otro nodo del mismo periodo. La razón está en
que, para este caso, la región factible del problema dual de todos los nodos de
un mismo periodo es la misma y, por consiguiente, las variables duales obtenidas
en cualquier nodo no son más que diferentes17 vértices de un mismo politopo.
Todos los cortes de los nodos de un periodo dado son aproximaciones lineales
exteriores de una misma función de recursos.
I.8.2 Multietapa
El problema lineal estocástico multietapa PLE-P se puede formular como:
P
ωp ωpT ω
min
ωp ∑ ∑p p p c xp p
xp
p =1 ωp ∈Ωp
ωp ωp −1 ω ω ω
B x
p −1 p−1 + Ap p x p p = bp p p = 1, …, P (7.8)
ωp
x p ≥0
B ω1
0 ≡0
17
El cambio en la función objetivo del dual (por cambio en la cota de las restricciones del
primal) es lo que hace que se obtenga un vértice diferente aun en el caso de nodos con un mismo
ascendiente.
112 04/11/2012
donde Ωp representa el conjunto de los posibles valores de ωp en la etapa p

ω
cada uno con probabilidad de ocurrencia pp p . Estas probabilidades son, en
general, dependientes de los valores tomados en la etapa anterior. Cuando se
ω ω
van a tomar las decisiones al comienzo de la etapa p los valores de cp p , Bp−p 1 ,
ω ω
Ap p , bp p son conocidos y se conocen las distribuciones condicionadas de los
vectores para las etapas futuras p + 1, …, P . Los diferentes valores de los
parámetros en las sucesivas etapas forman una estructura en árbol.
Si Ωp representa el cardinal de Ωp , el problema lineal completo tiene
P P ω
∑ p =1
m p Ωp restricciones y ∑ p =1
n p Ωp variables, siendo Ap p una matriz
m p × n p . La matriz de restricciones de un problema PLE-3 con dos escenarios en
la segunda etapa y cuatro en la tercera se representa en la siguiente figura.
A1
B11 A21
B12 A22
B21 A31
B22 A32
B23 A33
B24 A34
Figura 4.16 Estructura de la matriz de coeficientes de las restricciones en problemas lineales multietapa.
El método de resolución es descomposición anidada estocástica de Benders

(también denominado en la literatura programación dinámica dual estocástica
[Pereira]) teniendo presente que los subproblemas para cada etapa son
separables. La forma intuitiva de resolver este árbol de problemas lineales será:
resolver el maestro de la primera etapa, pasar las cotas a las restricciones de la
segunda etapa, resolver los Ω2 subproblemas. Pasar estos valores como cotas de
la tercera etapa, resolver estos subproblemas, etc. llegando hasta la última etapa
y luego ir retrocediendo hasta la primera etapa formando los cortes
correspondientes. El método converge cuando converge el subproblema de la
primera etapa.
La formulación del subproblema a resolver para la etapa p es la siguiente
por el método del monocorte:
04/11/2012 113
ωT ω ω
min
ω ω
cp p x p p + θp +p 1
x p p ,θp +p 1
ω ω ω ω a ( ω )l ω
Ap p x p p = bp p -Bp−p 1x p−1p : πp p
(7.9)
p p πpklT+1Bpk x p p + θp +p 1 ≥ q p p = (πpklT+1bpk+1 + ηpklT+1q pk +1 )
ω ω ω ω ωp ω
∑
k ∈d ( ωp )
∑
k ∈d ( ωp )
p : ηp p
ω
xp p ≥ 0
ω
donde θp +p 1 ∈ ℝ , a(ωp ) indica el ascendiente del subproblema y d(ωp ) indica sus
descendientes.
114 04/11/2012
I.9 Mejoras en las técnicas de descomposición

Los algoritmos explicados hasta el momento son los algoritmos matemáticos
básicos. Pero en los métodos de descomposición es importante el uso de técnicas
que permitan reducciones significativas de tiempo de cálculo. De hecho las
mismas técnicas de descomposición están orientadas a resolver más eficazmente
problemas de muy gran tamaño. Estas mejoras se pueden agrupar en cuatro
tipos:
1. Método de optimización de subproblemas

Como los subproblemas se resuelven numerosas veces es conveniente utilizar
el método de optimización (y el optimizador) más adecuado. En principio,
para el método de descomposición anidada estocástica de Benders el método
simplex dual parece el más indicado puesto que en cada iteración como
maestro se hace infactible al añadir un corte mientras que como subproblema
cambia la cota de las restricciones. Sin embargo, siempre que sea posible es
necesario no sólo comprobar otros métodos alternativos como el simplex
primal o punto interior sino diferenciar su uso para la primera iteración o
para las sucesivas.
También se deben utilizar como punto de partida las bases de soluciones
previas de los subproblemas (estrategia denominada trickling o bunching). El
uso de un método u otro dependerá de la interfaz con el optimizador que a
su vez puede estar condicionada por el lenguaje utilizado en su codificación.
Por ejemplo, en el caso del lenguaje GAMS se deben tener en cuenta los
tiempos de creación y de interfaz además del tiempo propio de optimización.
Antes de resolver el problema estocástico descompuesto es conveniente la
resolución del problema determinista con el valor esperado de ω u otro valor
determinista cualquiera. Esto permite disponer de una información doble:
una solución inicial para cada etapa y valores iniciales para las todas las
variables.
2. Protocolos de barrido
El protocolo de barrido de los subproblemas define el orden en que éstos son
resueltos. Se han desarrollado diversas estrategias [Abrahamson, Birge:80,
Gassmann:90, Morton:93, Jacobs, Morton:96, Scott, Wittrock] que se definen
a continuación aunque los mejores resultados se han obtenido generalmente
con la estrategia denominada pasada rápida (fast-pass) [Morton:93], que es
por consiguiente la generalmente utilizada. Este proceso acaba cuando se
consigue convergencia en la primera etapa, que las cotas inferior y superior
04/11/2012 115
9 MEJORAS EN LAS TÉCNICAS DE DESCOMPOSICIÓN
de dicha etapa difieren menos de una cierta tolerancia, independientemente

del protocolo de barrido.
Para definir las estrategias se definen previamente estas operaciones:
• Avanzar (u, v )
Se entiende por avanzar entre u y v a formar las cotas de las
restricciones para los problemas de la etapa u , resolverlos, formar las
cotas para los problemas de la etapa u + 1 , resolverlos, etc. Así hasta
alcanzar la etapa v y resolver todos los problemas.
• Retroceder (u, v )
Se entiende por retroceder obtener los cortes para la etapa v , solucionar
los problemas de esta etapa, pasar los cortes a la etapa anterior,
solucionar los problemas, etc. Así hasta alcanzar la etapa u y formar sus
cortes sin solucionar los problemas de esta etapa.
Las estrategias de barrido se definen conceptualmente como:
• pasada rápida
Avanzar desde 1 hasta P , retroceder desde P − 1 hasta 1. Iterar entre
ambos procesos hasta convergencia.
• remolona
Empieza avanzando desde 1 hasta P . Luego selecciona la etapa p con
mayor error entre cota superior e inferior y para ésta resuelve todos los
problemas de las etapas p + 1, …, P antes de pasar los cortes a la etapa
p − 1 . Cuando la etapa con mayor error es la primera comienza otra vez
a iterar. Esta estrategia resuelve más problemas en las últimas etapas.
• cautelosa
Esta empieza con una iteración de la estrategia pasada rápida para
obtener cotas superior e inferior. Para una etapa dada p nunca avanza a
menos que el error en convergencia en la etapa p − 1 sea suficientemente
pequeño. Como resultado esta estrategia resuelve más problemas en las
primeras etapas.
3. Estrategias de formulación y agregación de cortes

Cuando existe más de un subproblema descendiente de uno dado (PLE-2 o
PLE-P) los cortes de Benders se pueden formular mediante el método mono
o multicorte. El primero tiene menor número de variables y restricciones
pero probablemente la convergencia será más lenta ya que la información
dual correspondiente a los diferentes subproblemas descendientes es agregada
dentro de un único corte.
116 04/11/2012
La formulación habitual del corte de Benders es la denominada linealización

alrededor de un punto ya que sólo requiere el uso y almacenamiento de las
variables duales correspondientes a las restricciones donde aparecen las
variables de la etapa anterior.
4. Estrategias de agregación de subproblemas (multicoordinación)

Los nodos del árbol de probabilidad que representa la incertidumbre se
pueden agregar arbitrariamente formando subárboles conexos. Las ventajas
de resolver un árbol en subárboles pueden ser de diverso tipo. La principal es
que subdividiendo el árbol en el menor número de subárboles posible18 se
disminuye el número de iteraciones del algoritmo de descomposición. En el
caso de uso de GAMS se disminuyen tiempos accesorios como son el de
creación e interfaz con el optimizador y no se penaliza demasiado el tiempo
de optimización cuando se utiliza un método de punto interior. En el caso de
haber utilizado procesamiento distribuido [Consigli, Dempster, Meyer] se
disminuye el tiempo de paso de información entre procesadores.
Mediante un procedimiento automático se agregan los nodos del árbol en
subárboles hasta alcanzar un determinado tamaño, por ejemplo, de
restricciones.
Se han estudiado y probado diversas estrategias de agregación de nodos de
un árbol de probabilidad para formar subárboles de un tamaño máximo
preespecificado tanto en sentido ascendente (de las hojas a la raíz) como
descendente (de la raíz a las hojas) [Cerisola]. La estrategia que consume
menor tiempo total es la de agregación ascendente de nodos de uno en uno
sin necesidad de agregar todos los nodos de una etapa pertenecientes a un
mismo padre.
En un problema de optimización estocástica resuelto por el método de

descomposición anidada estocástica de Benders el tiempo dependerá del tamaño
del árbol de escenarios. Si se supone que se desea resolver el problema en 1 hora
y se necesitan aproximadamente 6 iteraciones y cada subproblema tarda por
término medio 3 segundos en ser resuelto (en el caso de un problema LP de
20000 x 25000 x 80000 en un ordenador personal con procesador a 1 GHz), el
tamaño del árbol resoluble será de 200 nodos. El problema determinista
equivalente tendría 4000000 x 5000000 x 16000000.
18
Este número mínimo estará asociado a la memoria principal disponible para la resolución
del subárbol de tamaño máximo.
04/11/2012 117
I.10 Simulación en optimización lineal estocástica

bietapa
Cuando el espacio de estados de ω es muy elevado es necesario recurrir a
simulación o bien al formar los cortes en el método monocorte
min cT1 x 1 + θ2
x1 ,θ2
A1x 1 = b1
(9.1)
θ2 − ∑ p ω f2ωl ≥ ∑p ω
π2ωlT B1ω (x 1l − x 1 ) l = 1, …, j
ω∈Ω ω∈Ω
x1 ≥ 0
o bien para calcular la función de recursos en el método multicorte.

min
ω
cT1 x 1 + ∑ p ω θ2ω
x1 ,θ2
ω ∈Ω
A1x 1 = b1
(9.2)
θ2ω − f2ωl ≥ π2ωlT B1ω (x 1l − x 1 ) ω ∈Ω l = 1, …, j
x1 ≥ 0
Estas ecuaciones se puede interpretar, además de como una enumeración del

conjunto de estados ω ∈ Ω , como una integral o sumatorio en el hiperespacio
definido por los parámetros aleatorios del sistema.
La determinación de los valores x 2ω para cada ω es computacionalmente
muy costosa, implica la solución de un problema de programación lineal. Por
consiguiente, se deseará muestrear pocos estados. Cada estado corresponde a
una combinación de valores de los parámetros aleatorios19 del sistema.
z 2ω = c2ωT x 2ω = f2ω (e1, …, eE ) (9.3)

siendo E el número de parámetros del sistema y f2ω es una variable aleatoria
que se distribuye según la función de cuantía conjunta p ω .
Si estos parámetros son independientes entre sí la función de cuantía
conjunta se puede expresar como producto de las funciones de cuantía de los
parámetros. Entonces,
p ω = p1(e1 )⋯ pE (eE ) (9.4)

siendo p1(e1 ) = prob(parámetro 1 = e1 ) .
19
Por ejemplo, en el caso de sistemas de energía eléctrica éstos pueden ser los elementos de
generación y red, los escenarios de hidraulicidad y los niveles de demanda.
04/11/2012 119
10 SIMULACIÓN EN OPTIMIZACIÓN LINEAL ESTOCÁSTICA BIETAPA
El muestreo se puede realizar de dos formas diferentes:
• Muestreo externo
Se toman las muestras para reducir el tamaño del problema y luego se aplica
el método de optimización correspondiente para resolverlo.
• Muestreo interno
Se toman las muestras al mismo tiempo que se aplica el método de
optimización que lo resuelve.
Como primera aproximación a la esperanza de la función objetivo para el

conjunto de estados se puede realizar un muestreo de Monte Carlo tomando una
muestra de n valores, obtenidos de acuerdo a la función de cuantía conjunta
p ω . La media muestral es un estimador centrado de la esperanza del conjunto
de estados y su varianza es igual a la varianza muestral dividida por el número
de muestras. El intervalo de confianza de la media es función de estos valores
[Law].
El número de muestras a realizar será el necesario para conseguir un
intervalo de confianza de la media de la función objetivo de la segunda etapa
suficientemente pequeño:
1 n ωT ω 1 n
z2 = ∑ c2 x 2 = ∑ f2ω (9.5)
n ω =1 n ω =1
Una vez alcanzado este número de muestras, se forman los cortes muestrales.
En el caso de utilizar el método monocorte o híbrido (convenientes cuando el
número de muestras es elevado), éstos dejan de ser exactos y se convierten en
aproximaciones de los cortes reales:
1 n ωlT ω 1 n ωlT ω
∑ π2 B1 x1 + θ2 ≥ n ∑ π2 b2 (9.6)
n ω =1 ω =1
o bien
1 n ωl 1 n ωlT ω l
θ2 − ∑ f2 ≥ ∑ π2 B1 (x 1 − x 1 ) (9.7)
n ω =1 n ω =1
Los cortes calculados por simulación dejan de ser necesariamente planos
soporte de la función objetivo convexa, pueden intersecar la función de recursos
de los subproblemas de la segunda etapa. En [Infanger] se hacen los siguientes
supuestos para su obtención:
• El error en los cortes reside únicamente en las cotas, no en los gradientes, y

su varianza es la misma que la de la función objetivo.
120 04/11/2012
• Se conserva la misma base en el problema maestro para todas las soluciones

independientemente de los valores de las cotas de los cortes.
• Los cortes generados en diferentes iteraciones son independientes y son
obtenidos a partir de diferentes muestras.
La utilización de simulación dentro del proceso de descomposición afecta

tanto al número de muestras requeridas como al criterio de convergencia del
algoritmo. Ahora tanto la cota superior como inferior son variables aleatorias
obtenidas por muestreo.
Como primera aproximación se realiza la optimización con un número de
muestras fijado de antemano. Si se utiliza la descomposición de Bd, se denomina
menor cota superior a la cota superior con menor media en todas las iteraciones.
La mayor cota inferior será la cota inferior de la última iteración. La varianza
de cada cota se obtiene a partir de la varianza de la función objetivo de los
subproblemas. Desde un punto de vista práctico se utiliza el siguiente criterio de
detención [Infanger, Law]: el intervalo de confianza de la diferencia entre las
cotas superior e inferior ha de contener al cero. Una vez detenido el algoritmo
de descomposición, el intervalo de confianza de la solución óptima se calcula a
partir del límite inferior del intervalo de confianza de la cota inferior y del límite
superior del intervalo de confianza de la cota superior. A su vez este intervalo
ha de ser inferior a un valor especificado. En caso contrario es necesario
aumentar el número inicial de muestras y resolver de nuevo el problema. En
[Infanger:94, Morton:93, Morton:94b] se encuentran reglas de detención que
permiten derivar el número de muestras necesarias en cada iteración para
garantizar que la probabilidad de detención incorrecta sea inferior a un umbral.
La simulación tiene la ventaja de ser básicamente independiente del número
de parámetros que intervienen pero la relación de convergencia
aproximadamente es 1 n . Esto es, se necesita multiplicar por cuatro el número
de muestras para dividir por dos el tamaño del intervalo de confianza. Cuando
los sucesos con poca probabilidad de ocurrencia son los que ocasionan los valores
de la función de recursos más elevados la varianza de la media será grande y
utilizando simulación de Monte Carlo serán necesarias muchas muestras.
Entonces, es necesario recurrir a técnicas de reducción de varianza.
I.10.1 Técnicas de reducción de varianza

Las técnicas de reducción de varianza (variance reduction techniques VRT)
permiten reducir el tamaño del intervalo de confianza de una media de una
variable que interesa estimar sin perturbar el valor de ésta para un mismo
número de muestras (mejorar la precisión de la estimación) o, alternativamente,
04/11/2012 121
conseguir la precisión deseada con menor esfuerzo de muestreo. Los libros de

referencia para estudiar estas técnicas son [Bratley, Law, Morgan, Rubinstein].
Para utilizar las VRT es conveniente tener presente los siguientes criterios
generales:
• Habitualmente, es imposible conocer de antemano cuál va a ser la reducción

de varianza conseguida o incluso si ésta se va a producir. Las VRT deben
experimentarse considerando el sistema real a analizar.
• Para hacer un uso adecuado de las VRT se debe conocer en detalle el modelo
que reproduce el comportamiento del sistema.
• Como se verá más adelante, la utilización de VRT se puede entender como
una forma de aprovechar información sobre el sistema en cuestión.
• El uso de VRT para conseguir reducción de la varianza implica un
sobrecoste computacional para realizar ciertos muestreos preliminares o
cálculos complementarios al proceso mismo de simulación.
Entre las más utilizadas cabe mencionar:
• Números aleatorios comunes o muestreo correlacionado o simulación

comparativa o cadenas (parejas) sincronizadas
Se realizan muestreos para configuraciones diferentes del sistema con el
mismo conjunto de números aleatorios siendo utilizado cada número
aleatorio para la misma función en los diferentes muestreos.
• Variables antitéticas
Se basa en la idea de introducir correlación negativa entre dos muestras
consecutivas. Consiste en la utilización de números aleatorios
complementarios en dos muestras sucesivas.
• Variable de control
La idea básica es utilizar los resultados de un modelo más sencillo para
predecir o explicar parte de la varianza del valor a estimar. Se necesita un
cálculo previo del valor esperado de la variable de control. Este cálculo debe
ser muy rápido frente al de la variable a estimar.
• Muestreo por importancia

Se reemplaza la variable aleatoria a estimar por otra con la misma media
pero menor varianza. Se modifica la función de densidad utilizada en el
muestreo de manera que esté centrada en la zona de interés. Se evita el
muestrear sucesos probables pero no interesantes.
122 04/11/2012
• Muestreo estratificado
La idea intuitiva de esta técnica es similar a la anterior pero en versión
discreta. Consiste en tomar más muestras de la variable aleatoria en las
zonas de mayor interés. La varianza se reduce al concentrar el esfuerzo de
simulación en los estratos más relevantes.
A continuación se describe cada una de ellas en más detalle y se dan los

fundamentos matemáticos que las sustentan.
I.10.1.1 Números aleatorios comunes
La idea básica de esta técnica es la comparación de alternativas bajo las
mismas condiciones experimentales y, por consiguiente, las diferencias en los
atributos del sistema se deben a las diversas configuraciones y no a las
condiciones experimentales. Es la única técnica que se emplea para comparar
dos o más configuraciones de un sistema (por ejemplo, instalación de nuevas
líneas en la red de transporte de energía eléctrica).
Se consideran dos configuraciones alternativas. Sean X1j y X 2 j las muestras
independientes j para las dos configuraciones. Se desea estimar
ξ = µ1 − µ2 = E (X1 j ) − E (X 2 j ) . Se efectúan n muestras y se calcula
Z j = X 1 j − X 2 j , entonces E (Z j ) = ξ y la media muestral Z (n ) es un estimador
centrado de ξ .
n
Z (n ) =
∑ j =1
Zj
(9.8)
n
La varianza muestral se calcula como
var(Z j ) var(X1 j ) + var(X 2 j ) − 2 cov(X1 j , X 2 j )
var(Z (n )) = = (9.9)
n n
Si las muestras para las dos configuraciones se hicieran independientemente
(es decir, con números aleatorios diferentes) X1j y X 2 j serían independientes y
la covarianza sería 0. Por el contrario, si las muestras están positivamente
correlacionadas la covarianza será positiva y la varianza de la media muestral
será menor. Esta técnica trata de inducir una correlación positiva entre los
muestreos hechos para las dos configuraciones utilizando los mismos números
aleatorios. Esto es posible por la característica de reproducibilidad de las
cadenas de números pseudoaleatorios. Dicho de otra forma, el uso de una cadena
cuya semilla inicial es incontrolablemente aleatoria, en general, impide el uso de
esta técnica.
04/11/2012 123
No existe garantía de que esta técnica producirá una reducción de la

varianza de la media del atributo a estimar. Incluso si así lo hiciera no se puede
conocer de antemano la magnitud de la reducción. Su eficacia depende
completamente de los sistemas a comparar. Si fuera posible sería conveniente la
realización de un estudio piloto para probar su eficacia, de la siguiente manera.
Se efectúan n muestras para cada configuración utilizando los mismos números
aleatorios. Se calculan las varianzas muestrales del atributo S12 (n ) y S 22 (n ) , se
calcula la varianza muestral de las diferencias SZ2 (n ) . Si el uso de esta técnica es
fructífero se observará que SZ2 (n ) < S12 (n ) + S 22 (n ) y la diferencia entre ambos
lados es una estimación de la reducción de la varianza que se puede esperar.
La implantación correcta de esta técnica requiere la sincronización de los
números aleatorios para las diferentes configuraciones. Es decir, un número
aleatorio utilizado con un fin en una configuración deberá ser utilizado
exactamente para el mismo fin en otras configuraciones. No es suficiente el
comenzar ambos muestreos con la misma semilla. Una posibilidad puede ser
desperdiciar números aleatorios en ciertos puntos. Otra posibilidad es utilizar
una cadena de números aleatorios para cada parámetro aleatorio teniendo
especial cuidado en evitar los solapes entre cadenas de números. Si la dificultad
para efectuar la sincronización es elevada se puede pensar en sincronizar para
las diferentes configuraciones únicamente un subconjunto de los parámetros
aleatorios y generar independientemente las muestras para los otros.
I.10.1.2 Variables antitéticas
Como la técnica anterior trata de inducir correlación entre muestreos
separados, pero en este caso negativa. La idea básica es realizar muestreos por
parejas de manera que una observación inferior a la media tienda a ser corregida
por otra superior a la media en la siguiente observación. Si se toma la media de
las dos observaciones como una muestra, se tenderá a estar más cercano a la
media que se quiere estimar que si se tomaran individualmente las dos
observaciones.
La implantación de esta técnica se hace mediante el uso de números
aleatorios complementarios. Es decir, si un número se utiliza para un fin en una
observación se toma el complementario para el mismo fin en la siguiente
observación, están sincronizados.
Se supone que se efectúan n parejas de muestras siendo X j(1) la muestra j
correspondiente a la primera observación y X j(2) la correspondiente a la segunda
observación con número aleatorio complementario al anterior, en total 2n
observaciones. Ambas observaciones son válidas y E (X j(1) ) = E (X j(2) ) = µ . Se
124 04/11/2012
X j(1) + X j(2)
define Z j = y la media muestral Z (n ) como estimador centrado de
2
µ . La varianza muestral se calcula como
var(Z j ) var(X j(1) ) + var(X j(2) ) + 2 cov(X j(1), X j(2) )

var(Z (n )) = = (9.10)
n 4n
Si los muestreos dentro de cada pareja fueran hechos independientemente la
covarianza sería 0. Si se puede inducir correlación negativa entre muestreos la
covarianza será negativa y esto reducirá la varianza muestral.
Igual que con la técnica anterior, no existe garantía de que se producirá una
reducción de la varianza de la media. Incluso si así lo hiciera no se puede
conocer de antemano la magnitud de la reducción. Su eficacia depende del
sistema a simular. El requisito fundamental, aunque no suficiente, para que
pueda ser eficaz es que la respuesta del sistema sea monótona frente a la
variación del número aleatorio. Es decir, que si el número aleatorio crece el
atributo aumente (o disminuya) siempre. La combinación de esta técnica con la
anterior no tiene por qué ser eficaz.
I.10.1.3 Variable de control
Sea una variable aleatoria X cuya media desea estimarse E (X ) = µ , es la
variable controlada. Se considera otra variable Y , correlacionada con X
positiva o negativamente, cuyo valor esperado E (Y ) = ν es conocido, es la
variable de control.
La idea básica es similar a la de las anteriores técnicas, aprovechar la
correlación entre la variable de control y la controlada para obtener una
reducción de la varianza. La variable de control explica cierta parte de la
varianza de la variable controlada. Implícitamente se supone que la evaluación
de una muestra de la variable de control cuesta mucho menos tiempo que la
evaluación de la variable controlada.
Se define la variable aleatoria Z = X − a(Y − ν ) . Tanto Z como X tienen
la misma media E (Z ) = E (X ) = µ . Por consiguiente, Z es un estimador
centrado de µ que es el valor a estimar. La constante a tiene el mismo signo
que la correlación entre X e Y , se utiliza para amplificar la desviación (Y − ν )
para ajustar X . La varianza de Z se calcula como
var(Z ) = var(X ) + a 2 var(Y ) − 2a cov(X ,Y ) (9.11)

De manera que si 2a cov(X ,Y ) > a 2 var(Y ) entonces la varianza de Z será
menor que la varianza de X . Esta condición depende del valor de a y de la
variable de control Y . El valor de a que minimiza la varianza de Z para una
variable Y dada se determina derivando con respecto a a e igualando a 0.
04/11/2012 125
∂ var(Z )
= 2a var(Y ) − 2 cov(X ,Y ) = 0 (9.12)
∂a
cov(X ,Y )
Por lo tanto, a * = . De esta expresión se deduce que si la
var(Y )
correlación entre X e Y es muy fuerte, covarianza elevada, a * también tendrá
un valor alto. Por consiguiente, producirá un ajuste fuerte a la desviación de Y
sobre su valor medio ν que nos da una indicación de la desviación de X sobre
su valor medio µ . Por otra parte, si la varianza de Y es pequeña también se
obtiene un valor elevado de a * ya que se tiene una mayor precisión en el valor
observado de Y .
Para este valor óptimo a * el valor de la varianza de Z es
cov2 (X ,Y )
var(Z * ) = var(X ) − 2
= (1 − ρXY ) var(X ) (9.13)
var(Y )
donde ρXY es el coeficiente de correlación entre X e Y . En caso de correlación
perfecta ρXY = ±1 entre ambas variables la varianza de Z es nula.
Dado que es imposible obtener el valor exacto de a * un método simple es
reemplazar a * por su estimador obtenido de n muestras
ˆ XY (n )
cov
aˆ* (n ) = (9.14)
SY2 (n )
n
∑ X j − X (n ) Yj −Y (n )
ˆ XY
siendo cov (n ) = 
j =1  .
n −1
Conviene mencionar que también se pueden utilizar varias variables de
control simultáneamente con correlación no sólo con la variable controlada sino
incluso entre ellas.
A la vista de las ecuaciones anteriores se puede decir que una buena variable
de control debe estar fuertemente correlacionada con la variable controlada para
proporcionar mucha información sobre ésta y hacer un buen ajuste. Además es
deseable que la variable de control tenga poca varianza. Para encontrar una
variable de control se puede hacer un análisis de la estructura del modelo o bien
recurrir a experimentación. Las variables de control pueden clasificarse en:
• Internas
Son parámetros del modelo o funciones de los mismos, como la media. Sus
valores medios son habitualmente conocidos. Estas variables de control
deben ser generadas en todo caso, luego no añaden esfuerzo de simulación y
por esta razón es importante probar aunque no reduzcan significativamente
la varianza.
126 04/11/2012
• Externas
Son aquellas variables aleatorias de salida derivadas de introducir
simplificaciones en el modelo que lo hacen inadecuado como tal pero
permiten utilizarlas como variable de control. En este caso se necesita un
muestreo adicional con números aleatorios comunes para la variable de
control. La correlación entre variable de control y controlada deberá ser
ahora mayor para que merezca la pena el esfuerzo adicional de simulación.
Un ejemplo de variables de control externas son las medidas de fiabilidad
sólo de generación o sólo de red para el cálculo de medidas de fiabilidad
compuesta (debida a fallos de generación y/o de red).
I.10.1.4 Muestreo por importancia
La idea intuitiva detrás del muestreo por importancia es aumentar la
frecuencia de aparición en el muestreo de los sucesos que tienen más peso en la
evaluación de la función. El método consiste en obtener muestras de una nueva
función cuya función de cuantía es diferente y deformada con respecto a la
original [Bratley, Dantzig:90]. Ambas funciones tienen la misma media pero la
varianza de la segunda es muy inferior. Aplicándolo a la función objetivo de los
subproblemas se tiene:
fω
zˆ2 = E ωc2ωT x 2ω = ∑ p ω f ω = ∑ p ωq ω = ∑ p ωq ω g ω = z 2′ (9.15)
ω∈Ω ω∈Ω qω ω ∈Ω
donde E ω indica la esperanza matemática con respecto a ω y g ω es una

variable aleatoria que se distribuye según la función de cuantía conjunta p ωq ω .
La varianza poblacional de la función original es:
2 2
var(z 2 ) = ∑ p ω ( f ω ) − (zˆ2 ) (9.16)
ω∈Ω
La varianza poblacional de la nueva función será:

2
ω 2 2 (f ω ) 2
var(z 2′ ) = ∑ p q (g ) − (zˆ2′ ) = ∑ p
ω ω ω
− (zˆ2 ) (9.17)
ω ∈Ω ω ∈Ω q ω
Veamos una aplicación concreta de muestreo por importancia [Infanger]. Se

define un punto base y se calcula el incremento del valor de la función en la
coordenada i cuando el parámetro i pasa del valor base ei0 a eii .
∆i (ei ) = f (e10, …, eii , …, eE0 ) − f (e10, …, ei0 , …, eE0 ) (9.18)

Los incrementos en cada coordenada serán positivos si se toma como punto
base aquél que proporcione el menor valor de la función objetivo (siendo ésta,
por ejemplo, los costes variables de explotación). En el caso de sistemas de
energía eléctrica éste sería aquél con todos los elementos de generación
04/11/2012 127
disponibles, mayor valor de hidraulicidad y menor valor de demanda, siendo

éstos los parámetros aleatorios.
Se llama ∆i al incremento medio del valor de la función objetivo en la
coordenada i al tomar el parámetro i sus diferentes valores.
eii
∆i = ∑ p (e )∆ (e )i i i i (9.19)
ei =ei0
La esperanza de la función objetivo de los subproblemas es:

zˆ2 = ∑ p ω f ω =
ω ∈Ω
= f0 + ∑ p ω ∆f ω
ω ∈Ω
E
∆f ω
= f0 + ∑ p ω E ∑∆ ω
i (ei ) (9.20)
ω ∈Ω ∑ i =1
∆ (e )
ω
i i i =1
E
= f0 + ∑ p ω g ω ∑ ∆iω (ei )
ω ∈Ω i =1
E E
= f0 + ∑ g ω ∏ pi (ei )∑ ∆iω (ei )
ω ∈Ω i =1 i =1
∆f ω
donde g ω = E será una función igual a uno o superior si hay efectos
∑
i =1
∆i
ω
(ei )
de sinergia positiva con la variación de los parámetros. El cálculo de la
esperanza de la función objetivo de los subproblemas es una integral
multidimensional en el hiperespacio de los parámetros del sistema.
Manipulando y desarrollando la ecuación anterior se obtiene:
E
∆iω (ei )pi (ei )
z 2′ = f0 + ∑ ∆i ∑ g ω ∏ pk (ek ) (9.21)
i =1 ω∈Ω ∆i k ≠i
Ahora la esperanza se convierte en tantas integrales unidimensionales de sus

incrementos marginales como parámetros hay.
De acuerdo con esta ecuación, la función de cuantía utilizada para el
muestreo de la función g ω en cada dimensión i es la original pk (ek ) para todos
los parámetros excepto el i y para éste la función de cuantía de importancia
∆iω (ei )pi (ei )
. Ésta se obtiene al ponderar la función de cuantía original por el
∆i
impacto relativo en la función objetivo que tiene la variación del parámetro i
hasta el valor eii con respecto al caso base. Esta función de cuantía conjunta es
también el producto de las funciones de cuantía de los parámetros.
En resumen, se ha pasado de una integral multidimensional a un sumatorio
de integrales unidimensionales y además, se muestrea en función del impacto en
128 04/11/2012
la función objetivo con respecto al caso base, deformando la función de cuantía

original. La expresión anterior está compuesta por un término constante f0 y un
sumatorio de E esperanzas matemáticas de la función g ω . Para su estimación
se utiliza otra técnica de reducción de varianza, el muestreo estratificado, que se
comenta en la siguiente sección.
Con este muestreo por importancia se obtiene mayor reducción en la
varianza cuantas más dimensiones (parámetros aleatorios) tenga el espacio
muestral. No se necesita muestrear en las coordenadas cuyos incrementos
medios son nulos ya que no afectan la media. Este caso aparece cuando el
sistema está sobredimensionado.
Por supuesto, se podría aplicar de nuevo el muestro por importancia a la
función g ω calculando de nuevo los incrementos en cada coordenada y demás.
Pero este paso no merece la pena si el número de coordenadas es elevado.
La varianza poblacional de la función objetivo de los subproblemas es:
2
var(z 2 ) = ∑ p ω ( f ω ) − zˆ22
ω ∈Ω
2
= ∑ p ( f0 + ∆f ω ) − zˆ22
ω
ω ∈Ω
2
= f02 + 2 f0 (zˆ2 − f0 ) − zˆ22 + ∑ p ω (∆f ω )
ω ∈Ω
2 2
= − (zˆ2 − f0 ) + ∑ p ω (∆f ω ) (9.22)
ω∈Ω
2
2 (∆f ω ) E
= − (zˆ2 − f0 ) + ∑ p ω E ∑ ∆ (e ) ω
i i
ω∈Ω ∑ i =1
∆ (e )
ω
i i i =1
E
2
= − (zˆ2 − f0 ) + ∑ p ωh ω ∑ ∆iω (ei )
ω∈Ω i =1
E E
2
= − (zˆ2 − f0 ) + ∑ h ω ∏ pi (ei )∑ ∆iω (ei )
ω∈Ω i =1 i =1
2
(∆f ω )
donde h =
ω
E .
∑
i =1
∆i
ω
(ei )
Manipulando y desarrollando la ecuación anterior se obtiene:
E
2 ∆iω (ei )pi (ei )
var(z 2 ) = − (zˆ2 − f0 ) + ∑ ∆i ∑ h ω ∏ pk (ek ) (9.23)
i =1 ω ∈Ω ∆i k ≠i
Esta ecuación, que calcula la varianza de la distribución original, tiene una

estructura similar a la del cálculo de la media y, por consiguiente, permite su
implantación simultánea.
Se calcula ahora la varianza poblacional de la nueva función:
04/11/2012 129
E 2  2
ω 2 ∆i (ei )pi (ei )
ω
var(z 2′ ) = ∑ ( i ) ∑
∆  (g )
∆i
∏ pk (ek ) − (gˆi ) 

(9.24)
i =1  ω∈Ω k ≠i 
El cociente entre ambas nos da la reducción conseguida mediante el
muestreo por importancia para el cálculo de la media.
Para aclarar la obtención de la función de cuantía de importancia veamos un
ejemplo. Se supone que se trata de evaluar los costes de explotación de un
sistema eléctrico con tres generadores. Su potencia disponible es una variable
aleatoria con distribución multinomial. Es nula con una probabilidad q1 = 0.05 ,
igual a su mínimo técnico con probabilidad q 2 = 0.05 e igual a su potencia
nominal con probabilidad p = 0.90 . En la tabla se presenta el impacto marginal
en la función objetivo (costes de explotación) con respecto al caso base del fallo
de cada generador y el valor de la probabilidad de importancia obtenida a partir
de estos valores y de las probabilidades originales. Se toma como caso base
aquél con todos los generadores totalmente disponibles.
Estado del grupo Probabilidad original Impacto marginal Probabilidad importancia

pi (ei ) en función objetivo
∆i (ei ) ∑ p (e )∆ (e )
i i i i
∆i
Generador 1
Funcionamiento 0.90 0 0.0
Fallo parcial 0.05 30 0.3
Fallo total 0.05 70 0.7
Generador 2
Generador 3
Tabla 4.1 Funciones de cuantía de las variables aleatorias.
Veamos ahora un caso completo de aplicación del muestreo por importancia.

Se supone que se tiene un sistema compuesto de tres generadores, cada uno con
probabilidad de funcionamiento igual a 0.95 y de fallo 0.05. Para cada estado
del sistema se representa en la siguiente tabla su probabilidad, el valor de los
costes de generación del sistema y los incrementos con respecto al caso base. Se
toma como case base aquél con todos los grupos disponibles. Los incrementos
medios en cada coordenada ∆i son respectivamente 1.25 [=0.05(225–200)], 5 y
2.5. Obsérvese que cuando la variable aleatoria tiene sólo dos estados, la función
de cuantía de importancia es 0 cuando funciona el generador y 1 cuando falla.
130 04/11/2012
Estado grupos Probabilidad Coste Incremento

pω fω ∆f ω
1 1 1 0.857375 200 0
1 1 0 0.045125 250 50
1 0 1 0.045125 300 100
1 0 0 0.002375 500 300
0 1 1 0.045125 225 25
0 1 0 0.002375 425 225
0 0 1 0.002375 350 150
0 0 0 0.000125 1075 875
Tabla 4.2 Estados del sistema.
Estado del grupo Probabilidad original Impacto marginal Probabilidad importancia

pi (ei ) en función objetivo
∆i (ei ) ∑ p (e )∆ (e )
i i i i
∆i
Generador 1
Generador 2
Generador 3
Tabla 4.3 Funciones de cuantía de las variables aleatorias.
La aplicación directa del método de Monte Carlo para un número

suficientemente elevado de muestras resultaría en la aparición de cada estado
con una frecuencia proporcional a su probabilidad. La media de los costes será:
zˆ2 = 0.857375 ⋅ 200 + 0.045125 ⋅ 250 + 0.045125 ⋅ 300
+0.002375 ⋅ 500 + 0.045125 ⋅ 225 + 0.002375 ⋅ 425
+0.002375 ⋅ 350 + 0.000125 ⋅ 1075
= 209.609375
Si ahora en lugar de Monte Carlo se utiliza el muestreo por importancia y de
nuevo tomando un número suficientemente elevado de muestras la frecuencia de
aparición de cada estado será proporcional a su probabilidad. La media se
calcula como:
04/11/2012 131
zˆ2′ = 200
 25 225 150 875 
+1.25  0.9025 + 0.0475 + 0.0475 + 0.0025 
 25 25 + 50 25 + 100 25 + 100 + 50 
 100 300 150 875 
+5  0.9025 + 0.0475 + 0.0475 + 0.0025 
 100 100 + 50 25 + 100 25 + 100 + 50 
 50 300 225 875 
+2.5  0.9025 + 0.0475 + 0.0475 + 0.0025 
 50 100 + 50 25 + 50 25 + 100 + 50 
= 200 + 1.25 ⋅ 1.1145 + 5 ⋅ 1.067 + 2.5 ⋅ 1.1525
= 209.609375
La varianza de la población original se calcula como:
var(z 2 ) = 0.857375 ⋅ 2002 + 0.045125 ⋅ 2502 + 0.045125 ⋅ 3002
+0.002375 ⋅ 5002 + 0.045125 ⋅ 2252 + 0.002375 ⋅ 4252
+0.002375 ⋅ 3502 + 0.000125 ⋅ 10752 - 209.6093752
= 983.0505371
Si ahora se aplica la ecuación (9.22) la varianza de la distribución será la
misma
var(z 2 ) = −(209.609375 − 200)2
 252 2252 1502 8752 
+1.25 0.9025
 + 0.0475 + 0.0475 + 0.0025 
 25 25 + 50 25 + 100 25 + 100 + 50 
 1002 3002 1502 8752 
+5  0.9025 + 0.0475 + 0.0475 + 0.0025 
 100 100 + 50 25 + 100 25 + 100 + 50 
 502 3002 2252 8752 
+2.5  0.9025 + 0.0475 + 0.0475 + 0.0025 
 50 100 + 50 25 + 50 25 + 100 + 50 
= 983.0505371
La varianza de la nueva función será:
132 04/11/2012
 2 2
 25   225 
var(z 2′ ) = 1.252 0.9025   + 0.0475  


 25   25 + 50 
2 2
 150   875  
2
+0.0475   + 0.0025 
  − 1.1145
 25 + 100   25 + 100 + 50  

 2 2
100   300 
+52  0.9025   + 0.0475  


100  100 + 50 
2 2 
 150   875  2
+0.0475   + 0.0025 
  − 1.067
 25 + 100   25 + 100 + 50  

2
  50 
2
 300 
+2.52  0.9025   + 0.0475  


 50  100 + 50 
2 2
 225   875  
2
+0.0475   + 0.0025 
  − 1.1525
 25 + 50   25 + 100 + 50  

= 4.05365742
Esto es, con el muestreo por importancia se ha conseguido en este ejemplo
una reducción de la varianza de dos órdenes de magnitud.
I.10.1.5 Muestreo estratificado
Los ejemplos clásicos de estratificación aparecen en el diseño de encuestas. Se
supone que se desea efectuar una encuesta a la población de una ciudad y hacer
un análisis estadístico sobre las respuestas. El método de Monte Carlo elegiría
una persona cualesquiera de forma inversamente proporcional al número de
habitantes. El método de estratificación dividiría la población en estratos, por
ejemplo según población de los distritos, y repartiría el número total de
encuestas de forma proporcional al tamaño de cada estrato. En este caso, dicho
tamaño se puede obtener con precisión a partir del censo. De esta manera se
reemplaza la variable aleatoria número de veces que dicho estrato se muestrea
por su valor esperado. Esto necesariamente reduce la varianza.
En la técnica de muestreo estratificado se divide el espacio muestral en
estratos disjuntos. La idea intuitiva de esta técnica es similar a la de muestreo
por importancia pero en versión discreta. Consiste en tomar más muestras de la
variable aleatoria en los estratos de mayor interés. La varianza se reduce al
concentrar el esfuerzo de simulación en los estratos más relevantes.
Por ejemplo, en el caso de un modelo de explotación que divide el alcance
del estudio en periodos de diferente duración. El tamaño del estrato se toma
proporcional al valor de la función objetivo (habitualmente costes variables de
explotación) en cada periodo. Una aproximación es hacerlo proporcional a la
04/11/2012 133
duración del periodo suponiendo implícitamente que el coste variable de

explotación por unidad de tiempo será similar en los diferentes periodos.
Para obtener la mínima varianza en la estimación se reparte el número total
de muestras a realizar en estratos proporcionales al producto de la función
objetivo por la desviación típica de la variable a estimar. Como ésta no es
conocida se efectúa un estudio piloto que determine la desviación estándar de la
variable en cada estrato.
Incluso tomando un número de muestras en cada estrato proporcional a su
tamaño se garantiza que la varianza obtenida por muestreo estratificado es
menor o igual que la que se obtendría por Monte Carlo.
Aplicando este muestreo al problema de optimización estocástica, el número
de muestras para cada coordenada es proporcional a su impacto medio en la
función objetivo20. Dado un número total de muestras a realizar n a cada
parámetro i se le dedica un número ni , aproximadamente proporcional a su
E
valor ∆i , teniendo presente que ni ≥ 1 y ∑ i =1 ni = n .
La media muestral de la función objetivo de los subproblemas si se aplica
muestreo por importancia y estratificado será:
E
∆i ni ω
z 2 = f0 + ∑ ∑g (9.25)
i =1 ni ω=1
La varianza muestral de la función original
E
∆i ni ω
var(z 2 ) = −(z 2 − f0 )2 + ∑ ∑h (9.26)
i =1 ni − 1 ω =1
La varianza muestral de la nueva función
2
E
(∆ )
i  ni ω 2 
 ∑ (g ) − (gi )2 
var(z 2′ ) = ∑ (9.27)
i =1 ni − 1  ω =1 
Y la varianza de la media muestral de la nueva función

2
E
(∆ )
i  ni ω 2 
 ∑ (g ) − (gi )2 
var(z 2′ ) = ∑ (9.28)
i =1 ni (ni − 1)  ω =1 
El proceso completo de muestreo conjunto por importancia y estratificado

está compuesto de dos fases.
20
Esto presenta el inconveniente de atender únicamente a una función objetivo. Si hubiera
varias éstas podrían entrar en conflicto.
134 04/11/2012
E
• La primera es la de preparación. En ella se resuelven 1 + ∑ i =1 (ei − 1)
 
subproblemas para obtener el valor de la función objetivo en el caso base
más los incrementos ∆i (ei ) en cada coordenada. Con estos datos se evalúan
las funciones de cuantía de importancia.
• En la segunda fase se muestrea en combinación:

muestreo estratificado para cada coordenada (número de muestras
proporcional al incremento medio de la función objetivo en la
coordenada) y
muestreo por importancia de acuerdo a la función de cuantía conjunta
obtenida previamente (función de cuantía de importancia para esa
coordenada y función original para el resto).
Para cada muestra de los valores de los parámetros se resuelve un
subproblema y se evalúa la función g ω (y h ω ) como cociente entre el
valor de la función objetivo del subproblema menos el valor en el caso
base (al cuadrado en el caso de h ω ) y la suma de los valores
incrementales calculados a partir de los datos de la fase preparatoria. A
partir de los valores obtenidos en las muestras se evalúa la media y su
varianza de acuerdo con las ecuaciones anteriores y su intervalo de
confianza.
Como muestra de la efectividad en la reducción de la varianza de este

muestreo combinado para un modelo de planificación de la expansión de la
generación y de la red con 107 escenarios posibles sólo fueron necesarias 100
muestras (además de las de la fase preparatoria) para conseguir un intervalo de
confianza de la función objetivo de tamaño aproximado del 0.2 % con
aceptabilidad del 95 % [Infanger].
04/11/2012 135
I.11 Referencias
Aardal, K. and Ari, A. “On the Resemblance between the Kornai-Liptak and
Cross Decomposition Techniques for Block-Angular Linear Programs“.
European Journal of Operational Research Vol 46, No 3, pp 393-398. June
1990.
Abrahamson, P.G. “A Nested Decomposition Approach for Solving Staircase
Linear Programs” Systems Optimization Laboratory. Department of
Operations Research. Stanford University. SOL 83-4. June 1983.
Alonso-Ayuso, A., E. Cerdá, L. F. Escudero, R. Sala (eds.) (2004) Optimización
bajo incertidumbre Tirant lo Blanch. Valencia, España.
Beale, E.M.L. “On Minimizing a Convex Function Subject to Linear
Inequalities” Journal of the Royal Statistical Society. Vol 17b, pp 173-184.
1955.
Benders, J.F. “Partitioning Procedures for Solving Mixed-Variables
Programming Problems” Numerische Mathematik. Vol 4, pp 238-252. 1962.
Birge, J.R. “Solution Methods for Stochastic Dynamic Linear Programs”
Systems Optimization Laboratory. Department of Operations Research.
Stanford University. SOL 80-29. December 1980.
Birge, J.R. “Decomposition and Partitioning Methods for Multi-Stage
Stochastic Linear Programs” Operations Research. Vol 33. pp 989-1007.
1985.
Birge, J.R. and Louveaux, F.V. “A Multicut Algorithm for Two-Stage
Stochastic Linear Programs”. European Journal of Operational Research.
Vol 34, No 3, pp 384-392. March 1988.
Birge, J.R. and Louveaux, F.V. Introduction to Stochastic Programming
Springer-Verlag. New York, USA. 1997.
Bratley, P., Fox, B.L. and Schrage, L.E. A Guide to Simulation. Second
Edition. Springer-Verlag. New York, USA. 1987.
Brooke, A., Kendrick, D., Meeraus, A. and Raman, R. (1998) GAMS A User’s
Guide. GAMS Development Co.
Cerisola, S. and Ramos, A. “Node Aggregation in Stochastic Nested Benders
Decomposition Applied to Hydrothermal Coordination” 6th PMAPS.
September 2000.
Consigli, G. and Dempster, M.A.H. “Solving Dynamic Portfolio Problems Using
Stochastic Programming” (internal technical report).
04/11/2012 137
11 REFERENCIAS
Dantzig, G.B. “Linear Programming Under Uncertainty” Management Science.

Vol 1, No 3-4, pp 197-206. April-July 1955.
Dantzig, G.B. Linear Programming and Extensions. Princeton University Press.
Princeton, USA. 1963.
Dantzig, G.B. “Planning Under Uncertainty Using Parallel Computing”
Stanford University. SOL 87-1. January 1987.
Dantzig, G.B. and Glynn, P.W. “Parallel Processors for Planning Under
Uncertainty” Annals of Operations Research. Vol 22, pp 1-21. 1990.
Dantzig, G.B. and Infanger, G. “Large-Scale Stochastic Linear Programs:
Importance Sampling and Benders Decomposition” Systems Optimization
Laboratory. Department of Operations Research. Stanford University. SOL
91-4. March 1991.
Dantzig, G.B., Ho, J.K. and Infanger, G. “Solving Stochastic Linear Programs
on a Hypercube Multicomputer” Systems Optimization Laboratory.
Department of Operations Research. Stanford University. SOL 91-10.
August 1991.
Dantzig, G.B. and Infanger, G. “Multi-Stage Stochastic Linear Programs for
Portfolio Optimization” Systems Optimization Laboratory. Department of
Operations Research. Stanford University. SOL 91-11. September 1991.
Dantzig, G.B. and Madansky, A. “On the Solution of Two-Stage Linear
Programs Under Uncertainty” Proceedings of the Fourth Berkeley
Symposium on Mathematical Statistics and Probability. University of
California Press. Berkeley. Vol I, pp 165-176. 1961.
Dantzig, G.B. and Thapa, M.N. Linear Programming Introduction. Springer,
1997.
Dantzig, G.B. and Wolfe, P. “Decomposition Principle for Linear Programs”
Operations Research. Vol 8, No 1, pp 101-111. January-February 1960.
Dempster, M.A.H. and Thompson, R.T. “Paralelization and Aggregation of
Nested Benders Decomposition” (internal technical report).
Duran, M.A. and Grossmann, I.E. “An Outer-Approximation Algorithm for a
Class of Mixed-Integer Nonlinear Programs” Mathematical Programming.
Vol 36, pp 307-339. 1986.
Entriken, R. “Using MINOS as a Subroutine for Decomposition” Systems
Optimization Laboratory. Department of Operations Research. Stanford
University. SOL 87-7. June 1987.
138 04/11/2012
Entriken, R. “A Parallel Decomposition Algorithm for Staircase Linear

Programs” Systems Optimization Laboratory. Department of Operations
Research. Stanford University. SOL 88-21. December 1988.
Entriken, R. “The Parallel Decomposition of Linear Programs” Systems
University. SOL 89-17. November 1989.
Entriken, R. and Infanger, G. “Decomposition and Importance Sampling for
Stochastic Linear Models” IBM Research Report. RC 15545 (#69107).
March 1990.
Ermoliev, Y. and Wets, R.J-B. (eds.) Numerical Techniques for Stochastic
Optimization. Springer-Verlag. Berlin, Germany. 1988.
Flippo, O. E. and Rinnooy Kan, A. H. G. “Decomposition in General
Mathematical Programming” Mathematical Programming. Vol 60, pp 361-
382. 1993.
Floudas, C.A., Aggarwal, A. and Ciric, A.R. “Global Optimum Search for
Nonconvex NLP and MINLP Problems” Computer in Chemical
Engineering. Vol 13, No 10, pp 1117-1132. 1989.
Fushimi, M. (1989) “Random Number Generation on Parallel Processors”
Proceedings of the 1989 Winter Simulation Conference. pp 459-461.
Gassmann, H.I. “MSLiP: A Computer Code for the Multistage Stochastic
Linear Programming Problem” Mathematical Programming. Vol 47, pp 407-
423. August 1990.
Gassmann, H.I. and Ireland, A.M. “Scenario Formulation in an Algebraic
Modelling Language” Annals of Operations Research. Vol 59, pp 45-75.
1995.
Gassmann, H.I. and Ireland, A.M. “On the Formulation of Stochastic Linear
Programs using Algebraic Modelling Languages” Annals of Operations
Research. Vol 64, pp 83-112. 1996.
Geoffrion, A.M. “Elements of Large-Scale Mathematical Programming”
Management Science. Vol 16, No 11, pp 652-691. July 1970.
Geoffrion, A.M. “Large-Scale Linear and Nonlinear Programming” in
Optimization Methods for Large Scale Systems with Applications. McGraw-
Hill. New York, USA. 1971.
Geoffrion, A.M. “Generalized Benders Decomposition” Journal of Optimization
Theory and Applications. Vol 10, No 4, pp 237-260. October 1972.
Geoffrion, A.M. “Lagrangian Relaxation for Integer Programming”
Mathematical Programming Study. Vol 2, pp 82-114. 1974.
04/11/2012 139
11 REFERENCIAS
Gill, P. E., Murray, W. and Wright, M.H. Numerical Linear Algebra and
Optimization. Addison Wesley, 1991.
Hilger, J., Harhalakis, G. and Proth, J.M. “Generalized Cross-Decomposition
Method: Algorithm and Implementation” Institut National de Recherche en
Informatique et en Automatique. Report 1055. July 1989.
Infanger, G. “Monte Carlo (Importance) Sampling within a Benders
Decomposition Algorithm for Stochastic Linear Programs. Extended
Version: Including Results of Large-Scale Problems” Systems Optimization
Laboratory. Department of Operations Research. Stanford University. SOL
91-6. March 1991.
Infanger, G. Planning Under Uncertainty: Solving Large-Scale Stochastic Linear
Programs Boyd Fraser. 1994.
Jacobs, J, Freeman, G., Grygier, J., Morton, D. P., Schultz, G., Staschus, K.
and Stedinger, J. “SOCRATES: A System for Scheduling Hydroelectric
Generation under Uncertainty” Annals of Operations Research. Vol 59, pp
99-133. 1995.
Kall, P. and Wallace, S.W. Stochastic Programming. Wiley, 1994.
Kim, S., Cho, S.-C. and Um, B.-S. “A Simplified Cross Decomposition
Algorithm for Multiple Right Hand Choice Linear Programming” Journal of
Operations Research Society of Japan. Vol 32, No 4, pp 441-449. December
1989.
Krishna, A. “Enhanced Algorithms for Stochastic Programming”. PhD Thesis.
Stanford University. September 1993.
Lasdon, L.S. Optimization Theory for Large Systems. The MacMillan Company.
New York, USA. 1970.
Law, A.M. and Kelton, W.D. Simulation Modeling and Analysis. Second Edition
McGraw-Hill. New York, USA. 1991.
Massachusetts Institute of Technology “Electric Generation Expansion Analysis
System. Volume 1: Solution Techniques, Computing Methods, and Results.
Volume 2: Details of Solution Techniques, Data of Test Systems, and
Glossary of Terms” Electric Power Research Institute. EPRI EL-2561.
August 1982.
Mathiesen, L. “Pricing of Electricity in the Presence of Water Uncertainty”
February 1992.
The MathWorks, Inc. MATLAB User's Guide. August 1992.
The MathWorks, Inc. Optimization Toolbox for Use with MATLAB User's
Guide. August 1992.
140 04/11/2012
Mayer, J. (1998) Stochastic Linear Programming Algorithms. A Comparison

Based on a Model Management System. Gordon and Breach Science
Publishers.
Meyer, R.R. and Zakeri, G. “Multicoordination Methods for Solving Convex
Block-Angular Programs” SIAM Journal on Optimization. Vol 10, No 1, pp
121-131. 1999.
Moré, J. J. and Wright, S. Numerical Optimization Algorithms and Software.
Software Guide and Course Notes SIAM Tutorial. May 1992.
Morgan, B.J.T. Elements of Simulation. Chapman and Hall. London, UK. 1984.
Morton, D. P. “Algorithmic Advances in Stochastic Programming” Systems
University. SOL 93-6. July 1993.
Morton, D. P. “An Enhanced Decomposition Algorithm for Multistage
Stochastic Hydroelectric Scheduling” Annals of Operations Research. Vol
64, pp 211-235. 1996.
Morton, D. P. “Stopping Rules for a Class of Sampling-Based Stochastic
Programming Algorithms” Naval Postgraduate School. NPS-OR-94-003.
January 1994.
Murtagh, B.A. and Saunders, M.A. “MINOS 5.1 User's Guide” Systems
University. SOL 83-20R. December 1983. Revised January 1987.
Mulvey, J.M. “Introduction to financial optimization: Mathematical
Programming Special Issue” Mathematical Programming Vol 89 (2001) 2,
pp. 205-216
Paules, G.E. and Floudas, C.A. “APROS: Algorithmic Development
Methodology for Discrete-Continuous Optimization Problems” Operations
Research. Vol 37, No 6, pp 902-915. November-December 1989.
Pereira, M.V.F. and Pinto, L.M.V.G. “Multi-Stage Stochastic Optimization
Applied to Energy Planning” Mathematical Programming. Vol 52, pp 359-
375. August 1991.
Pérez Arriaga, J.I., Ramos, A. y Latorre, G. “Volumen 1: Guía de usuario del
programa de planificación estática de la red de transporte a largo plazo.
PERLA. Volumen 2: Casos ejemplo” Instituto de Investigación Tecnológica.
Julio 1991.
Prékopa, A. Stochastic Programming. Kluwer Academic Publishers, 1995.
Ramos, A. “Integrated Model for Capacity Planning for Manufacturing
Systems” March 1992.
04/11/2012 141
11 REFERENCIAS
Ramos, A. “Modelos de explotación a medio plazo de la generación eléctrica.

Aplicaciones para el sistema español” Tesis doctoral. Universidad
Politécnica de Madrid. Noviembre 1990.
Ramos, A., A. Alonso-Ayuso, G. Pérez (eds.) Optimización bajo Incertidumbre
Universidad Pontificia Comillas. Madrid, España 2008
Rockafellar, R.T. and Wets, R.J.B. “Scenarios and Policy Aggregation in
Optimization Under Uncertainty” Mathematics in Operations Research. Vol
16, pp 119-147.
Rubinstein, R. Y. Simulation and the Monte Carlo Method John Wiley Sons.
New York, USA. 1981.
Ruszczynski, A. “Parallel Decomposition of Multistage Stochastic Programming
Problems” Mathematical Programming. Vol 58, pp 201-228. 1993.
Schweitzer, E. “Modifying MINOS for Solving the Dual of a Linear Program”
Stanford University. SOL 90-11. August 1990.
Scott, D.M. “A Dynamic Programming Approach to Time-Staged Convex
Programs” Systems Optimization Laboratory. Department of Operations
Research. Stanford University. SOL 85-3. March 1985.
Stanford University “Decomposition Techniques for Multi-Area Generation and
Transmission Planning Under Uncertainty” Electric Power Research
Institute. EPRI EL-6484. August 1989.
Vanderbei, R.J. Linear Programming: Foundations and Extensions. Kluwer
Academic Publishers, 1996.
Van Roy, T.J. “Cross Decomposition for Mixed Integer Programming”
Mathematical Programming. Vol 25, pp 46-63. 1983.
Van Slyke, R.M. and Wets, R.J.-B. “L-Shaped Linear Programs with
Application to Optimal Control and Stochastic Programming” SIAM
Journal on Applied Mathematics. Vol 17, No 4, pp 638-663. July 1969.
Wittrock, R.J. “Advances in a Nested Decomposition Algorithm for Solving
Staircase Linear Programs” Systems Optimization Laboratory. Department
of Operations Research. Stanford University. SOL 83-2. March 1983.
142 04/11/2012
I.12 Biblioteca de problemas

PROBLEMA 1
Resolver el siguiente problema donde todas las variables son positivas o
nulas y se supone una solución inicial con todas las variables iguales a 0.
min a − b − 3c + 3d + 20e + 30 f + 7g + h + i
a −c +d ≤3
b −c −d ≤4
3a −b −3c +2d +e +2 f −h +i ≤1
e +f +g +h ≤1
2e +f −2g +i ≤ 2
PROBLEMA 2
Demostrar que los cortes de Benders son siempre diferentes.
PROBLEMA 3
El método de descomposición de Benders aplicado a un problema es
precisamente el método de descomposición de DW aplicado a su problema dual.
PROBLEMA 4
Resolver el siguiente problema de optimización por el método de DW y RL
partiendo del punto (0, 0, 0, 0) .
04/11/2012 143
12 BIBLIOTECA DE PROBLEMAS
x 1 
 
x 2 
min −2 −1 −3 −1 x 
( )
xi  3 
x 
 4 
x 1 

1 1 1 1x 2  6
    
 1 2 1 x 3  ≤ 4
     
x 4 
1 1  x 1  6
     
 1  x 2  2
    ≤  
 −1 1  x 3  3
     
  

 1 1 x 4  5
xi ≥ 0
PROBLEMA 5
Deducir matemáticamente los cortes de acotamiento para el método de
descomposición de RL cuando el maestro resulta no acotado.
PROBLEMA 6
Demostrar matemáticamente que si en lugar de resolver el problema PL-2 en
su forma primal se resuelve su dual, con el método de descomposición Primal-
Dual resulta la siguiente simetría si t = j : el problema maestro del dual es el
dual del subproblema del primal y el subproblema del dual es el dual del
maestro del primal.
PROBLEMA 7
Programar en Matlab el método de descomposición Primal-Dual.
PROBLEMA 8
Formular un ejemplo de programación lineal estocástica bietapa en
planificación de sistemas de energía eléctrica.
La planificación de la expansión de la generación cuando se considera la
demanda como un parámetro aleatorio. Las decisiones de inversión son las
decisiones de la primera etapa que han de tomarse a priori, sin conocer el
futuro. Las decisiones de la segunda etapa son las de operación del sistema, que
se toman una vez conocido el valor del parámetro aleatorio (la demanda de
144 04/11/2012
energía eléctrica). La función objetivo es la minimización de los costes totales,

costes de inversión y valor esperado de los costes de operación.
PROBLEMA 9
Definir qué condición deben cumplir los incrementos en cada coordenada con
respecto al caso base.
PROBLEMA 10
Comprobar el resultado eligiendo como caso base el estado con todos los
grupos indisponibles.
PROBLEMA 11
Partiendo de la información sobre el universo de estados de un sistema
implantar el método de Monte Carlo y las técnicas de reducción de varianza.
04/11/2012 145
04/11/2012 147

Optimizacion Estocastica

Cargado por

Copyright:

Formatos disponibles

Optimizacion Estocastica

Cargado por

Información del documento

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Optimizacion Estocastica

Cargado por

Copyright:

Formatos disponibles

OPTIMIZACIÓN ESTOCÁSTICA

I.1 INTRODUCCIÓN ........................................................................................... 3

ocurrido un suceso aleatorio. Por ejemplo, se toman hoy un conjunto de

$title Planificación óptima de la expansión de la generación

COST .. COSTE =E= sum(i, F(i) * X(i))

El problema se resuelve en primer lugar de forma determinista para cada

Determinista 1 Determinista 2 Determinista 3 Medio Estocástico

Coste total 262 346.67 437.33 355.73 362.47

Para la solución de problemas estocásticos se ha recurrido históricamente a

i) Una es solucionar el problema para cada posible escenario. Ésta se

El paso de estas aproximaciones a la optimización estocástica supone un

La optimización estocástica permite incorporar actitudes aversas al riesgo,

Figura 1.1 Árbol de probabilidad.

Figura 1.2 Árbol de probabilidad para agregación de escenarios.

La estructura de la matriz de restricciones del problema estocástico se

Figura 1.3 Estructura de la matriz de coeficientes de las restricciones.

La función objetivo del problema estocástico (recourse problema RP) es el

almacenamiento, del tamaño del depósito y de la demanda en cada periodo. La

Escenario Coste del gas (€) Demanda Probabilidad

$title Planificación óptima de la compra de gas

COSTS .. COSTE =E= CVS('nr') * (X('cyv')+X('cya'))

BALDEM2 .. X('ayv') + Y('cyv') =E= DEM ;

GSTDEP .. X('cya') =E= X('ayv') ;

* escenario de demanda media

Los resultados para cada escenario determinista, el escenario medio y el

Compra y Compra y Depósito a Compra y Compra y Coste

I.1.1 Reseña histórica

parámetros. Esta estructura es la que se rompe mediante el uso de las técnicas

• Planificación a largo plazo de la expansión de la red [Perez, Stanford]

• Coordinación hidrotérmica [Morton:93, Jacobs, Morton:96, Pereira]

Otros ámbitos de aplicación pueden ser:

• Planificación de capacidad de una línea de producción [Ramos]

I.2 Optimización lineal bietapa y multietapa

donde A1 ∈ ℝ m1×n1 y A2 ∈ ℝ m2×n2 . Las dimensiones de los demás vectores y

Figura 2.1 Estructura de la matriz de coeficientes de las restricciones en PL-2.

Los dos ejemplos del capítulo anterior corresponden a problemas lineales

Bp−1x p−1 + Ap x p = bp p = 1, …, P (2.2)

El tamaño del problema crece linealmente con el número de etapas. Aunque

I.3 Técnicas de descomposición

• simulación y reducción de varianza

Una técnica de descomposición para la solución de problemas lineales

En la segunda se optimiza sobre una región factible mayor e incrementando las

el tiempo de cálculo se redujo en 37.5 veces. Este decremento depende del

Figura 3.1 Matrices de restricciones diagonales por bloques.

Por ejemplo, un problema multiperiodo con esta estructura en su matriz de

Figura 3.2 Matriz de restricciones de un problema multiperiodo.

En el caso de un problema multidivisional existen restricciones que

Figura 3.3 Matriz de restricciones de un problema multidivisional.

conjunto de periodos. Se define como periodo a la unidad relevante de tiempo

Figura 3.4 Matriz de restricciones del problema de coordinación hidrotérmica.

En la figura 3.5 se observa claramente la estructura en escalera de la matriz

Figura 3.5 Estructura de la matriz de restricciones.

El planteamiento matemático para su resolución se puede hacer utilizando la

Bajo el punto de vista de descomposición de Bd las variables de producción

Figura 3.6 Matriz de restricciones del maestro y subproblema de Benders.

Bajo el prisma de descomposición de DW las restricciones de movimiento de

Figura 3.7 Matriz de restricciones del maestro y subproblema de Dantzig-Wolfe.

La descomposición anidada de Benders es la aplicación recursiva de Benders

Figura 3.8 Matriz de restricciones del maestro y subproblema de Benders anidado.