ARTÍCULO ORIGINAL

Operating costs optimization for steady state regional aquifers

]]>

Lemuel Carlos Ramos Arzola ^I, Eric Cabrera Estupiñán ^II, David Ernesto Marón Domínguez ^III, Armando Orestes Hernández Valdés ^I

^I Centro de Investigaciones Hidráulicas (CIH), Instituto Superior Politécnico José Antonio Echeverría (Cujae), Habana.
^II Universidad Laica «Eloy Alfaro» de Manabí, extensión Chone, Ecuador.
^III Centro de Estudios Matemáticos (CEMAT), Instituto Superior Politécnico José Antonio Echeverría (Cujae), Habana.

RESUMEN ]]> El objetivo de este estudio es desarrollar un modelo matemático para la optimización no lineal de la explotación de un acuífero confinado a escala regional en régimen permanente. Se utiliza AQÜIMPE (en régimen permanente) como modelo de simulación. Se logra la vinculación entrelos modelos de simulación y optimización mediante el enfoque matriz respuesta. Finalmente el modelo de optimización es resuelto por medio de la función «quadprog» del asistente matemático MATLAB. Los resultados obtenidos, en un acuífero hipotético, muestran el esquema de explotación óptima del acuífero que minimiza los costos de bombeo, al mismo tiempo que garantiza restricciones de demanda y de niveles del acuífero.

Palabras clave: AQÜIMPE, costo de bombeo en acuíferos, enfoque matriz respuesta, MATLAB, optimización no lineal.

ABSTRACT

The aim of this study is to develop a mathematical model for nonlinear optimization of well extractions in confined aquifers at a regional scale and steady state flow. AQÜIMPE is used as the simulation model in steady state mode. Simulation and optimization models are linked,which is accomplished by means of there sponse matrix approach. Finally, the optimization model is solved using the MATLAB function «quadprog». Results for a hypothetic aquifer show the optimum well extractions scheme in which pumping costs are minimized while demand and hydraulic head constrains are guaranteed as well.

Keywords: AQÜIMPE, groundwater pumping costs, response matrix approach, MATLAB, nonlinear optimization.

  ]]>  

INTRODUCCIÓN

Las demandas crecientes sobre los recursos de agua subterránea están creando la necesidad de mejorar las técnicas de análisis científico para comprender y administrar mejor los sistemas de agua subterránea (Barlow 2005). Desde la década de los 60, los modelos de simulación numérica han sido importantes herramientas para la valoración y administración de los sistemas de flujo de agua subterránea. Comúnmente, estos modelos eran usados para probar específicas políticas de administración del recurso agua, o, mediante un enfoque de prueba y error, para seleccionar la política que mejor satisfacía los objetivos de la administración (Barlow 2005). El uso de tal enfoque esquiva la formulación rigurosa de los objetivos de la administración del agua subterránea y no puede considerar importantes restricciones físicas y operacionales (Gorelick 1983).

Debido a la compleja naturaleza de los sistemas de agua subterránea y al gran número de factores que a menudo afectan su administración, el proceso de seleccionar la mejor política puede ser extremadamente difícil. Es por ello que para abordar esta dificultad, los modelos de simulación del agua subterránea han sido «conectados» con técnicas de optimización, para determinar la mejor (óptima) estrategia de administración (Barlow 2005). Los modelos que vinculan un modelo de simulación con un modelo de optimización para identificar estrategias de administración óptimas se conocen como modelos de simulación-optimización o modelos de administración (Barlow et al. 1996). Según Ahlfeld and Mulligan (2000), estos modelos explícitamente tienen en cuenta los objetivos y las restricciones de la administración del recurso agua. Los principales enfoques utilizados para llevar a cabo la vinculación entre el modelo de simulación y el modelo de optimización han sido el enfoque Matriz Respuesta (MR) y el enfoque Embebido (Peralta et al. 1991). En el primero mediante el uso del principio de superposición se obtiene una matriz que permite representar la respuesta del acuífero bajo diferentes patrones de explotación. Por su parte el segundo enfoque, «incrusta» todas las ecuaciones del modelo de simulación como restricciones del modelo de optimización y de esta manera el modelo de simulaciónqueda completamente embebido dentro del problema de optimización. El presente trabajo se fundamenta en el enfoque MR.

En Cuba desde la década de los 80 se dispone de un modelo matemático para la simulación del flujo subterráneo, AQÜIMPE (Martínez 1989), con amplias aplicaciones dentro y fuera del territorio nacional. Esta tecnología ha sido ampliada recientemente con un modelo matemático para la administración óptima de la explotación de acuíferos llamado MADA (Cabrera 2009). Según Cabrera y Dilla (2011), MADA está basado en la inclusión de AQÜIMPE dentro de un algoritmo de programación lineal, y por tanto solo está capacitado para resolver problemas donde la función objetivo y las restricciones sean funciones lineales. La función que minimiza los costos de operación de la explotación regional de acuíferos es una función no lineal (cuadrática), es por ello que en la presente contribución se propone un modelo matemático para la optimización no lineal de la explotación de un acuífero confinado a escala regional en régimen permanente. El modelo de administración propuesto utiliza el enfoque MR y resuelve el problema de optimización mediante la función «quadprog» del asistente matemático MATLAB. Finalmente se aplica el modelo de administración en un acuífero hipotético donde se obtiene el esquema de explotación óptima que garantiza, además, restricciones de demanda y de niveles en el acuífero.

MODELO DE SIMULACIÓN

El modelo de simulación AQÜIMPE se encuentra detallado en Martínez (1989) y más recientemente en Cabrera y Dilla (2011). En el presente trabajo se utiliza una modificación de dicho modelo considerando régimen permanente. La ecuación diferencial para este caso es la siguiente:

donde: h es la carga hidráulica [L], b espesor saturado del acuífero [L], K conductividad hidráulica [LT^-1], y W término fuente de entrada o salida vertical [LT^-1]. El término W puede venir dado por la expresión (2).

]]> con Q_s gasto superficial o recarga [LT^-1], Q_pj gasto puntual (o explotación) en el punto j [L³T^-1], y d es la función delta de Dirac [L^-2].

La ecuación (1), sujeta a un conjunto de condiciones de frontera, es resuelta numéricamente mediante el Método del Elemento Finito con triángulo cuadrático y utilizando la aproximación de Galerkin. La solución de esta ecuación permite obtener la carga hidráulica en cualquier punto (x,y) del acuífero.

ENFOQUE MATRIZ RESPUESTA

La vinculación entre el modelo de simulación y el modelo de optimización se lleva a cabo mediante el enfoque de Matriz Respuesta (MR). En el caso del régimen permanente, el enfoque MR utiliza el principio de superposición en el espacio para obtener la respuesta de la carga hidráulica al bombeo (Peralta et al. 1991). Según Riefler and Ahlfeld (1996) en un acuífero confinado la relación entre la carga y el bombeo es lineal, por tanto la MR es una representación exacta de la respuesta del acuífero. A continuación se describen los principales aspectos de este enfoque.

Ahlfeld and Mulligan (2000) plantean que la respuesta de la carga en un punto del acuífero producto de un cambio en el bombeo, puede ser descrita usando la expansión de la serie de Taylor de primer orden:

donde: h_i(Q) es la carga en el punto i para un nuevo vector de gastos Q[L], h_i⁰ es la carga en el punto i cuando no existe gasto en los puntos que se van a optimizar [L], dh_i/dQ_j es el cambio de la carga en el puntoipor un cambio del gasto en el punto j [L/L³T^-1], Q es un nuevo vector de extracciones con componentes Qj [L³T^-1], yn es el número de puntos con gasto a optimizar.

Los términos dh_i/dQ_j se denominan coeficientes respuesta del acuífero (Gorelick 1983), (Ahlfeld and Mulligan 2000). Los coeficientes respuesta proveen información de la respuesta del flujo del agua subterránea a cambios en el gasto de extracción, por tanto el vínculo entre el modelo de simulación y el modelo de optimización es mediante estos coeficientes (Ahlfeld and Mulligan 2000). El método más utilizado para calcular los coeficientes respuesta se denomina método de la perturbación (Ahlfeld et al. 2005), y consiste en utilizar una diferencia finita adelantada, según la expresión (4), para aproximar la derivada de la carga respecto al gasto.

donde: ΔQj es una perturbación o incremento del gasto en el punto j [L³T^-1], h_i(Q_ΔQ) es la carga en el punto i luego de la perturbación en el punto j, y Q_ΔQ es un vector de gastos que tiene todas sus componentes nulas exceptuando la componente j-ésima que es igual a la perturbación ΔQj. Es bueno señalar que como la respuesta de la carga al gasto es lineal (en acuíferos confinados) la derivada parcial en (3) es constante y la expresión (4) solo está sujeta a errores de redondeo. ]]> La matriz formada por los coeficientes respuesta es la llamada MR del acuífero. El cálculo de esta matriz requieren +1 corridas del modelo de simulación (en este caso, de AQÜIMPE). En la primera corrida, considerando que no existen extracciones en los puntos a optimizar, se obtienen las cargas h_i⁰en cada uno de los puntos donde se desea controlar la carga. Luego, en las restantes ncorridas, se determinan las cargas h_i(Q_ΔQ) y se calculan todos los coeficientes respuesta según (4) para conformar la MR. Esta MR, como se verá más adelante, permitirá transformar el modelo de optimización en uno donde las variables de decisión son solo los gastos en los puntos a optimizar.

MODELO DE PROGRAMACIÓN CUADRÁTICA

El asistente matemático MATLAB está dotado con la función «quadprog» que permite resolver un problema de programación cuadrática como el siguiente:

donde: x es el vector de las variables de decisión, H es la matriz Hessiana de la función objetivo, f es un vector que representa los términos lineales de la función objetivo, A_ineq es la matriz de las restricciones lineales de desigualdad, b_ineq es el vector de los términos independientes de las restricciones lineales de desigualdad, A_eq es la matriz de las restricciones lineales de igualdad, b_eq es el vector de los términos independientes de las restricciones lineales de igualdad, y los vectores lb y ub son las cotas inferiores y superiores, respectivamente, de las variables de decisión.

La función «quadprog» es usada para resolver el modelo de administración propuesto y su sintaxis puede presentar la siguiente estructura:

donde: x o pt puede ser un punto de mínimo local para problemas no convexos y, para problemas convexos, ser un punto de mínimo global. El argumento de salida fval es el valor mínimo local o global de la función objetivo.

Por otra parte, la matriz Hessiana H permite determinar la naturaleza del óptimo encontrado durante la optimización. Si H tiene la diagonal predominante, quiere decir que el valor absoluto de cada elemento de la diagonal principal de la fila i es mayor que la suma de los valores absolutos de los restantes elementos de esa fila, y los elementos de la diagonal principal son positivos, se puede decir que el problema es convexo. Entonces si el problema de programación cuadrática es convexo, el mínimo es global (Ahlfeld and Mulligan 2000). Por último, la matriz Hessiana se define como la matriz de las segundas derivadas de la función objetivo F, y cada elemento de la matriz tiene la forma:

]]> FORMULACIÓN DEL MODELO DE ADMINISTRACIÓN

Función objetivo

La función objetivo de los costos de bombeo o de operación es una función que depende del producto de la carga y el gasto, por tanto se trata de una función cuadrática. Esta función objetivo se muestra a continuación:

donde: F es la función objetivo de los costos de operación [$/T], c es el coeficiente de costo del bombeo, expresado como el producto del peso específico del agua [FL^-3] y la tarifa o costo unitario de la potencia eléctrica [$/FL], y Lj son las distancias entre la superficie del terreno y el fondo del acuífero en cada pozo [L].

Como se dijo anteriormente, con el enfoque matriz respuesta se puede transformar la función (8) en una función que solo dependa del gasto. Esto se lleva a cabo sustituyendo la expresión (3) en (8) como se muestra a continuación:

Se observa en (9) cómo ahora F solo es función del gasto Q. Además, gracias a que F solo depende de Q, es relativamente sencillo obtener la matriz Hessiana de esta función. Los elementos de la matriz Hessiana de F tienen la forma:

Ahlfeld and Mulligan (2000) plantean que en la práctica los mayores valores de Hkj ocurren cuando k=j. Esto se debe a que el mayor abatimiento ocurre en el propio punto donde se bombea, y por tanto la derivada respecto al gasto es grande (en términos absolutos, ya que dh/dQ< 0). Entonces es bastante probable que la matriz Hessiana tenga la diagonal predominante (y positiva) y F sea una función convexa, garantizando que el óptimo sea global. ]]> Restricciones

Las restricciones utilizadas en el modelo que se propone pueden ser de tres tipos: restricción de demanda, restricción de capacidad instalada y restricción de niveles. A continuación se desarrolla cada restricción.

La restricción de demanda consiste en que la suma de los gastos en todos los puntos a optimizar debe igualar o superar una demanda prefijada Q_dem. Esta restricción genera una inecuación de la forma:

o en notación matricial:

Según Cabrera y Dilla (2011) el gasto de extracción en cada punto a optimizar no debe sobrepasar la capacidad máxima instalada, Q_j^u, en cada pozo, aunque puede entregar una demanda prefijada Q_j^l. Estas restricciones generan inecuaciones de la forma

Estos dos tipos de restricciones, ecuaciones (11) y (13), están en función dela variable de decisión de gasto, por lo que su incorporación en el modelo de optimización es bien sencilla. Sin embargo, como se verá a continuación, en las restricciones de niveles es necesario el uso de la MR para lograr que estas restricciones solo sean funciones del gasto.

Los puntos donde se controlará la carga deben cumplir que estén por encima de un valor prefijado h_j^l. Este límite inferior puede servir para controlar caídas excesivas de la carga y evitar fenómenos de subsidencia, para controlar los niveles de un cuerpo superficial en comunicación con el acuífero o los niveles próximos a la costa en los acuíferos costeros. Esta restricción genera las siguientes inecuaciones: ]]>

Debe observarse que estas restricciones no están limitadas a los puntos que se van a optimizar, sino que pueden ser aplicadas tanto en puntos a optimizar como en puntos cualesquiera. Como se planteó arriba es necesario transformar las restricciones de niveles para que sean funciones del gasto. Nuevamente, esto se lleva a cabo mediante la MR sustituyendo (3) en (14) y luego acomodándose de la siguiente manera:

Esta nueva expresión (16) genera un conjunto de inecuaciones lineales, que escrito en notación matricial resulta:

Modelo de administración

El modelo de administración propuesto consiste en la formulación del problema de programación cuadrática compuesto por la ecuaciones (9), (11), (13) y (16), que es resuelto por la función «quadprog» de MATLAB. Este modelo solo es válido para acuíferos confinados (o libres que puedan ser tratados como confinados), flujo bidimensional y régimen permanente. La secuencia de pasos para implementar el modelo de administración se muestra en la figura 1 y luego se aplica el modelo propuesto en un acuífero hipotético.

CASO DE ESTUDIO

En esta sección se presenta un modelo hipotético de un acuífero tomado de Martínez et al. (2000) y es sometido a varias modificaciones. Se trata de un acuífero confinado de aproximadamente 1,5 km² de área, con dos zonas de propiedades hidrogeológicas que coinciden con dos zonas de infiltración. Las condiciones de contorno son (ver figura 2A): un contorno AB de carga conocida (95 m), un contorno CD que representa un río, y un contorno EF que descarga al exterior con un gasto lateral fijo (100 m²/d). En aras de obtener una representación gráfica del problema de optimización, se considerarán solamente dos pozos de bombeo como puntos de optimización, nodos 522 y 662, ver figura 2A. ]]>

En la figura 2A se muestra también la malla de triángulos utilizada en la discretización (384 triángulos y 817 nodos entre principales y secundarios). También se presenta en la figura 2B un perfil en la zona de los pozos a optimizar. En las tablas 1 y 2 se dan otros datos adicionales para la reproducción de este modelo hipotético.

RESULTADOS NUMÉRICOS

El caso estudiado consiste en encontrar la explotación óptima en los puntos a optimizar que minimice el costo de bombeo, garantizando una demanda de 800 L/s y no se produzcan niveles inferiores a 50 m en el río modelado por los nodos 273, 548 y 734 (ver figura 2A).

Además se dispone de una capacidad instalada en cada pozo de 500 L/s.Siguiendo el diagrama de flujo presentado en la figura 1 se realizó una primera simulación del acuífero sin considerar gasto en los puntos a optimizar. Los resultados de esta corrida permitieron obtener las cargas h_i⁰ en los nodos a optimizar y en los nodos a controlar la carga (273, 548 y 734) y se muestran en la tabla 3.

Luego se realizaron dos corridas adicionales que permitieron calcular los coeficientes respuesta en los nodos a optimizar y en los nodos a controlar la carga. Los coeficientes respuesta se muestran en la tabla 4 a continuación. Obsérvese que los mayores valores de los coeficientes respuesta se obtienen en los puntos que se optimizan.

]]> A partir de estos resultados se puede formular el problema de optimización, ecuaciones (18 a 23). Es importante señalar que en las ecuaciones (18-23) todas las unidades han sido convertidas a metros y días según corresponda. Como se trata de dos variables de decisión, es posible una representación de la región factible y del óptimo global (figura 3). Este aspecto será demostrado más adelante a partir de la matriz Hessiana.

Sujeto a las restricciones:

Finalmente se resolvió el problema de optimización (18-23) con la función «quadprog» y se obtuvieron los siguientes resultados: Q₅₂₂=469, 43 L/s, Q₆₆₂ = 330, 57 L/s y F = 33684, 77 $/año. A partir de estos resultados se muestra en la figura 4 un mapa de hidroisohipsas. En este mapa se observa que los niveles en el río permanecen iguales o superiores a los niveles mínimos especificados. También se puede comprobar que la suma de los gastos óptimos permite suplir la demanda declarada anteriormente. Por último, la matriz Hessiana del problema se muestra a continuación en (24) observándose que tiene la diagonal predominante y positiva, y por tanto el óptimo encontrado es global.

]]>

CONCLUSIONES

Se concibe e implementa (de forma muy básica) un modelo matemático para la administración de la explotación regional de acuíferos que minimiza la función objetivo no lineal de los costos de operación. Este modelo utiliza el enfoque Matriz Respuesta para la vinculación entre el modelo de simulación y el modelo de optimización, y resuelve el problema de optimización mediante la programación cuadrática utilizando el asistente matemático MATLAB. Además, se tienen en cuenta restricciones de demanda, de capacidad instalada y de niveles.

Se obtiene el patrón óptimo de la explotación de los pozos de bombeo en un acuífero hipotético, donde se garantizó que el óptimo encontrado era global. Además, el modelo de manera satisfactoria cumplió con todas las restricciones impuestas.

El modelo propuesto es aplicable solamente en acuíferos hidráulicamente confinados, bajo flujo regional bidimensional y régimen permanente.

]]>

REFERENCIAS BIBLIOGRÁFICAS

Ahlfeld D. P. and Mulligan A. E. «Optimal management of flow in groundwater systems», Academic Press (ISBN9780120448302), San Diego, California, USA. (2000).

Ahlfeld D. P., Barlow P. M. and Mulligan A. E. «GWM-A ground-water management process for the US Geological Survey modular ground-water model (MODFLOW-2000)», Open-File Report 2005-1072:124, ISBN9781288716388, U.S. Geological Survey. USA. (2005).

Barlow P. M. «Use of simulation-optimization modeling to assess regional ground-water systems», Fact Sheet 2005-3095, U.S. Geological Survey. USA. (2005).

Barlow P. M., Wagner B. J. and Belitz K. «Pumping strategies for management of a shallow water table: The value of the simulation-optimization approach», Groundwater, (1996), 34(2): 305-317, ISSN 1745-6584, National Ground Water Association, USA. ]]>

Cabrera E. «Sistema para la administración de la explotación del agua subterránea», Tesis de doctorado, Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico José Antonio Echeverría (Cujae), La Habana, Cuba. (2009).

Cabrera E. y Dilla F. «Modelo de administración de acuíferos: MADA», Tecnología y Ciencias del Agua, (2011), 2(4): 5-24, ISSN 2007-2422, Instituto Mexicano de Tecnología del Agua, México.

Gorelick S. M. «A review of distributed parameter groundwater management modeling methods», Water Res. Research, (1983), 19(2): 305-319, ISSN 1944-7973, American Geophysical Union, USA.

Martínez J. B. «Simulación matemática de cuencas subterráneas: flujo impermanente bidimensional», Monografía, Facultad de Ingeniería Civil, Instituto Superior Politécnico José Antonio Echeverría (Cujae), Ciudad de la Habana, Cuba. (1989).

Martínez J. B., Llanusa H., HernándezA. O. y Dilla F. «Manual del usuario del sistema AQÜIMPE», Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico José Antonio Echeverría (Cujae), La Habana, Cuba. (2000). ]]>

Peralta R. C., Azarmnia H. and Takahashi S. «Embedding and response matrix techniques for maximizing steady-state ground-water extraction: computational comparison» Groundwater, (1991), 29(3): 357-364, ISSN 1745-6584, National Ground Water Association, USA.

Riefler R. and Ahlfeld D. P. «The impact of numerical precision on the solution of confined and unconfined optimal hydraulic control problems», Hazardous Waste & Hazardous Materials, (1996), 13(2): 167-176, ISSN 0882-5696, Mary Ann Liebert, USA.

]]> Recibido: 7 de junio de 2015.
Aprobado: 22 de noviembre de 2015.

Lemuel Carlos Ramos Arzola, Eric Cabrera Estupiñán, David Ernesto Marón Domínguez, Armando Orestes Hernández Valdés. Centro de Investigaciones Hidráulicas (CIH), Instituto Superior Politécnico José Antonio Echeverría (Cujae), Habana, Universidad Laica «Eloy Alfaro» de Manabí, extensión Chone, Ecuador, Centro de Estudios Matemáticos (CEMAT), Instituto Superior Politécnico José Antonio Echeverría (Cujae), Habana, e-mail: lemuel@cih.cujae.edu.cu, e-mail: ecabrerae@gmail.com, e-mail: dmaron@cemat.cujae.edu.cu, e-mail: ahernandez@cih.cujae.edu.cu ]]> 2000 2005 2005 1996 34 2 2 305-317 2009 2011 2 4 4 5-24 1983 19 2 2 305-319 1989 2000 1991 29 3 3 357-364 1996 13 2 2 167-176