Las fuentes de energía no convencionales, también llamadas energías alternativas, han venido alcanzando un notable desarrollo y aceptación en las últimas décadas y su empleo ha sido muy beneficioso para el medio ambiente, al disminuir el uso de combustibles fósiles que aseveran el efecto invernadero y a su vez el cambio climático. La energía eólica constituye una de esas fuentes no convencionales.
Varios autores en Cuba en la primera década del siglo XXI trabajaron en la temática de la energía eólica, fundamentalmente en la dirección de lograr la caracterización del viento para ser empleado como recurso energético. Se pueden señalar los trabajos de (Proenza, 2004), (Soltura et al., 2006), (Carrasco, 2008) y (Álvarez et al., 2009b), (Álvarez et al., 2009c), (Álvarez et al., 2009a), (León, 2010), (Roque et al., 2010), (Perdigón and Rodríguez, 2011) y (Martínez, 2012). Más recientemente se obtuvo, en el 2014, como parte de un proyecto conjunto con Canadá, el Atlas Eólico de Cuba con una resolución de 2 km y niveles de altura de 10 m, 30 m, 50 m, y 100 m (Roque and Yu, 2014).
Actualmente, según (Roque and Yu, 2014), el país tiene en su proyección instalar para el año 2030 más de 300 MW. Un total de 13 parques eólicos asumirán esa producción energética ((Rosell, 2015) citado por (Roque et al., 2015b) ) y se prevé para el 2020 que la generación eléctrica a partir de energía eólica represente el 6 % en el sistema eléctrico cubano.
Dada la intermitencia y variabilidad del recurso eólico, el Sistema Electroenergético Nacional debe estar preparado para asumir la generación de energía eléctrica a partir de fuentes convencionales en momentos en que esta no se pueda generar a partir de la fuente eólica. La forma en que esto puede garantizarse es contando con pronósticos del viento en los lugares de instalación de los parques eólicos. En el caso específico de tener que realizar operaciones en la red eléctrica en tiempo real, garantizando el balance y la operación óptima de las reservas, los pronósticos a muy corto plazo, que van desde pocos segundos hasta 30 minutos, son muy útiles.
En el ámbito internacional las investigaciones dirigidas a pronosticar a muy corto plazo la rapidez del viento con fines energéticos o la potencia eólica propiamente son en número menores que aquellas enmarcadas en el corto plazo. Aun así los métodos utilizados para pronosticar en ambos plazos son los mismos, basados principalmente en series de tiempo o en inteligencia artificial, o la combinación de ambos en métodos híbridos, pueden citarse entre los más recientes los trabajos de (Jiang et al., 2013), (Xiaodan et al., 2013), (Carpinone et al., 2015), (A.Sapronova et al., 2015),(Filik., 2016), (Li et al., 2016); también puede encontrarse la aplicación de modelos de pronóstico numéricos, el método de persistencia, método del vecino más cercano en (Appice et al., 2015), wavelet y redes neuronales en (Senkal and Ozgonenel, 2013) y otros métodos, pero en menor medida.
En Cuba, relacionados a las predicciones con fines eólicos se cuenta con tres estudios anteriores, el de (Hinojosa, 2015), donde utilizando un modelo ARMA y redes neuronales artificiales el autor obtuvo pronósticos a muy corto plazo de la variable viento, pero en dos zonas de ubicación de torres eólicas, no propiamente para el lugar de emplazamiento de un parque eólico en funcionamiento e incorporado al Sistema Electroenergético Nacional; y (Roque et al., 2015b) y (Roque et al., 2015a), donde se generaron pronósticos a corto plazo de la rapidez del viento con el modelo numérico Weather Research and Forecast (WRF) para sitios de emplazamiento de torres eólicas y de la potencia eólica para los parques eólicos Gibara y Los Canarreos, respectivamente.
Lo anterior pone de manifiesto la problemática de que no existen pronósticos de la potencia eólica en la escala de muy corto plazo para ningún parque eólico del país, situación que generó la siguiente interrogante: ¿Cómo pronosticar a muy corto plazo la potencia eólica esperada del Parque Eólico de Gibara I?
En este artículo se pretende, utilizando un enfoque estadístico, responder la interrogante anterior, teniendo como objetivo general obtener a través de métodos estadísticos pronósticos energéticos a muy corto plazo de la potencia eólica que generará el Parque Eólico Gibara I y como específicos: a)obtener un modelo ARMA (Modelo Autorregresivo de Medias Móviles) para pronosticar la rapidez del viento a muy corto plazo en el Parque Eólico Gibara I, en el nivel de 50 m de altura; b) evaluar la calidad de los pronósticos emitidos por el modelo a partir del cálculo de diversos estadísticos; c) combinar el pronóstico de rapidez del viento obtenido con la curva de potencia del parque para obtener el pronóstico energético.
El modelo AR(p) (Casimiro., 2009) es una aproximación natural al modelo lineal general expresado por:
Donde es el valor de la variable en el momento , en el momento y así sucesivamente.
Este es un proceso mediante el cual se expresa en función de sus valores pasados hasta el retardo más una innovación contemporánea :
El modelo de medias móviles finito MA(q) es también una aproximación natural al modelo lineal general y se obtiene truncando el modelo de medias móviles de orden infinito. Expresa el valor de en función de la innovación contemporánea y de su pasado hasta el retardo q:
A través de los procesos autorregresivos de medias móviles ARMA (p,q) se puede determinar en función de sus valores pasados hasta el retardo p, de la innovación y de su pasado hasta el retardo q de la siguiente forma:
y son los parámetros del modelo.
Este modelo comparte las características de los modelos AR (p) y MA (q) debido a que contiene ambas estructuras a la vez. El modelo ARMA (p, q) tiene media cero, varianza constante y finita y una función de autocovarianzas infinita. La función de autocorrelación es infinita decreciendo rápidamente hacia cero, pero sin truncarse.
(Sclove, 1987) citado por (Díaz, 2011) plantéo que la aplicación de criterios para la selección de modelos podría ayudar a tratar má expeditamente problemas que se tratan mediante secuencias de test de hipótesis. En modelos donde se trabaja con parámetros estimados por métodos de máxima verosimilitud se pueden utilizar estos criterios de información, y su aplicación es extensible a la situaciones en las que la meta es la selección de modelos adecuados, siendo el caso los modelos de series temporales.
El criterio básico entre aquellos basados en la información estadística es el criterio de Akaike, que en sus comienzos de desarrolló en series temporales. El mismo no es más que un estimador muestral de la esperanza de la log-verosimilitud. Su expresión es:
Es decir, depende de la función de verosimilitud de las observaciones y del número de parámetros .
Para seleccionar el modelo el criterio de Akaike se basa en la maximización de la log-verosimilitud de un modelo esperado utilizando la metodología de la máxima verosimilitud. Se selecciona el modelo que alcance el menor valor del criterio.
Para comprobar la existencia de autocorrelaciones en una serie de tiempo existen varias pruebas estadísticas basadas en los residuos. Estas pruebas pueden ser para determinar atocorrelaciones de primer orden, entre las que se hallan la de Durbin -Watson, la razón de Von Neumann y la prueba de Berenblut & Webb. Para órdenes superiores pueden mencionarse la de (Godfrey, 1978) y (Godfrey, 1988) y la de (Box and Pierce, 1970) citados por (Ramírez), cuyo estadístico se contrasta con los valores críticos de la tabla de Chi-cuadrado bajo la hipótesis nula de ausencia de autocorrelación. Este estadístico fue refinado por Ljung y Box en 1979 obteniendo la variante:
La zona de Gibara en la provincia de Holguín fue seleccionada como región de estudio, específicamente el sitio donde se halla enclavado el mástil de prospección eólica de Los Cocos (Figura 1). El mismo se halla a 300 metros de la línea costera y presenta una altura sobre el nivel del mar de 3 m (Roque., 2015).
La selección de esta torre de prospección se basó en la cercanía que presenta la misma al Parque Eólico Gibara I, el cual no posee mástiles de prospección instalados dentro de su perímetro, por lo que las condiciones físico-geográficas son similares al igual que el comportamiento del régimen de vientos sinópticos y locales. Por tanto, la selección realizada contribuirá a que el pronóstico represente lo más cercanamente posible el comportamiento futuro de la rapidez del viento en las zonas donde se hallan instalados los aerogeneradores.
Los datos utilizados corresponden a observaciones cada 10 minutos en el período desde el 1ero de julio de 2015 a las 00:00:00 hasta el 9 de diciembre de 2015 a las 09:50:00, lo que equivale a un total de 23244 observaciones. El nivel de medición empleado fue el de 50 m de altura. Se emplearon 11680 datos (hasta el 28 de septiembre de 2015 a las 23:50:00 horas) como serie de tiempo (en adelante serie de tiempo) para la obtención del modelo ARMA. Los datos restantes (desde el 29 de septiembre de 2015 a las 00:00:00 horas hasta el 9 de diciembre de 2015 a las 09:50:00 horas, en lo adelante período de comparación) ayudaron a la comparación de la salida de los modelos con información real para su validación.
Para la obtención del modelo ARMA se trabajó en Matlab 8.3. El primer paso consistió en representar gráficamente el comportamiento promedio de la rapidez del viento cada 10 minutos con la data real en la zona del Parque Eólico Gibara I, para su posterior comparación con las salidas del modelo. Una vez realizado esto, se procedió a obtener los correlogramas de la función de autocorrelación fac y la función de autocorrelación parcial facp con un retardo de hasta 50, utilizando las funciones autocorr y parcorr, con el objetivo de identificar los posibles órdenes p y q del modelo ARMA. Una vez identificados estos posibles valores se escribió un código con el que, utilizando las funciones armax y aic se obtuvieron con la primera de ellas modelos ARMA (p,q) con p y q desde 1 hasta 50 y con la segunda el criterio de Akaike para cada uno de esos modelos, generándose en total 2500 modelos y el valor del criterio para cada uno de ellos. Los valores del criterio de Akaike se almacenaron en una matriz de datos donde las columnas representan el valor de p y las filas el valor de q. Dentro de la matriz se identificó el valor más pequeño y los órdenes de p y q correspondientes. La función armax de Matlab utiliza para la estimación de los parámetros del modelo el método de predicción-error descrito en (Ljung, 1999).
Para valorar el cumplimiento de los supuestos de normalidad y de no existencia de correlación entre los residuos del modelo se obtuvo, utilizando la función probplot , la gráfica de probabilidad con ajuste a una distribución normal y nuevamente con la función autocorr la de autocorrelación de los mismos.
La prueba Ljung-Box para aproximadamente un 95% de intervalo de confianza, utilizando un retardo también hasta 50, fue realizada al modelo ajustado para evaluar si los residuos fueron independientes entre sí, utilizando la función lbqtest.
Con la instrucción modelo. ParameterVector./sqrt(diag(modelo. CovarianceMatrix)) se obtuvo un vector con los valores del estadístico de student, utilizado para verificar la hipótesis de que los parámetros del modelo fuesen significativamente diferentes de cero, con un 95% de nivel de confianza.
En un segundo código nuevamente con la función armax se obtuvo la ecuación del modelo ARMA (p, q) con los valores de p y q obtenidos para el criterio de Akaike más pequeño. Dentro de ese mismo código se obtuvieron además los residuos del modelo para comprobar posteriormente si se comportaron como ruido blanco, cuestión que de no ser así podría conllevar a plantearse un nuevo modelo.
Después, usando la función compare de Matlab 8.3 se graficaron los valores pronosticados del modelo cada 10 minutos. Esta función tiene como datos de entrada la serie real con la que se va a comparar y el modelo identificado, y genera la serie pronóstico para su comparación con la real.
Los pronósticos cada 10 minutos con los modelos ARMA fueron llevados a una hoja de cálculo en Excel 2010 donde se promediaron los mismos horarios del día para obtener un único gráfico pronóstico promedio cada 10 minutos para el período de comparación. Todos los datos promedio, tanto reales como los pronósticos, fueron representados gráficamente utilizando Excel 2010.
Para evaluar los pronósticos realizados y siguiendo la metodología de trabajo con los promedios explicada en el párrafo anterior, se calculó el Error Medio Absoluto (EMA), así como el sesgo o BIAS para determinar la tendencia del modelo a subestimar o sobrestimar en los pronósticos. La ecuación para el EMA:
Donde es el valor pronosticado, el valor real y N la cantidad de datos empleados.
La ecuación para el cálculo del BIAS:
Se obtuvo además el coeficiente de correlación entre los datos reales y los pronosticados con la función Coef. de correl. de Excel.
Con el pronóstico de rapidez del viento realizado se calculó, utilizando la curva de potencia del parque eólico, el pronóstico energético para el mismo. La curva de potencia de los aerogeneradores del parque eólico (Roque et al., 2015a) suministrada por el fabricante se muestra en la Figura 2. En la ecuación representa la potencia y la rapidez del viento.
La rapidez del viento promedio cada 10 minutos en 50 m de altura en el período de comparación mostró valores que superaron los 6 m/s en el 93.75% de los horarios del día analizados. Se observó una tendencia a la disminución de la rapidez del viento en horas de la madrugada y hasta las 07:00:00 hora local, momento en el cual la rapidez del viento comenzó a ascender, produciéndose el valor más elevado a las 13:00:00 hora local y pasadas las 17:00:00 una tendencia a su disminución, comportamiento similar al estudiado por (Martínez, 2012) para torres eólicas costeras del país. Los valores obtenidos y la persistencia del viento en las horas donde mayores valores se alcanzaron confirman estudios anteriores que establecen la zona de estudio como de las más aprovechables en el país en cuanto a la disponibilidad del recurso viento.
La fac con retardo hasta 50 para la región de estudio en el nivel de medición analizado se muestra en la Figura 4. Se observó cómo la fac fue disminuyendo suavemente a medida que aumentó el retardo, mostrando la existencia de diferencias leves en la autocorrelación de los distintos retardos hasta el valor de 50.
En cuanto a la facp se apreció que existían varios coeficientes diferentes de cero, por tanto el comportamiento de la fac y facp llevó a considerar un modelo ARMA (p,q) a estimar con órdenes superiores al ARMA (1,1).
Una vez determinados los valores del criterio de Akaike para los distintos modelos con p y q variando valores desde 1 hasta 50 para la serie de datos del nivel de altura en cuestión, estos se almacenaron en una matriz de 50x50. Dentro de la matriz de cada nivel se identificó el modelo con el menor valor del criterio de Akaike, resultando que el modelo ARMA (23,17) el de menor valor del criterio de Akaike, coincidiendo esto con lo evidenciado por el comportamiento de la fac y facp , que sugerían la utilización de modelos de órdenes superiores al ARMA (1,1) para el pronóstico a muy corto plazo de la variable rapidez del viento.
La gráfica de probabilidad para la distribución normal de los residuos del modelo ARMA (23,17) se muestra en la Figura 5 En ella la curva negra discontinua representa la distribución normal. Si los residuos tuviesen este tipo de distribución, los puntos de la curva azul hubiesen coincidido o estado muy próximo a la línea discontinua, que representa el ajuste a la distribución normal. Como se observa la curva de los residuos tomó una forma de S, lo que sugiere que los mismos, para el modelo en cuestión, no se distribuyeron normalmente.
Sin embargo, la fac de los residuos, que se puede observar en la Figura 6, evidenció que no existieron valores significativos de autocorrelación para los residuos del modelo ajustado, con menos del 5% de ellos fuera de las bandas de Bartlett que indican aproximadamente el 95% de intervalo de confianza. Por tal motivo se concluyó la no existencia de correlación entre los residuos, y por tanto que los mismos asemejaron a un ruido blanco.
En la Tabla 1 se muestran los valores calculados en Matlab 8.3 para el estadístico Q con un 95% de intervalo de confianza y un retardo de 50 en el nivel de medición en cuestión. En el caso de Matlab 8.3, el valor del estadístico Q puede ser 0 o 1. El valor nulo indica aceptar la hipótesis nula de no existencia de correlación entre los residuos, y el valor 1 lo contrario.
Se observa que el estadístico no rechazó la hipótesis de la no existencia de correlación entre los residuos de cada uno de los modelos identificados, lo que se reafirmó por el Pvalor que fue superior a 0,05. Por tanto, se aceptó que el modelo es adecuado para la finalidad del pronóstico a muy corto plazo.
La prueba de significación de los parámetros del modelo, efectuada a través del estadístico t-student, arrojó los valores que se muestran en la Tabla 2 para el modelo obtenido. Para que los parámetros sean significativamente diferentes de cero deben ser, en valor absoluto, mayores que 1.96. Como se aprecia, existieron varios parámetros que no cumplieron esta condición, lo que sugirió debía rechazarse el modelo identificado.
Hasta el momento, los residuos del modelo solo cumplieron el supuesto de no correlación entre ellos, no así el de ajuste a una distribución normal. Por otra parte, hubo varios parámetros del modelo que no fueron significativamente diferentes de cero. Todo lo anterior se consideró debió estar relacionado a la no linealidad del proceso. Tales resultados pudieron haber conllevado a la decisión de rechazar el modelo ARMA (23,17) para el pronóstico a muy corto plazo de la rapidez del viento en el sitio en cuestión. No obstante, se consideró el hecho de que los procesos en la atmósfera no tienen un comportamiento normal y además, que las relaciones entre los mismos son no lineales, y en esta investigación se está representando el proceso de forma lineal mediante el modelo autorregresivo. Por tanto, se decidió analizar su ajuste a los datos reales, para valorar si a pesar del incumplimiento de los supuestos anteriores, se obtenía un pronóstico aceptable.
En la Figura 7 se muestra el pronóstico a muy corto plazo de la rapidez del viento realizado por el modelo ARMA (23,17) así como la serie de datos reales del período de comprobación en el Parque Eólico Gibara I. Es apreciable que existió buena similitud entre ambas curvas, sobre todo entre las 13:00:00 y 17:10:00 hora local, horario este último en el que persisten las velocidades más altas del viento como promedio y por tanto el de mayor generación de potencia eólica.
El EMA calculado por el modelo en el pronóstico de la rapidez del viento (Figura 8) osciló entre 0.83 y 1.93 m/s, estando solo en el 10.4% de los pronósticos efectuados por debajo de 1.0 m/s, en el 77.8% entre 1.0 y 1.5 m/s y en el restante 11.8% superior a ese último valor. Coincidiendo con lo expresado en el párrafo anterior, los mayores EMA (superiores a 1.5 m/s) se sitúan en horarios de la madrugada y mañana y posteriores a las 17:50:00 horas.
En el 90.9% de las previsiones realizadas por el modelo ARMA (23,17) los datos pronosticados tuvieron respecto a los reales coeficientes de correlación superiores a 0.8, lo que se puede discernir en la Figura 9 a continuación, demostrándose la buena correlación existente entre el pronóstico y la realidad.
En el 38.9% de los casos el modelo subestimó en el pronóstico, mientras que en 33.9% sobrestimó y en el 27.8% de los casos el modelo prácticamente no tuvo sesgo (Figura 10), por lo que se concluye en este aspecto que no existe una marcada diferencia en el comportamiento del BIAS que permita decidir la tendencia en el sesgo del modelo durante el pronóstico para los próximos 10 minutos.
El pronóstico energético se obtuvo aplicando la curva de potencia de los aerogeneradores instalados en el Parque Eólico Gibara I a los valores de rapidez del viento pronosticados por el modelo ARMA (23,17) en los próximos 10 minutos, para el nivel de 50 m de altura, y se compararon con la potencia eólica obtenida de aplicar la misma curva a los valores reales de rapidez del viento en el período de comprobación. La Figura 11 muestra los resultados obtenidos para la potencia eólica y su respectivo EMA en el nivel de 50 m de altura, donde se pudo constatar que tal como sucedió con el pronóstico de la rapidez del viento, existió similitud entre la curva de potencia real y la pronosticada.
El EMA de la potencia (Figura 11b) osciló en el 6.2% de los horarios entre 300-400 kW, en el 24.3% entre 400-500 kW, en el 42,4% entre 500-600 kW, en el 25.7% entre 600-700 kW y en solo el 1.4% rebasó los 700 kW. Considerando que la dependencia entre la potencia y la rapidez del viento es exponencial, lo que provoca que a medida que aumente el valor de esta última variable errores bajos en su pronóstico conllevan a errores elevados en la predicción energética, y que se realizó la predicción en un lugar donde como promedio la rapidez del viento superó los 6 m/s, poseer la mayoría de los pronósticos (72.9%) con errores inferiores a 600 kW puede ser valorado como un buen pronóstico, y es consistente con los resultados obtenidos por (Roque et al., 2015a) para este mismo parque eólico en el pronóstico a corto plazo con el modelo WRF.
Se obtuvo un Modelo Autorregresivo de Medias Móviles para pronosticar a muy corto plazo la rapidez del viento en el Parque Eólico Gibara I, con valores de p= 23 y q=17. El EMA de los pronósticos realizados osciló entre 0.83 y 1.93 m/s, estando solo en el 10.4% de los pronósticos efectuados por debajo de 1.0 m/s, en el 77.8% entre 1.0 y 1.5 m/s y en el restante 11.8% superior a ese último valor, lo que se registraron en horarios de la madrugada y la mañana y posteriores a las 17:50:00 horas.
El coeficiente de correlación r entre el pronóstico de rapidez del viento dado por el modelo y los datos reales fue superior a 0.8 en más del 90.0 % de los pronósticos realizados.
En el análisis del BIAS no se logra identificar una tendencia del modelo a subestimar, sobrestimar o no presentar sesgo en el pronóstico, visto de manera general. Por tanto, este tipo de análisis deberá realizarse para cada pronóstico por separado, a fin de identificar la tendencia del modelo en cada horario y considerarlo entonces en la toma de decisiones.
El pronóstico energético obtenido para el parque eólico Gibara I tuvo un EMA inferior a 600 kW en el 72.9% de las previsiones realizadas, lo que comparado con otros resultados para este parque eólico puede considerarse aceptable.
Dados los resultados de este trabajo se recomienda la aplicación de modelos ARMA para pronosticar la rapidez del viento a muy corto plazo en este período del año en el Parque Eólico Gibara I.
Aplicar modelos ARMA para pronosticar la rapidez del viento a muy corto plazo en períodos de tiempo disímiles a los de este resultado, a fin de comprobar su aplicación en distintas épocas del año.
Utilizar modelos capaces de representar la no linealidad del proceso, tales como aquellos basados en redes neuronales artificiales, para el pronóstico a muy corto plazo en el Parque Eólico Gibara I y comparación con los resultados obtenidos en el presente trabajo, o una combinación de ambos tipos de modelos para tales efectos.
Considerar los resultados alcanzados en el momento de realizar cualquier operación en tiempo casi real en el Parque Eólico Gibara I.