+ All Categories
Home > Documents > DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

Date post: 12-Nov-2021
Category:
Upload: others
View: 5 times
Download: 0 times
Share this document with a friend
18
DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES DOCUMENTO DE TRABAJO - DW-DT-018-003 Jesús María Velásquez Bermúdez DecisionWare Ltd., Colombia [email protected] Caracas, 1981 (Revisado 2012)
Transcript
Page 1: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN

SISTEMAS CON MÚLTIPLES EMBALSES

DOCUMENTO DE TRABAJO - DW-DT-018-003

Jesús María Velásquez Bermúdez

DecisionWare Ltd., Colombia

[email protected]

Caracas, 1981

(Revisado 2012)

Page 2: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

2

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN

SISTEMAS CON MÚLTIPLES EMBALSES

Jesús María Velásquez Bermúdez

DecisionWare Ltd., Colombia

[email protected]

RESUMEN

El artículo presenta los fundamentos teóricos para la determinación de la distribución óptima

del recurso hídrico en un sistema con múltiples embalses. Este enfoque ha sido denominado

operación con base en "metas hidrológicas", que corresponden a superficies

multidimensionales en las que el sistema esta en un punto de equilibrio óptimo de acuerdo con

el agua, o a la energía, almacenada agregadamente entre todos los embalses.

La metodología propuesta se fundamenta en el análisis estocástico del conjunto de parámetros

que definen los puntos críticos de operación del sistema a lo largo del tiempo, en lo que se

refiere a déficits y a superávits. Estos parámetros se analizan por medio de modelos de

optimización estocástica que definen niveles de los embalses que en términos de valores

esperados optimizan la operación del sistema. El análisis se realiza separadamente para

déficits y para superávits, una vez se ha determinado las regiones críticas para cada caso. Los

conceptos que se presentan pueden aplicarse a sistemas de recursos hidráulicos con múltiples

propósitos; como la generación de hidroelectricidad, el riego, el abastecimiento de agua

potable, y/o a control de contaminación.

ABSTRACT

The article presents the theoretical bases for the determination of the optimum water resource

distribution in a multireservoir system. This approach has been designated as operation based

on " hydrological targets", that correspond to multidimensional surfaces in which the system is

in an optimum equilibrium point for the water, or the energy, stored in the reservoirs.

The methodology proposed is based on the stochastic analysis of the set of parameters that

define the critical operation points of the system in the medium and long term, with respect to

deficits and to surpluses. These parameters are analyzed by means of stochastic optimization

models to define the amount of water in the reservoirs that in mean values terms optimize the

operation of the system. The analysis is accomplished severally for deficits and for surpluses,

once it has been determined the critical regions for each case. The concepts that are presented

can be applied to water resource systems with multiple purposes, such as hydroelectricity

generation, irrigation, drinkable water supply, and/or control of pollution.

Page 3: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

3

1. INTRODUCCION

A continuación, se presenta la política de operación de sistemas de recursos hídricos con

múltiples embalses denominada "metas hidrológicas", que se asocia a la distribución óptima de

agua, o de energía, almacenada en los diferentes embalses. Las "metas hidrológicas" definen

superficies multidimensionales sobre las que el sistema se encuentra en un punto de equilibrio

óptimo. La formulación inicial de las metas hidrológicas, realizada en 1975 por Durán et al.

[2], presenta simplificaciones que han sido revisadas en varios trabajos posteriores de

Velásquez [6], Velásquez y Durán [7], y Velásquez [8]. La metodología básica que se presenta

y revisa en este artículo fue utilizada en el diseño y la planificación de la operación del sistema

de riego del río Guárico en Venezuela (PDC Ingeniería [4]).

La determinación de las metas hidrológicas se fundamenta en el estudio de las características

hidrológicas del sistema de múltiples embalses, dejando de lado su valoración económica. Esto

es consonó con el hecho probado que la distribución óptima del recurso hídrico es

independiente de su valor de oportunidad como recurso competitivo con otros recursos

(Velásquez [9]). Las metas hidrológicas se implementan en dos pasos: inicialmente se determina

la distribución óptima del recurso hídrico, esta fase se denomina diseño; posteriormente las

metas hidrológicas se utilizan para simular o para operar el sistema, esta fase se denomina

aplicación.

La metodología propuesta se fundamenta en el análisis estocástico del conjunto de parámetros

que definen los puntos críticos de operación del sistema a lo largo del tiempo, en lo que se

refiere a déficits y a superávits. Estos parámetros se analizan por medio de modelos de

optimización estocástica que definen estados que en términos de valores esperados optimizan la

operación del sistema. El análisis se realiza separadamente para déficits y para superávits, una

vez se ha determinado las regiones críticas para cada caso.

El punto de partida de las metas hidrológicas es el uso de la teoría de procesos estocásticos

aplicada al diseño y a la operación de embalses, uno de los tópicos más desarrollados en

hidrología. Para el caso de sistemas con un embalse, el problema ha sido analizado

profundamente, llegándose a soluciones analíticas para series de tiempo estocásticas que

cumplen con estrictas hipótesis con respecto a sus características probabilísticas. En el caso de

series multivariadas se ha recurrido a técnicas numéricas para realizar el análisis. La teoría de

metas hidrológicas es elaborada a partir del análisis estocástico de la hidrología de cada uno

de los embalses, el cual produce información que es condensada de acuerdo a la topología de

las redes, en parámetros que definen zonas críticas en lo referente a desbordamientos y a

déficits. El paso final está basado en un proceso de optimización que define distribución óptima

del agua en el sistema, metas hidrológicas, y sus características estadísticas bajo condiciones

de optimalidad.

Las metas hidrológicas no solo producen los números apropiados para la operación, si no que

aumentan el conocimiento del proceso hidrológico que afecta a un sistema de manejo de

recursos hídricos.

Page 4: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

4

2. CONCEPTUALIZACIÓN DE UN SISTEMA CON MÚLTIPLES EMBALSES

A continuación, se sintetiza la descripción matemática de la topología del sistema. Los

componentes a considerar son: embalses, canales, redes de embalses y puntos de demanda de

agua a los embalses y/o a las redes. Los embalses se consideran como las unidades primarias

del sistema, en tanto que las redes están formadas por la unión de embalses y/u otras redes. A

continuación, se presenta los parámetros básicos del sistema, en ellos el índice i esta asociado a

los embalses y el j a las redes.

¡Error! Marcador no definido.PARÁMETROS BÁSICOS DE UN SISTEMA CON MÚLTIPLES EMBALSES

NOMBRE DESCRIPCIÓN

CONJUNTOS QUE DESCRIBEN LA TOPOLOGÍA DEL SISTEMA

EMB Conjunto de embalses en consideración.

RED Conjunto de redes en consideración.

RE(i) Conjunto de las redes que contienen al embalse i.

ER(j) Conjunto de los embalses contenidos en la red j.

RR(j) Conjunto de las redes contenidas en la red j.

PR(j) Conjunto de las redes que contienen a la red j.

ERP(j) Conjunto de embalses que no están contenidos en la red j.

CAPACIDADES DE ALMACENAMIENTO Y UTILIZACIÓN

CAi(t) Capacidad de almacenamiento del embalse i en el período t.

CUi(t) Capacidad de utilización del embalse i durante el período t.

CURj(t) Capacidad de utilización de la red j durante el período t.

DEMANDA DE AGUA A REDES Y EMBALSES

DEi(t) Demanda propia del embalse i durante el período t, que solo puede atenderse con agua proveniente del embalse i

DRi(t) Demanda propia de la red j durante el período t. que puede ser atendida indistintamente con agua proveniente de

cualquiera de los componentes de la red.

RANGOS DE VARIACIÓN DE LA DEMANDA

DEmaxi(t) Máxima demanda de agua al embalse i durante el período t

DEmaxi(t) = DEi(t) + jRE(i) DRj(t)

DEmini(t) Mínima demanda de agua al embalse i durante el período t

DEmini(t) = DEi(t)

DRmaxj(t) Máxima demanda de agua a la red j durante el período t

DRmaxj(t) = DRminj(t) + kPR(j) DRk(t)

DRminj(t) Mínima demanda de agua a la red j durante el período t

DRminj(t) = DRi(t) + iER(j) DEi(t) + rRR(j) DRr(t)

DEMANDAS MAXIMAS EFECTIVAS

SDEmaxi(t) Máxima demanda efectiva que puede ser satisfecha con agua proveniente del embalse i en el tiempo t

SDEmaxi(t) = MIN{Demaxi(t), CUi(t)}

SDEmini(t) Mínima demanda efectiva que puede ser satisfecha con agua proveniente del embalse i en el tiempo t

SDEmini(t) = MIN{Demini(t), CUi(t)}

SDRmaxj(t) Máxima demanda efectiva que puede ser satisfecha con agua proveniente de la red j en el tiempo t

SDRmaxj(t) = t DRmaxj(t)

SDRminj(t) Mínima demanda efectiva que debe ser satisfecha con agua proveniente de la red j en el tiempo t

SDRminj(t) = t DRminj(t)

3. DISEÑO DE LA POLÍTICA DE METAS HIDROLÓGICAS

El análisis matemático está orientado a la determinación de la distribución óptima del recurso

hídrico para un instante dado en el tiempo, debiéndose repetir el proceso para todos los

momentos en los que se debe realizar, o simular, la operación.

Page 5: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

5

3.1. ANÁLISIS Y SÍNTESIS DEL PROCESO HIDROLÓGICO

El análisis y la síntesis estocástica se realiza partiendo del análisis detallado individual de cada

embalse, e integrando esta información a lo largo de cada una de la redes del sistema.

3.1.1. ANÁLISIS INDIVIDUAL DE CADA EMBALSE

Si se desea operar un embalse durante un período que comienza en el instante 0 y se extiende t

unidades de tiempo, las posibilidades de déficit, de uso y de desbordamiento de agua se pueden

definir con base en parámetros que definen la región de factibilidad de uso del agua y

funciones de desbordamiento dependientes del nivel inicial y la cantidad de agua utilizada en el

período. Estos parámetros, presentados en la siguiente tabla, son variables aleatorias derivadas

del proceso hidrológico estocástico Qi(t) que aporta agua al embalse i.

PARÁMETROS QUE DEFINEN LAS POSIBILIDADES DE OPERACIÓN DE UN EMBALSE

NOMBRE DESCRIPCIÓN

Tmaxi(t) Máximo volumen de agua que puede utilizarse en el embalse i durante el período de longitud t.

Nmini(t) Mínimo volumen de agua que debe estar almacenado en el embalse i al comienzo del período para poder utilizar Tmaxi(t).

Tmini(t) Mínimo volumen de agua que se debe extraer del embalse i durante el período para evitar desbordamientos, cuando el volumen

inicial de agua almacenada es cero.

Nmaxi(t) Máximo volumen inicial de agua que puede estar almacenada en el embalse i al comienzo del período, por encima del cual se

producen desbordamientos adicionales a los inevitables, cuando se utiliza Tmaxi(t),

i(t) Capacidad de acumulación del embalse i durante el período t.

Ai(t) Alivio inevitable en el embalse i durante el período t.

NIEcriti(t) Mínimo volumen de agua que debe estar almacenado en el embalse i al comienzo del período con el fin de poder satisfacer la

mínima demanda al embalse -SDEmini(t)-.

SDEFi(t) Déficit inevitable en la satisfacción de la demanda propia del embalse i durante el período.

Tmaxi(t) y Nmini(t) determinan la potencia de uso del agua del embalse y se calculan en forma

simultánea. Nmini(t) determina el nivel mínimo a partir del cual el embalse puede ser utilizado

a su máxima potencia, su valor esta determinado por relación entre su tendencia a secarse y su

tendencia a aliviar. Nmini(t) se calcula con base en la integral de la diferencia entre el aporte

Qi(t) y la máxima demanda al embalse SDEmaxi(t), ajustada por la capacidad de

almacenamiento CAi(t). Nmini(t) es una función creciente con respecto a t, acotada por la

capacidad del embalse en el instante 0, CAi(0). Tmaxi(t) determina la máxima potencia

utilizable del embalse y se determina por la relación entre la tendencia a aliviar, la tendencia a

secarse, la capacidad de almacenamiento y el volumen almacenado al comienzo del período.

Las siguientes figuras presentan el proceso de determinación de Nmini(t) y la zona de

posibilidades de uso del agua como función del nivel NIi(0).

Page 6: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

6

PROCESO DE ESTIMACION DE Nmin(t)

POTENCIA NO

UTILIZABLE

Nmin(t1)

t

t2 t3

Nmin(t3)

t0

CAPACIDADA

INSTALADA

t1

Nmin(t)

X(t)

Nmin(t3)

Nmin(t1)

t

0

N(t)

0

C(0)

C(0)

T(t)

Tmax (t)

Nmin (t)

45

o

N(0)

REGIÓN DE POSIBILIDADES DE USO DE AGUA

Nmini(t) determina el volumen por encima del cual todo incremento inicial de agua almacenada

no influye en la capacidad de utilización del embalse. Para valores inferiores a Nmini(t) todo

decremento en el nivel inicial implica decremento de la máxima cantidad de agua que es

posible utilizar. En términos de ecuaciones esta zona se define como:

0 Ti(t) Tmaxi(t)

0 Ti(t) NIi(0) + Tmaxi(t) - Nmini(t)

0 NIi(0) CAi(0)

donde Ti(t) representa la cantidad de agua utilizada durante el período de longitud t, y NIi(0) el

volumen de agua almacenada inicialmente. Tmini(t), Nmaxi(t), i(t) y Ai(t) están relacionados

con las posibilidades de pérdida de agua por razones de desbordamiento del embalse. Nmaxi(t)

está determinado por la tendencia a aliviar cuando el embalse es utilizado a su máxima

potencia. Nmaxi(t) es una función decreciente con respecto a t, cuyo máximo valor es CAi(0).

Tmini(t) determina la mínima utilización que debe hacerse del embalse con el fin de no perder

agua durante el período. Las siguientes figuras presentan el proceso de cálculo de Tmini(t) y de

Nmaxi(t).

Page 7: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

7

PROCESO DE ESTIMACION DE N max(t)

X(t)X(t)

C(0)

0.t

L

AGUA PERDIDA

Nmax(t)

t0

Nmax(t0)

0.t

Nmax(t0)

t1 t2

Tmin (t) SQ(t)

0

C(t)

N(t)

PROCESO DE ESTIMACION DE Tmin(t)

Tmin(t)

t

DESBORDAMIENTO

INEVITABLE

i(t) representa la capacidad de acumulación del embalse y determina la capacidad útil que

esta disponible para propósitos de regulación. El alivio inevitable Ai(t) establece los volúmenes

de agua que no pueden ser controlados. Conjuntamente Nmaxi(t), Tmini(t) y Tmaxi(t)

determinan la zona de posibilidades de alivio evitable Di(t), adicional a Ai(t):

Di(t) = Max{0., NIi(0) - Nmaxi(t) + Tmaxi(t) - Ti(t), NIi(0) + Tmini(t) - Ti(t)}

Las figuras siguientes presentan la región y las isolíneas de desbordamiento adicional, como

función del agua utilizada Ti(t) y del nivel inicial NIi(0).

D(t)

C(0) -N(t)

T(t)

REGION DE DESBORDAMIENTO

ISOLINEAS DE DESBORDAMIENTO

T(t)

Tmax(t)

Tmin(t)

Nmin(t)

N(o)

C(0)Nmax(t)

D(t) = 0

SDEFi(t) y NIEcriti(t) determinan las posibilidades de atender la demanda propia del embalse

i. SDEFi(t) determina el déficit inevitable que se debe a razones diferentes a la operación del

sistema, ya que su origen son las condiciones del proceso hidrológico y/o las capacidades de las

componentes del sistema. SDEFi(t) se define como:

SDEFi(t) = Tmaxi(t) - SDEmini(t)

Page 8: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

8

NIEcriti(t) determina el máximo volumen que puede estar almacenado al comienzo del período

con el fin de no incurrir en déficits evitables, superiores a SDEFi(t), imputables a la operación

del sistema, y cumple con la siguiente definición:

NIEcriti(t) = MAX[0., Nmini(t) - MAX{0., Tmaxi(t) - SDEmini(t)}]

Como paso inicial para el cálculo de los anteriores parámetros consideremos los rangos de

variación de nivel de un embalse, los cuales se definen en la siguiente tabla:

VARIABLES QUE DEFINEN EL RANGO DE VARIACIÓN DEL EMBALSE i

VARIABLE DESCRIPCIÓN FÓRMULA

Qi(t) Caudal afluente al embalse i en el instante t.

SQi(t) Caudal afluente al embalse i acumulado hasta el instante t. SQi(t) = Qi(t) t t Qi(t)

Xi(t) Variación mínima del embalse i en el tiempo t (caudal

descontando la demanda máxima). Xi(t) = {Qi(t) - SDEmaxi(t)}t t {Qi(t) - SDEmax(t)}

Si(t) Serie acumulada de Xi(t) hasta el tiempo t. Si(t) = Xi(t) t t Xi(t)

Mi(t) Máximo valor de Si(t) hasta el tiempo t. Mi(t) = Max{0., Si(1), Si(2), ... , Si(t)} = Max{0., Mi(t-1), Si(t)}

mi(t) Mínimo valor de Xi(t) hasta el tiempo t. mi(t) = Min{0., Si(1), Si(2), ... , Si(t)} = Min{0., mi(t-1), Si(t)}

Ri(t) Rango de variación de Xi(t) hasta el tiempo t. Ri(t) = Mi(t) - mi(t)

Ri*(t) Rango de variación ajustado de Xi(t) hasta el tiempo t

(limitado por la capacidad de almacenamiento del embalse. Ri*(t) = Min{Ri(t), Ci(t)}

Con base en una relación recursiva, el proceso de cálculo de los parámetros básicos es:

CALCULO RECURSIVO DE LOS PARÁMETROS FUNDAMENTALES PARA UNA EMBALSE

PARÁMETRO CONDICIÓN INICIAL RELACIÓN RECURSIVA

Tmini(t) Tmini(0) = 0. Tmini(t) = Max{Tmini(t-1), SQi(t) - Ci(t)}

Nmini(t) Nmini(0) = 0. Nmini(t) = Nmini(t-1) + Min{Ci(t) - Ri*(t-1), MAX[0., mi(t-1) - mi(t)]}

Nmaxi(t) Nmaxi(0) = Ci(0) Nmaxi(t) = Nmaxi(t-1) - Min{Ci(t) - Ri*(t-1), MAX[0., Mi(t) - Mi(t-1)]}

Tmaxi(t) Tmaxi(0) = 0. Tmaxi(t) = Tmaxi(t-1) + Min{SDEmaxi(t), CUi(t), i(t-1) + Qi(t)}

i(t) i(0) = Ci(0) i(t) = i(t-1) + Min{0., Qi(t) - [Tmaxi(t) - Tmaxi(t-1)]}

Ai(t) Ai(0) = 0. Ai(t) = Ai(t-1) + Max{0., Qi(t) - [Tmaxi(t) - Tmaxi(t-1)] - i(t-1)}

3.1.2. ANÁLISIS AGREGADO DEL SISTEMA

Para garantizar la confiabilidad del sistema debe analizarse el comportamiento conjunto de

todos los embalses con el fin de determinar las regiones críticas de operación en un espacio

multidimensional. Los parámetros que definen las regiones críticas para una red son:

PARÁMETROS QUE DEFINEN LAS REGIONES CRÍTICAS PARA UNA RED

PARÁMETRO DESCRIPCIÓN

TAmaxj(t) Máximo volumen de agua proveniente de la red j que puede

utilizarse durante el período de longitud t, Tamaxj(t) = Min[SDRmaxj(t), iERP(j) Tmaxi(t)

+ kRR(j) TAmaxk(t)

NAminj(t) Mínimo volumen inicial de agua requerido en la red j para poder

utilizar Tamaxj(t) entre todos los embalses de la red.

NAminj(t) = Max[0., iER(j) Nmini(t) + kRR(j)

NAmink(t)

+ Tamaxj(t) - iER(j) Tmaxi(t)]

NAcritj(t) Mínimo volumen inicial de agua requerido en la red j para

satisfacer la mínima demanda de agua a la red -SDRminj(t)-, NAcritj(t) = i Tmini(t) + k TAmink(t)

NAmaxj(t) Máximo volumen inicial de agua que puede estar almacenado en

la red j para pasar el período con mínimo alivio posible, Namaxj(t) = Max{0., SDRmaxj(t) - TAminj(t)}

DRFj(t) Déficit inevitable en la red j durante el período t. DRFj(t) = Max{0., SDRminj(t) - TAmaxj(t)}

Page 9: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

9

Los parámetros anteriores definen la región de no déf icit como función del volumen de agua

almacenada inicialmente. Para cualquier combinación de niveles iniciales NIi(0) que viole los

niveles críticos individuales de los embalses, o los agregados de las redes, se producirán

déficits. Las figuras siguientes presentan el caso de una red de dos embalses.

REGION DE DEFICITS PARA DOS EMBALSES

REGION

DE NO

DEFICIT

N2(0)

NAcrit(t)

Ncrit2(t)

Ncrit1(t) NAcrit(t)

N1(0)

DEF(t)=0

DEF (0)

N2(0)N1(0)

NIVELES DE DEFICIT PARA DOS EMBALSES

En la determinación de los niveles críticos es necesario considerar dos aspectos: i) cuando la

red ó embalse actúa en su nivel más elemental, es decir no tiene en cuenta las redes que la

contienen, ii) cuando la red ó embalse actúa dentro de una red. Las expresiones que definen

los niveles críticos del sistema son:

PARÁMETROS QUE DEFINEN LAS REGIONES CRÍTICAS PARA UNA RED

PARÁMETRO DESCRIPCIÓN FÓRMULA

NIEcriti(t) Nivel crítico para el embalse i cuando este es considerado

aisladamente (ya definido) Max[0., Nmini(t) - Max{0., Tmaxi(t) - SDEmini(t)}]

NAEcriti(t) Nivel crítico para el embalse i al considerarlo en todas las

redes a las que pertenece.

Max(0., Nmini(t) - Max[0., Tmaxi(t) -

MaxjRE(i){ 0., SDRminj(t) - TAmaxj(t) + Tmaxj(t)}])

NEcriti(t) Nivel crítico del embalse i, definido como el máximo de las

dos expresiones anteriores. NEcriti(t) = Max[NIEcriti(t), NAEcriti(t)]

NAcritj(t) Nivel crítico en la red j Max{0., NAminj(t) - Max[0., TAmaxj(t)-SDRminj(t)],

iER(j) NEcriti(t)+ kRR(j) NAcritk(t)}

DEFT(t) Déficit total en el sistema. Maxj(0., Max[Max{0., NIRcritj(t) - iER(j) Nii(0)}, Maxi{0.,

NIEcriti(t)} - NIi(0)])

3.1.3. ANÁLISIS ESTOCÁSTICO DE LOS PARÁMETROS

Los parámetros del sistema definidos en las secciones anteriores son funciones de las variables

aleatorias que afectan al sistema: los aportes a los embalses y las demandas a embalses y a

redes. Por lo tanto, se deben caracterizar mediante funciones de distribución multivariadas.

Existen dos vías para tal fin: i) distribuciones analíticas, y ii) distribuciones numéricas. Dada la

complejidad de las expresiones que definen los parámetros y la fuerte estructura de correlación

debida a la estacionalidad de los caudales y de las demandas, se descarta la vía analítica,

realizándose el análisis con base en la generación sintética de datos multivariados. Por

simplicidad solo se estudia el tratamiento estocástico de la hidrología.

Page 10: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

10

Si se consideran NH series sintéticas es posible asociar a cada una de ellas un conjunto

conformado por los valores de los parámetros fundamentales. Adicionando el índice h que se

refiere a uno de los NH posibles valores que puede tomar la hidrología, se tiene:

▪ Parámetros para los embalses: {Tmaxi,h(t), Tmini,h(t), Nmaxi,h(t), Ai,h(t), SDEFi,h(t),

NEcriti,h(t)}

▪ Parámetros para las redes: {TAmaxj,h(t), TAminj,h(t), NAminj,h(t), NAmaxj,h(t), NAcritj,h(t)}

Para series equiprobables, la probabilidad de ocurrencia para cada conjunto de parámetros es

igual a 1/NH.

El análisis estadístico multivariado para los parámetros determina dos áreas:

▪ Zona con probabilidad de déficits igual a "cero" y,

▪ Zona con probabilidad de déficits diferente de cero.

La frontera de la primera zona está definida por los siguientes parámetros:

NCcriti(t) = MAXh(NEcriti,h(t))

NAcritj(t) = MAXh(NAcritj,h(t))

NCcriti(t) y NAcriti(t) constituyen los niveles críticos del sistema desde el punto de vista de

déficits una vez se ha integrado la información acerca de su variabilidad estocástica.

3.2. OPTIMIZACIÓN ESTOCÁSTICA

3.2.1. OPERACIÓN EN LA ZONA DE NO DÉFICIT

Cuando se está operando en la zona de no déficits el problema principal está relacionado con

la minimización de las pérdidas del recurso hídrico lo que se consigue resolviendo un problema

de optimización que minimice los alivios (desbordamientos) totales en el sistema.

3.2.1.1. MODELAJE MATEMÁTICO

Consideremos el problema de optimización estocástica PP(t,NA) que minimiza el valor

esperado de los alivios en el sistema, para un volumen agregado de agua almacenada NA, con

un horizonte de planificación de operación que se extiende t unidades de tiempo, cuando se

consideran NH escenarios hidrológicos equiprobables:

PP(t,NA): =

Minimizar

Z(t,NA) = [ iEMB h=1,NH Di,h]/NH

Sujeto a:

Ti,h - NIi Tmaxi,h - Nmini,h

0 Ti,h Tmaxi,h

SDEmini,h - DEFi,h Ti,h

SDRminj,h - DRFj,h iER(j) Ti,h SDRmaxj,h - DRFj,h

iEMB, h=1,NH

iEMB, h=1,NH

iEMB, h=1,NH

iEMB, h=1,NH

Ti,h + Di,h - NIi Tmini,h

Ti,h + Di,h - NIi - Nmaxi,h + Tmaxi,h

iEMB, h=1,NH

iEMB, h=1,NH

Page 11: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

11

NCcriti NIi

Nacritj iER(j) NIi

iEMB

jRED

iEMB NIi = NA

0 NIi Ci iEMB

Ti,h, Di,h 0 iEMB, h=1,NH

donde NIi(0) representa el nivel inicial del embalse i, y Di,h el alivio mínimo en el embalse i si

ocurre la serie hidrológica h, Ti,h la cantidad óptima de agua a utilizar en cada embalse.

PP(t,NA): resuelto para valores de NA que varían entre NAcritNR y la capacidad agregada

total ( iEMB CAi(0)) proporciona la distribución óptima de agua entre los embalses (NI*) en la

zona de no déficits, como función del volumen agregado de agua almacenada al comienzo del

período NA, para un horizonte de planificación de longitud t.

3.2.1.2. MODELAJE DE GRAN ESCALA

En forma matricial PP(t,NA) se puede formular como:

PP(t,NA): { Minimizar h=1,NH CTXh A Xh + B NI [:] bh D NI [:] b , 0XhUh , LNNI UN

}

donde [:] representa un operador de relación que puede ser , =, ó según sea el caso. El

vector Xh está referido al conjunto de variables que simulan agregadamente la operación del

sistema para serie sintética h, y bh y Uh los parámetros de la misma serie. El vector NI contiene

los volúmenes iniciales de los embalses. Los vectores LN y UN corresponden a los vectores de

niveles críticos y de capacidades. La estructura matricial de PP(t,NA) es dual-angular y

permite utilizar la teoría de Benders [1] [3] para la solución del problema.

ESTRUCTURA MATRICIAL DE PP(NA)

X1 X2 . . . XNH NI

A 0 . . . 0 B : b1

0 A . . . 0 B : b2

. . . . . . .

. . . . . . .

. . . . . . .

0 0 . . . A B : bNH

0 0 . . . 0 D : B

Con base en Benders la solución de PP(t,NA) se realiza en dos niveles. En el primer nivel de

resuelve un problema coordinador CNI: que controla el valor de los niveles iniciales NI, y en el

segundo, nivel de descentralización, se resuelve un problema relacionado con las NH series

hidrológicas. Teniendo en cuenta que la estructura matricial del problema sobre las variables

Xh es del tipo diagonal, el problema de optimización sobre Xh se puede descomponer en NH

subproblemas SPh(NI): todos ellos función de NI y de la forma:

Page 12: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

12

SPh(NI): = { Minimizar ZX(NI) = h=1,NH CT XhA Xh : bh - B NI , Uh Xh 0 }

El modelo coordinador CNI: se formula como

CNI: = { Minimizar Z = Q

Q h=1,NH Qh

Qh (Whk)T [bh - B NI] h=1,NH, kITF(h) ;

(Whk)T [bh - B NI] 0 h=1,NH, kITN(h) ;

D NI b ; LN NI UN }

donde Whk representa las variables duales obtenidas en la k-ésima solución del subproblema

SPh(NIk):, ITF(h) el conjunto total de iteraciones para el subproblema h, e ITN(h) el conjunto

de las iteraciones en que no se ha conseguido la factibilidad. El siguiente diagrama presenta un

diagrama del esquema del algoritmo de Benders.

COORDINACIÓN

Coordinador CNI

NIk W1k NIk W2

k NIk Whk NIk WNH

k

SP1(NIk) SP2(NIk) . . . SPh(NIk) . . . SPNH(NIk)

D E S C O M P O S I C I O N

ALGORITMO DE PARTICIÓN Y DESCOMPOSICIÓN DE BENDERS.

El proceso de solución se puede interpretar como:

▪ En el nivel de coordinación se proponen valores de NI como potencialmente óptimos: y

▪ En el nivel de descomposición se evalúa el impacto de la solución propuesta por el

coordinador para cada posible serie hidrológica, y se envía información al coordinador a

través de las variables duales asociadas.

Para cualquier solución a SPi(NIk) se generarán valores factibles de Wh que aportan

información para restringir la zona de optimalidad de NI. En caso de que no exista solución

factible a SPh(NIk), la información generada afecta la zona de factibilidad de NI por razones

de sus relaciones con la zona de factibilidad de las Xh, debiéndose generar un corte de acuerdo

al Lema de Farkas para zonas factibles lineales [3].

Dado que PP(t,NA): debe resolverse paramétricamente para diferentes valores de NA, es

conveniente analizar su estructura con el fin de acelerar la solución conjunta de todos los

problemas. Si se tiene en cuenta que la generación de variables Whk factibles al dual del

problema agregado es independiente de los valores NA, un mecanismo de aceleración consiste

en que la primera solución para NI, al comenzar el proceso para un valor de NA que no es el

primero, es obtenido teniendo en cuenta todos los valores de Whk factibles en iteraciones

Page 13: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

13

previas, independientemente del valor de NA, y de que la restricción hubiera sido o no activa.

Adicionalmente, al considerar los valores óptimos de los niveles iniciales como función de NA,

NIi* = NIi*(NA)

se puede demostrar un comportamiento continuo y monotónicamente creciente de NIi*(NA):

NIi*(NA) NIi*(NA+NA)

La anterior conclusión genera un corte adicional para problemas PP(t,NA): posteriores al

primero, modificando la restricción de cota inferior para NIi:

Max [NCcriti, NIi*(NA-NA)] NIi*

3.1.2.3. SÍNTESIS ESTOCÁSTICA

La solución de PP(t,NA): por el método de descomposición propuesto presenta no solo ventajas

computacionales en tiempo y en memoria, si no que permite un análisis teórico conceptual

desde el punto de vista estocástico. Al considerar que la solución a SPh*[NI] da la respuesta

óptima del sistema condicionada en la distribución de agua en los embalses, y en la serie

sintética h, es posible, integrando la información sobre todas las series sintéticas, determinar la

función de distribución de Z (alivios del sistema integrado), así como rangos de variación e

intervalos de confianza, lo que permite analizar márgenes de riesgo para las soluciones

propuestas.

Otra variable importante es la suma de las variables duales de las tres restricciones donde

aparece el nivel inicial (NIk), que en adelante llamaremos hk, que cumple con

E [Z[NI]/NI | NI=NIk , ESCENARIO h] = h

k

infiriéndose que la unión de los diferentes valores de hk determina la función de distribución

para el gradiente de la función objetivo con respecto al nivel del embalse, lo que en la práctica

da la tendencia del sistema a aliviar como función de los niveles de los embalses. Dada la

estructura matricial de SP*h[NI] los posibles valores que puede tomar las componentes del

vector hk son 0 ó 1 lo que implica que el comportamiento de estas variables aleatorias es del

tipo Bernoulli y por lo tanto su valor esperado puede interpretarse como la probabilidad de que

una variación NIi sea aliviada por el sistema, información de gran utilidad para la operación

en tiempo real. Esta probabilidad se puede definir como

E [Z(NI)/NI | NI=NIk ] = h=1,NH hk

3.2.2. DISEÑO DE LA OPERACIÓN EN LA ZONA DE DÉFICIT

Page 14: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

14

Cuando el sistema se encuentra en la zona de déficit, su operación depende fundamentalmente

del equilibrio entre los déficits causados en el momento de la operación y el compromiso que se

adquiere con los posibles déficits futuros. Desde este punto de vista el diseño de las metas

hidrológicas debe orientarse a estimar, para diferentes niveles agregados de agua almacenada,

los mínimos déficits esperados en el sistema, como función de la distribución de agua. Teniendo

en cuenta las definiciones dadas previamente los anteriores objetivos se alcanzan al resolver el

problema:

PD(t,NA): = {

Minimizar Z = [ h=1,NH Max(0., Max(0., NAcritj(t) - iER(j) NIi), i Max(0., NIEcriti(t) -

NI)]/NH

| i=1EMB NIi = NA ; 0 NIi(0) Ci(0) }

PD(t,NA) debe resolverse para valores de NA entre 0 y el nivel crítico agregado NACcrit(t),

proporcionando el valor del mínimo déficit esperado como función de NA, y de la distribución

óptima de agua. El modelo planteado corresponde uno de programación lineal que se resuelve

siguiendo los mismos lineamientos utilizados para el caso de operación en la zona de no déficit.

4. APLICACIÓN DE LA POLÍTICA DE METAS HIDROLÓGICAS

Este numeral considera la aplicación de las metas hidrológicas en un modelo de planificación

operativa orientado a determinar las decisiones de operación (reales o simuladas) para un

periodo [t,t+t]. Se debe tener en cuenta que las metas hidrológicas sintetizan toda la

información con respecto al futuro en el mediano y en el largo plazo. Bajo este enfoque la

operación óptima esta orientada a: i) determinar la cantidad agregada de recurso hídrico que

se va a utilizar durante el periodo, y ii) llevar los embalses a un nivel predeterminado al final

del período como consecuencia de las decisiones tomadas en el paso previo. Estas decisiones

deben ser el resultado de un estudio que contemple los siguientes aspectos:

▪ Costos de los déficits de demanda causados.

▪ Valor esperado de los déficits de demanda en el futuro.

▪ Costos y beneficios del uso del recurso hídrico

▪ Ajuste del nivel de los embalses a las metas hidrológicas.

▪ Operación en la zona no crítica.

Para lograr lo anterior se puede formular un modelo de optimización multiobjetivo con base en

los siguientes supuestos:

▪ Que en el instante t, en que se va a planificar la operación del período [t,t+t], se conoce el

valor de los siguientes parámetros:

▪ NIi Nivel inicial del embalse i

▪ Qi(t) Predicción de la escorrentía en al embalse i el período [t,t+t]

▪ Las restricciones correspondientes al modelo de operación del sistema real que se resumen

como XNI,Q ▪ La función objetivo se maneja mediante un esquema de optimización jerárquico

multiobjetivo de acuerdo con los siguientes criterios: en primera instancia se tienen en

Page 15: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

15

cuenta los aspectos económicos relacionados con los déficits presentes y futuros, en la

segunda fase se determina los movimientos de embalses teniendo en cuenta que estos se

ajusten a sus metas hidrológicas y respeten sus niveles críticos.

En la optimización se deben considerar funciones objetivos diferentes de acuerdo a cada

instancia del proceso. En el esquema jerárquico, cada vez que se resuelve un nivel superior el

valor de la función objetivo se convierte en restricción para los niveles inferiores, de manera

tal que se respete el nivel logrado para los objetivos de mayor nivel en la de optimización de los

niveles inferiores. Los niveles propuestos son:

▪ minimización de los déficits en el período [t,t+t] más la ponderación de los déficits en

futuros períodos

▪ ajuste de los embalses a sus niveles críticos al final del período de planificación

▪ minimización de los alivios del sistema en el período [t,t+t]

▪ ajuste de los embalses a sus metas hidrológicas de acuerdo al nivel agregado final del

sistema al final del período [t,t+t]

Los anteriores criterios pueden cambiar según el sistema.

Consideremos la función objetivo en cada caso:

▪ Minimización del costo por déficit

Para determinar las decisiones optimas el problema a resolver es:

{ Minimizar w1 = CDEF(DEF(t)) + CDF(NA(t+t))

| XNI,Q NA(t+t) = iEMB NFi(t+t) }

donde CDEF(t) representa los costos por déficit en el periodo t, DEF(t) el déficit en el

periodo t y CDF(NIVA(t+t)) la función de costo esperado de los déficits futuros

dependiente del nivel agregado final, que se obtiene en la fase de diseño de la política de

operación al solucionar la familia de problemas PD(t,NA).

▪ Ajuste de los embalses a los niveles críticos

El ajuste de los embalses a sus niveles críticos, desde el punto de vista de déficits, se

consigue resolviendo un problema que incluye como restricción la definición del nivel

agregado final de cada red y se fija como restricción el valor de la función objetivo

obtenida el paso anterior (w1*).

{ Minimizar w2 = iEMB Max(0, NCriti(t+t) - NFi(t+t))

+ rRED Max(0, NRCritr(t+t) - NFRr(t+t))

| w1* = CDEF(DEF(t)) + CDF(NA(t+t)) ,

XNI,Q NFRr(t+t) = iER(j)NFi(t+t) NA(t+t) = iEMB NFi(t+t) }

▪ Minimización de los vertimientos.

Se deben minimizar los alivios en el sistema, y fijando como restricción el valor de la

función objetivo obtenida en el paso anterior (w2*):

Page 16: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

16

{ Minimizar w3 = iEMB VEi(t,t+t) |

w1* = CDEF(DEF(t)) + CDF(NA(t+t)) ,

w2* = iEMB Max(0, NCriti(t+t) - NFi(t+t)) + rRED Max(0, NRCritr(t+t) - NFRr(t+t)) ,

XNI,Q NFRr(t+t) = iER(j)NFi(t+t) NA(t+t) = iEMB NFi(t+t) }

▪ Ajuste de los embalses a sus metas hidrológicas.

El ajuste de los embalses a sus metas de operación lo que se consigue resolviendo un

problema en el que se fija el valor de la función objetivo en el paso anterior (w3*):

{ Minimizar | iEMB i | NFi - NDi(t+t,NA)| |

w1* = CDEF(DEF(t)) + CDF(NA(t+t)) ,

w2* = iEMB Max(0, NCriti(t+t) - NFi(t+t)) + rRED Max(0, NRCritr(t+t) - NFRr(t+t)) ,

w3* = m VE m(t,t+t)

XNI,Q NFRr(t+t) = iER(j)NFi(t+t) NA(t+t) = iEMB NFi(t+t) }

donde NDi(t,NA) representa la meta hidrológica para el embalse i en el instante t dado un

nivel agregado NA, y i corresponden a pesos para la violación de la metas hidrológicas

los cuales pueden evaluarse a partir de las variables duales del problema de determinación

de las metas hidrológicas que, como se indicó, están relacionados con las probabilidades

de alivio y de déficit. En los casos más triviales pueden asumirse iguales a 1. es una

potencia de ajuste, para la cual se proponen dos valores1 o 2. En el primer caso el

problema se mantiene lineal (si todas las restricciones son lineales) y en el segundo el

problema se convierte en uno de programación cuadrática. Por simplicidad del análisis en

adelante solo consideraremos el caso lineal. En este caso el valor absoluto introducido en

la función objetivo puede reemplazarse como la suma de dos variables positivas, esto es

NFi - NDi (t+t,NA) = Yi+ - Yi

-

y se reformula la función objetivo como:

Minimizar iEMB i (Yi+ + Yi

-)

Esta transformación, dadas las características de colinealidad entre los vectores Yi+ y Yi,

representa exactamente al problema original ya que el producto Yi+ por Yi

- es igual a cero.

5. IMPLEMENTACIÓN EN EL SISTEMA DE RIEGO GUARICO-ORITUCO

La metodología de las metas hidrológicas se utilizó para evaluar diferentes alternativas de

desarrollo de un sistema de riego y de abastecimiento de agua potable en las cuencas de los ríos

Guárico y Orituco en el zona central de Venezuela. El proyecto tenía como objetivo determinar

la conveniencia de:

▪ Realizar un proyecto en el río Orituco (embalse y/o canal de trasvase al Guárico) y

determinar las características de diseño del mismo;

Page 17: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

17

▪ Construir los embalses de la cuenca media del río Guárico: Taguay, Tinapuy y Portachuelo;

▪ Construir los embalses de la cuenca alta del río Orituco; y

▪ Aumentar la capacidad instalada para bombeo de agua subterraea en la zona de riego del

Guárico.

La siguiente gráfica presenta la conectividad del máximo sistema que se podía desarrollar.

MEMO GUANAPITOPIEDRA

PINTADA

TAGUAYTINAPUY

CAMATAGUA

PORTACHUELO

GUARICOORITUCO

AGUA

SUBTERRANEA

SISTEMA DE RIEGO GUARICO-ORITUCO

La metodología de análisis se basó en la evaluación de los costos y rendimientos marginales,

comparados con los costos y rendimientos de un sistema de referencia constituido por los

embalses de Guárico, Camatagua y Guanapito, asumiendo que la demanda del sistema se ha

estabilizado y desarrollado a su máxima capacidad.

El proceso de simulación del sistema Guárico-Orituco se realizó con base en la coordinación de

tres modelos:

▪ GENSIN: modelo de generación sintética de caudales basado en la metodología propuesta

por Valencia&Schake[5]

▪ METHID: modelo que determina las metas hidrológicas

▪ SIMMET: modelo que simula la operación con base en las metas hidrológicas.

Los tres modelos se utilizan coordinadamente de acuerdo al diagrama que se presenta a

continuación.

Page 18: DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN …

DISTRIBUCIÓN ÓPTIMA DEL RECURSO HÍDRICO EN SISTEMAS CON MÚLTIPLES EMBALSES

18

MODELO DE SIMULACIÓN DEL SISTEMA CON MÚLTIPLES EMBALSES

SERIES

SINTETICAS

# 1

METHIDMETHIDDISEÑO

POLÍTICA

OPERACION

OPTIMA

METAS

HIDROLOGICAS

DE

OPERACIÓN

TOPOLOGÍA Y

PARÁMETROS DE

CONTROL

SIMMETSIMMETSIMULACIÓN

POLITICA

OPERACIÓN

ÓPTIMA

COSTOS

RENDIMIENTO

SERIES

SINTETICAS

# 2

GENSINGENSINGENERADOR

SINTETICO

CAUDALES

REFERENCIAS

[1] Benders, J.F. "Partitioning Procedures for Solving Mixed Variables Programming

Problems". Numer. Math 4, 238-252 (1962).

[2] Durán, H., et. al. "A Model for Planning Hydrothermal Power Sustyems”. P.I.C.A.

Conference USA (1975).

[3] Lasdon, L. “Optimization Theory for Large Systems”. MacMillan Publishing Co. (1970).

[4] P.D.C. Ingeniería S.A. “Selección de la Alternativa Optima de Abastecimiento

Superficial al Embalse del Guárico”. Documento DGI/IT/147. Ministerio del Ambiente y

de los Recursos Naturales Renovables. Caracas, Venezuela (1981).

[5] Valencia., D. and Shaake, J. "A Disaggregation Model for Time Series Analysis and

Synthesis". MIT Report No. 149 (1972).

[6] Velásquez, J., "Modelo de Evaluación de Alternativas de Expansión de un Sistema

Hidrotérmico, No-interconectado, con Múltiples Embalses". Tesis de Grado de Magister

en Ingeniería Industrial, Universidad de Los Andes, Bogotá, COLOMBIA (1975).

[7] Velásquez, J. y Durán, H., "Operación de Sistemas Hidrotérmicos con Múltiples

Embalses". Seminario sobre Análisis de Sistemas de Recursos Hídricos de la Conferencia

Internacional de la O.N.U. sobre el Agua, Mar de Plata, ARGENTINA (1977).

[8] Velásquez, J., "Operación de Sistemas Interconectados con Múltiples Embalses. Un

Análisis Integral". Trabajo de Ascenso a Profesor Agregado, Universidad Simón Bolívar,

Caracas, VENEZUELA (1981).

[9] Velásquez, J., "El Efecto de las Tasas de Descuento en la Distribución Optima del

Recurso Hídrico". Revista ENERGÉTICA No. 24, Universidad Nacional de Colombia

Sede Medellín, Colombia (2000).


Recommended