Modelo predictivo de radiación solar mediante técnicas de machine learning: aplicación a la isla de Gran Canaria

Departamento de Ingenier´ıa El´ectrica Modelo Predictivo de Radiaci´on Solar mediante t´ecnicas de Machine Learning. Aplicaci´on a la isla de Gran Canaria Solar Radiation Forecasting Model using Machine Learning techniques. Application to Gran Canaria Island Tesis Doctoral Luis Mazorra Aguiar Las Palmas de Gran Canaria Octubre 2015 Departamento de Ingenier´ıa El´ectrica Programa de doctorado: Tecnolog´ıa Insdustrial Modelo Predictivo de Radiaci´on Solar mediante t´ecnicas de Machine Learning. Aplicaci´on a la isla de Gran Canaria Solar Radiation Forecasting Model using Machine Learning techniques. Application to Gran Canaria Island Tesis Doctoral Autor Director Director Luis Mazorra Aguiar Felipe D´ıaz Reyes Philippe Lauret Las Palmas de Gran Canaria Octubre 2015 Agradecimientos Aprovecho este momento para agradecer la ayuda y el apoyo recibido durante la elaboraci´on de esta tesis. En primer lugar me gustar´ıa agradecer la dedicaci´on y apoyo recibida de mis directores Felipe D´ıaz Reyes y Philippe Lauret durante la elaboraci´on de este trabajo. En especial querr´ıa agradecer al Dr. Felipe D´ıaz la confianza depositada en mi para escribir esta tesis. Su ayuda ha sido constante desde la organizaci´on del trabajo hasta la b´usqueda de financiaci´on para poder realizar la estancia en la Universidad de la Reuni´on. No s´olo ha sido un apoyo acad´emico durante la investigaci´on sino un amigo en quien confiar durante estos meses de duro trabajo. Menci´on especial tambi´en a la inestimable ayuda del Dr. Philippe Lauret, no s´olo durante mi estancia en la Universidad de la Reuni´on sino siempre que la diferencia horaria nos lo permitiera. Durante mi estancia ofreci´o sin condiciones toda su experiencia y trabajo en este campo y gui´o mi investigaci´on para conseguir este resultado final independientemente del tiempo y esfuerzo que conllevara. Adem´as quisiera agradecerle su compa˜n´ıa durante mi estancia, su tiempo para ense˜narme su maravillosa isla y que me recibiera con los brazos abiertos en su casa. No quisiera olvidarme de todas las personas a las que conoc´ı en el Laboratoire de Physique et Ing´enierie Math´ematique pour l’Energie et l’environnement (PIMENT), a Philippe Lauret, Mathieu David, Emeric Tapaches, Thierry Mara y Malik Mamode por compartir su tiempo y sus conocimientos conmigo y tener la paciencia de hablar en ingl´es. En especial agradecer a Mathieu David su ayuda con los datos del ECMWF y sus consejos y a Thierry Mara por su compa˜n´ıa durante mi estancia. Por otra parte, quisiera agradecer su apoyo a mis compa˜neros del Departamento de Ingenier´ıa El´ectrica y a los miembros del grupo de investigaci´on del IUSIANI. Quisiera destacar sobre todo su generosidad al aceptarme en el equipo y sus consejos a la hora de elaborar esta tesis. Quiero agradecer tambi´en la capacidad de trabajo, ayuda con los c´odigos y aportaci´on de ideas a Brais Pereira y Raquel P´erez. iii iv No puedo olvidarme de los datos de radiaci´on sin los que no podr´ıa haber realizado esta tesis. Quiero agradecer la ayuda al Instituto Tecnol´ogico de Canarias por los datos de radiaci´on de las estaciones de medida de las Islas Canarias, en especial a Antonio Orteg´on por su atenci´on y su trabajo sin condiciones. Agradecer tambi´en al Dr. Mathieu David y al PIMENT por los datos extra´ıdos del modelo num´erico ECMWF y al Dr. Philippe Blanc, Responsable des activit´es de recherche sur l’´evaluation des ressources ´energ´etiques renouvelables of MINES ParisTech / ARMINES por la base de datos satelitales. Agradecer tambi´en el apoyo econ´omico recibido de la C´atedra Endesa Red de la Universidad de Las Palmas de Gran Canaria en su Convocatoria del curso 2014-2015 de ayudas para realizar tesis doctorales y al Vicerrectorado de Internacionalizaci´on y Cooperaci´on de la Universidad de Las Palmas de Gran Canaria, en el marco de la Convocatoria 2014/15 del Programa de apoyo a PFC, TFT y Tesis Doctorales de la ULPGC definidas en el ´ambito de la Cooperaci´on Internacional para el Desarrollo. Por ´ultimo pero no por ello menos importante agradecer el cari˜no de mis hijos Irene y Pablo. De manera muy especial quiero agradecer a Elena por su paciencia, su apoyo sin condiciones a mi trabajo, su ayuda infinita durante estos meses en los que he estado ausente y su compa˜n´ıa. Tambi´en agradecer la ayuda y apoyo recibida de mis padres para llegar hasta aqu´ı. Esta tesis ha sido desarrollada en el marco del siguiente proyecto subvencionado: Integraci´on de nuevas metodolog´ıas en simulaci´on de campos de vientos, radiaci´on solar y calidad del aire. Subvencionado por la Convocatoria 2014 −Proyectos I+D+I −Programa Estatal de Investigaci´on, Desarrollo e Innovaci´on orientada a los retos de la sociedad. Referencia: CTM2014-55014-C3-1-R ´ Indice general Abstract 1 1. Introducci´on 3 2. Fundamentos y datos de radiaci´on solar 15 2.1. Fundamentos de radiaci´on solar . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.1. Naturaleza de la radiaci´on solar . . . . . . . . . . . . . . . . . . . . . 15 2.1.2. Principios del movimiento solar . . . . . . . . . . . . . . . . . . . . . 17 2.1.3. Radiaci´on exoatmosf´erica sobre superficie horizontal . . . . . . . . . 23 2.2. Modelo de Cielo despejado . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3. Datos terrestres de radiaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.1. Estaciones de medida en la isla de Gran Canaria . . . . . . . . . . . 28 2.3.2. Equiposdemedida............................ 32 2.3.3. Tratamiento de los datos de radiaci´on solar . . . . . . . . . . . . . . 33 2.3.4. An´alisis de los datos de radiaci´on solar . . . . . . . . . . . . . . . . . 35 2.4. Adquisici´on de datos del ECMWF . . . . . . . . . . . . . . . . . . . . . . . 43 2.4.1. Datos de radiaci´on solar del ECMWF . . . . . . . . . . . . . . . . . 44 2.4.2. Datos de viento y humedad del ECMWF . . . . . . . . . . . . . . . 46 3. Radiaci´on solar satelital 49 3.1. Datos obtenidos del modelo de radiaci´on Heliosat . . . . . . . . . . . . . . . 50 3.2. An´alisis de los datos de radiaci´on satelitales . . . . . . . . . . . . . . . . . . 52 3.3. An´alisis espacio temporal de los datos . . . . . . . . . . . . . . . . . . . . . 54 3.3.1. Correlaci´on con el ´ındice de cielo despejado . . . . . . . . . . . . . . 56 3.3.2. Correlaci´on con la variaci´on del ´ındice de cielo despejado . . . . . . 63 v vi ´ INDICE GENERAL 4. Modelos de predicci´on de radiaci´on solar 67 4.1. Introducci´on.................................... 67 4.2. Modelos estad´ısticos de referencia . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2.1. Modelo de predicci´on Persistente . . . . . . . . . . . . . . . . . . . . 68 4.2.2. Modelo de predicci´on Smart-Persistence . . . . . . . . . . . . . . . . 69 4.2.3. Modelo de predicci´on Climatol´ogico . . . . . . . . . . . . . . . . . . 69 4.3. Modeloslineales ................................. 70 4.3.1. Modelo lineal Autorregresivo AR . . . . . . . . . . . . . . . . . . . . 71 4.3.2. Modelo lineal autorregresivo de medias m´oviles ARMA . . . . . . . 71 4.3.3. Estudio de la complejidad de los modelos lineales . . . . . . . . . . . 72 4.4. Redes Neuronales Artificales . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.4.1. Introducci´on a las redes neuronales . . . . . . . . . . . . . . . . . . . 76 4.4.2. Aplicaciones................................ 78 4.4.3. Fundamentos te´oricos de las redes neuronales . . . . . . . . . . . . . 79 4.4.4. Laneuronasimple ............................ 83 4.4.5. El Perceptr´on Multicapa . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.4.6. Regla de aprendizaje. Backpropagation . . . . . . . . . . . . . . . . 89 4.4.7. Algoritmo de gradiente conjugado escalado . . . . . . . . . . . . . . 94 4.4.8. T´ecnicas de optimizaci´on de la arquitetura de la red . . . . . . . . . 95 4.5. Redes Neuronales Bayesianas . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.5.1. Aproximaci´on probabilista del aprendizaje en RNA . . . . . . . . . . 97 4.5.2. Optimizaci´on bayesiana de los par´ametros de control . . . . . . . . . 99 4.5.3. Selecci´on bayesiana de la arquitectura de una RNA . . . . . . . . . . 100 4.5.4. T´ecnica bayesiana para la selecci´on autom´atica de entradas relevantes101 5. Aplicaci´on de los modelos de predicci´on 105 5.1. Introducci´on....................................105 5.2. Predicci´on a partir de datos terrestres y del ECMWF . . . . . . . . . . . . 107 5.2.1. Modelos lineales de predicci´on . . . . . . . . . . . . . . . . . . . . . 108 5.2.1.1. Estudio de la complejidad de los modelos lineales . . . . . 108 5.2.1.2. Resultados obtenidos con los modelos lineales . . . . . . . . 115 5.2.2. Modelos de predicci´on basados en Redes Neuronales Artificiales . . . 120 5.2.2.1. Estudio de la complejidad de las RNAs . . . . . . . . . . . 122 5.2.3. Resultados obtenidos y comparaci´on entre los modelos . . . . . . . . 126 ´ INDICE GENERAL vii 5.2.3.1. Resultados anuales y trimestrales . . . . . . . . . . . . . . 130 5.2.3.2. Resultados seg´un el tipo de d´ıa . . . . . . . . . . . . . . . . 140 5.3. Predicci´on a partir de datos satelitales y otras variables . . . . . . . . . . . 153 5.3.1. Selecci´on de los datos de radiaci´on Sat´elite . . . . . . . . . . . . . . 154 5.3.2. Estudio de la predicci´on utilizando otros datos meteorol´ogicos del ECMWF .................................156 5.3.3. Resultados obtenidos y comparaci´on entre los modelos . . . . . . . . 158 5.3.3.1. Resultados anuales y trimestrales . . . . . . . . . . . . . . 159 5.3.3.2. Resultados seg´un el tipo de d´ıa . . . . . . . . . . . . . . . . 168 5.4. Producci´on de energ´ıa el´ectrica . . . . . . . . . . . . . . . . . . . . . . . . . 186 6. Conclusiones y l´ıneas futuras 189 Conclusions 195 Listado de Figuras 199 Listado de Tablas 207 Bibliograf´ıa 211 6CAP´ ITULO 1. INTRODUCCI ´ ON con resultados precisos [Hammer99, Perez10]. Estas im´agenes satelitales se encuentran disponibles con resoluciones espaciales entre 1 km y 5 km y temporales entre 15 min y 30 min para los actuales sat´elites. Los modelos num´ericos NWP se utilizan para realizar predicciones desde las 6 horas hasta los d´ıas de horizonte temporal. Los modelo estad´ısticos de predicci´on se muestran en general eficaces para todos los horizontes temporales de predicci´on. As´ı, se obtienen buenos resultados a muy corto plazo utilizando modelos para series temporales basadas ´unicamente en datos terrestres medidos en la localidad en cuesti´on. Tambi´en se obtienen predicciones fiables de para horizontes diarios de predicci´on utilizando m´etodos estad´ısticos para refinar los resultados previos de un modelo num´erico de predicci´on [Diagne13]. Existe un gran potencial de mejora en las predicciones realizando simulaciones con la combinaci´on de datos diferentes en conjunto con diversos modelos estad´ısticos basados en el aprendizaje. El objetivo de esta tesis es el estudio de distintos m´etodos de predicci´on de la irradiancia solar global horizontal en el horizonte temporal horario, de 1 a 6 horas en adelanto con intervalos horarios. Como ya se ha comentado, seg´un el horizonte temporal de predicci´on existen diversas t´ecnicas en constante desarrollo. En general, los modelos de predicci´on estad´ısticos, como los m´etodos Autorregresivos (AR), los m´etodos Autorregresivos de Media M´oviles (ARMA) y las Redes Neuronales Artificiales (RNAs), son los m´as apropiados para horizontes de predicci´on entre 5 min. y 6 h [Lorenz12]. Durante esta tesis se trabajar´a con estos modelos y datos de entrada provenientes de estaciones de medida terrestres, de im´agenes satelitales y del modelo num´erico ECMWF. Estado del arte En el presente trabajo se analizar´an principalmente diferentes modelos estad´ısticos de predicci´on de radiaci´on solar, as´ı como m´etodos para la mejora de la predicci´on utilizando datos provenientes de otro modelos de predicci´on. En el primer caso nos centraremos en los trabajos realizados con modelos lineales y RNAs en la predicci´on de radiaci´on solar. Los datos ex´ogenos utilizados ser´an los obtenidos por medio de m´etodos num´ericos de predicci´on y datos de radiaci´on basados en im´agenes satelitales. Los modelos lineales de predicci´on (AR y ARMA) has sido ampliamente utilizadas desde los a˜nos setenta en el campo de la radiaci´on solar. Los modelos de series lineales se han utilizado para describir su comportamiento, para generar series sint´eticas de valores o para definir A˜nos Meteorol´ogicos T´ıpicos [Mazorra10], Typical Meteorological Year (TMY) en ingl´es, y para realizar predicciones de radiaci´on solar, a escala tanto horaria como diaria. J. Boland, [Boland95, Boland08], realiz´o una descripci´on de la radiaci´on solar diaria y horaria, un estudio del modelado de series temporales de radiaci´on solar y un m´etodo de estimaci´on de la radiaci´on solar difusa en Australia a partir de modelos Autorregresivos (AR), de medias m´oviles (MA) y Autorregresivos de medias m´oviles (ARMA). Estos trabajos se basan en la Metodolog´ıa de Box-Jenkins, [Box98], que describe un proceso iterativo para identificar el modelo ´optimo y luego utilizarlo en las predicciones. G. Reikard [Reikard09] se bas´o en modelos Autoregresivos Integrados de medias m´oviles (ARIMA) para seis estaciones en los Estados Unidos. P. Bacher [Bacher09] por su parte, evalu´o para estaciones en Dinamarca el comportamiento de modelos Autoregresivos (AR) simples, 7 con datos pasados de la serie temporal y a˜nadiendo datos ex´ogenos (ARX) provenientes de modelos n´umericos de predicci´on meteorol´ogica, (NWP) por sus siglas en ingl´es. S. SAFI, [Safi02], estudi´o distintos modelos de medias m´oviles (MA) para localidades en Marruecos. Aguiar & Collares-Pereira realizaron trabajos con modelos ARMA a partir de datos del ´ındice de claridad horario [CP89]. En los ´ultimos a˜nos nos encontramos con numerosas metodolog´ıas basadas en modelos estad´ısticos autorregresivos (AR) a partir de datos de radiaci´on solar, tanto de estaciones de medida como provenientes de sat´elites [Dambreville14b, Dambreville14a, Zagouras15]. Por otro lado, en [Lauret12] se propone un m´etodo h´ıbrido entre redes neuronales y ARMA basados en t´ecnicas de decisi´on Bayesianas, mientras que en [David14] se propone el uso de un modelo ARMA recurrente para la predicci´on de radiaci´on solar horaria en una isla de clima Tropical, Isla Reuni´on. En esta tesis las Redes Neuronales Artificiales utilizadas para predecir la radiaci´on solar es el Perceptr´on Multicapa. Desde que en 1957 Frank Rosenblatt [Rosenblatt58] desarroll´o este tipo de redes y en 1986 Rumelhart, Hinton y Williams presentan su trabajo en el que se desarrolla un algoritmo de aprendizaje conocido como retropropagaci´on (backpropagation) [Rumelhart86] para redes neuronales multicapa, el n´umero de trabajos sobre RNA se han multiplicado. Se han desarrollado un gran n´umero de aportaciones en los m´etodos de aprendizaje y tipo de estructuras. Las RNAs es una herramienta capaz de reproducir una relaci´on no lineal entre un conjunto de datos de entrada y otro de salida [Bishop95], lo cual lo hace una herramienta muy atractiva. Las RNAs se han utilizado con ´exito para realizar predicciones con las series temporales de radiaci´on solar en forma del ´ındice de cielo despejado. Estos modelos estad´ısticos pueden trabajar ´unicamente bas´andose en datos hist´oricos de radiaci´on solar. Se han realizado trabajos de predicci´on con RNAs para datos de irradiancia solar desde horizontes temporales horarios, como en los trabajos de [Hontoria02, Lauret06b, Mellit10, Inman13], hasta para predecir la irradiaci´on solar con 24 horas de antelaci´on [Bosch08, Rehman08]. Los modelos estad´ısticos basados en RNAs permiten a˜nadir otro tipo de variables como datos de entrada. As´ı, se puede trabajar con RNAs combinando datos de radiaci´on hist´oricos y otras variables meteorol´ogicas. De esta manera podemos encontrar a Rehman [Rehman08], utilizando datos de temperatura y humedad relativa terrestre para predecir la irradiaci´on diaria. Kemmoku [Kemmoku99] predice la irradiaci´on diaria con una aplicaci´on en serie de varias RNAs partiendo de datos de presi´on atmosf´erica y otros datos meteorol´ogicos. Sfetsos & Coonick [Sfetsos00] introducen datos de temperatura, velocidad del viento y presi´on, adem´as de los datos de radiaci´on solar, para predecir valores de irradiancia global horaria. De la misma manera, se pueden encontrar diversos trabajos de predicci´on con diferentes combinaciones de datos meteorol´ogicos como la longitud del d´ıa (horas de sol), temperatura media, humedad relativa, latitud y longitud para obtener datos tanto de irradiancia horaria como de irradiaci´on diaria [Mohandes98, Ghanbarzadeh09, Mellit10]. Adem´as de las RNAs existen otros tipos de modelos estad´ısticos no lineales ampliamente utilizados en la predicci´on de la radiaci´on solar. Estos modelos tambi´en se basan en t´ecnicas de aprendizaje por lo que requieren de una base datos hist´oricos. En esta tesis se utilizar´an las RNAs basadas en las t´ecnicas probabilistas bayesianas seg´un los trabajos de [MacKay03]. Estas t´ecnicas permiten entre otros aspectos mejorar el proceso de aprendizaje, adem´as de estudiar la complejidad el modelo a definir [Lauret08, BS10]. 8CAP´ ITULO 1. INTRODUCCI ´ ON En los ´ultimos a˜nos se han desarrollados otras t´ecnicas de Machine Learning con las que se est´an obteniendo buenos resultados de predicci´on. De esta manera encontramos trabajos de predicci´on basados en Support Vector Machine [Zeng13, FJ13, Wolff13] en los que se exponen los resultados obtenidos con esta t´ecnica, y trabajos basados en los Gaussian Process [Sun14, Lauret15]. Como variante a las RNAs convencionales se han desarrollado en los ´ultimos a˜nos redes neuronales con funciones de activaci´on basadas en la funci´on de ondas (wavelet function) [Mellit10, Cao08]. Estas redes se conocen como Wavelet Neural Networks (WNNs) o Wavelet Networks (WNs). Los modelos num´ericos de predicci´on son adecuados para predecir distintas variables atmosf´ericas hasta un horizonte temporal de unos 15 d´ıas. Las variaciones en el estado de la atm´osfera se modelan en base a unas ecuaciones diferenciales que describen los procesos f´ısicos. Los modelos NWP globales se encuentran actualmente operados por 15 diferentes agencias meteorol´ogicas mundiales. Como ejemplo tenemos el Global Forecaste System (GFS) utilizado por la US National Oceanic and Atmospheric Administration (NOAA) y el Integrated Forecast System (IFS) operado por la European Centre for Medium-Range Weather Forecasts(ECMWF). Por otro lado los modelos de mesoescala se encuentran disponibles ´unicamente para algunas zonas del globo terrestre pero ofrecen una resoluci´on espacial mayor que los globales. Entre estos modelos encontramos el MM5 desarrollado por la Pennsylvania State University y el National Centre for Atmospheric Research (NCAR) o el modelo WRF dise˜nado como un modelo. En los ´ultimos a˜nos se han realizado estudios de comparaci´on con las predicciones realizadas por estos modelos locales [Heinemann06b, Perez11]. La precisi´on de estos modelos de predicci´on var´ıa seg´un la escala temporal utilizada y zona geogr´afica en la que se trabaje. Heinemann [Heinemann06b] muestra que se pueden obtener datos de radiaci´on para cielos despejados sin pr´acticamente ninguna desviaci´on. Una comparaci´on de los resultados obtenidos por estos modelos para estaciones de Estados Unidos, Canad´a y Europa se describen en [Perez07, Perez10, Perez13] mostrando errores en la predicci´on de radiaci´on horaria de 38 % rRMSE. En Europa, se han obtenido resultados del orden del 40 % de error para estaciones del centro de Europa y del 30 % de error para estaciones en Espa˜na. Adem´as han sido analizados los resultados con respecto a diferentes propiedades relevantes para su aplicaci´on en la producci´on fotovoltaica [Lorenz09b, Lorenz09a, Lorenz11]. Las series temporales de radiaci´on solar tienen una parte determinista diaria y anual, pero tambi´en existe una componente aleatoria debido entre otros aspectos a la presencia de nubes en el cielo. As´ı, determinar la nubosidad de un determinado lugar podr´a ofrecer una informaci´on valiosa para predecir la radiaci´on. Para horizontes de predicci´on superiores a las horas, el cambio de la nubosidad est´a fuertemente influenciado por el movimiento de las nubes. Las im´agenes satelitales y las im´agenes hemiesf´ericas del cielo proporcionan la posibilidad de predecir la presencia de nubes, extrapolando el movimiento de las mismas en las horas anteriores al horizonte temporal deseado. Los modelos de predicci´on basados en im´agenes satelitales e im´agenes hemiesf´ericas del cielo detectan el movimiento de las nubes utilizando t´ecnicas de seguimiento de los vectores del movimiento Cloud Motion Vectors [Lorenz12]. Las predicciones a corto plazo basadas en im´agenes del cielo es un campo relativamente nuevo con diversos trabajos para horizontes temporales intrahorarios principalmente [Chow11, Urquhart13]. 9 Las predicciones basadas en los vectores de movimiento de las nubes obtenidas a partir de im´agenes satelitales logran unos resultados del orden del 17 %rRMSE para predecir el ´ındice de nubosidad en horizontes temporales de 30 minutos y un 30 %rRMSE hasta 2 horas [Hammer99]. Por otro lado en [Lorenz09a] se realiz´o una comparaci´on de varios m´etodos basados en los vectores de movimiento de las nubes a partir de las im´agenes del Meteosat para predecir irradiancia solar con horizontes temporales superiores a una hora. Mientras que [Perez10] muestra los resultados de las predicciones de irradiancia basados en las im´agenes del Geostationary Operational Enviromental Satellite (GOES). Los datos satelitales utilizados en esta tesis se obtuvieron de la base de datos del Helioclim-3, en particular, de la versi´on 5 (HC3v5). Los datos de radiaci´on se estiman a partir de las im´agenes satelitales del Meteosat con el m´etodo del Heliosat-2 [Rigollier04, Blanc11b]. Esta versi´on le proporcion´o al Helioclim una mejor resoluci´on temporal (15 minutos) y espacial (3 km nadir). Las im´agenes se obtienen en tiempo real en la estaci´on de recepci´on del METEOSAT y se realizan los c´alculos cada 15 minutos. La base de datos HelioClim estima la radiaci´on a cielo despejado con el modelo de McClear, que usa los datos de AOD, Ozono y vapor de agua del proyecto MACC [Lefevre13]. Con las ´ultimas actualizaciones se ha mejorado el error obtenido por el modelo Heliosat-2 tanto para d´ıas despejados como para d´ıas nubosos [Eissa15] y se ofrece, adem´as de la estimaci´on de la radiaci´on ocurrida, una predicci´on para las pr´oximas horas [Thomas15]. En los ´ultimos a˜nos han surgido diversos estudios basados en la utilizaci´on conjunta de datos hist´oricos de radiaci´on medidos en estaciones terrestres, datos obtenidos por alg´un modelo num´erico y datos de radiaci´on obtenidos a partir de datos satelitales. Estos datos se utilizan como variables de entrada en diversos modelos estad´ısticos de predicci´on. As´ı, con los datos provenientes de im´agenes satelitales nos encontramos con modelos autorregresivos (AR) [Dambreville14b, Dambreville14a, Zagouras15], redes neuronales artificiales [Marquez13, MA15] y algoritmos gen´eticos para elegir las informaci´on relevante de toda la base de datos del sat´elite [Zagouras15]. La elecci´on de los p´ıxeles de inter´es para mejorar la predicci´on de los datos terrestres es unos de los puntos m´as importantes del trabajo. En este caso nos encontramos con [Dambreville14a] que propone utilizar las correlaciones entre las series temporales de la variaci´on del´ındice de cielo despejado cada 15 minutos, tanto de los datos terrestres como satelitales. Por otra parte [Zagouras15, MA15] proponen utilizar las correlaciones entre las series temporales del ´ındice de cielo despejado horarias. En ambos casos las correlaciones se realizan con los datos satelitales desfasados con respecto a los terrestres. Por otro lado, para tiempos de predicci´on superiores a 6 h los m´etodos num´ericos de predicci´on (NWP) post-procesados con un modelo estad´ıstico y datos terrestres de radiaci´on muestran los mejores resultados [Diagne13]. Objetivos y metodolog´ıa Con los antecedentes descritos, el conocimiento de la radiaci´on se revela fundamental para diversas actividades relacionadas con la producci´on y la gesti´on de la energ´ıa el´ectrica, sobre todo en regiones insulares. Por lo tanto, contar con modelos de predicci´on de radiaci´on solar fiables para ser utilizados en diversos campos resulta de clara aplicaci´on. F. D´ıaz [DR13] desarrolla un modelo basado en las consideraciones geom´etricas de la isla 10 CAP´ ITULO 1. INTRODUCCI ´ ON de Gran Canaria para establecer los niveles de radiaci´on solar en cada punto. Los datos de radiaci´on solar utilizados se basaban en datos hist´oricos y datos obtenidos a partir de un modelo num´erico de predicci´on meteorol´ogica (MM5). En esta tesis, para continuar con la l´ınea de investigaci´on, se ha decidido mejorar las predicciones de radiaci´on utilizadas por el modelo geom´etrico descrito por F. D´ıaz [DR13], adem´as de conseguir un horizonte temporal de predicci´on horario. Para cumplir este objetivo se han trabajado los siguientes aspectos: Un modelo de radiaci´on a cielo despejado que contemple las condiciones f´ısicas de atm´osfera y la geometr´ıa solar. Recopilaci´on de los datos de radiaci´on terrestres de las estaciones de medida en la isla de Gran Canaria y tratamiento de los mismos. Se ha realizado un filtrado de los datos seg´un el modelo SERI QC y un ´angulo cenital concreto. Las medias horarias se han calculado a partir de los datos recogidos cada minuto ´unicamente en aquellas horas en las que se dispon´ıa del 50 % de los datos. An´alisis de los datos de radiaci´on solar realizando una distribuci´on de los tipos de d´ıas en funci´on de la media diaria del ´ındice de cielo despejado y la desviaci´on est´andar de la variaci´on del mismo. Estudio de la cuadr´ıcula de datos de radiaci´on solar obtenidos alrededor de la isla de Gran Canaria a partir de im´agenes satelitales, en concreto del HelioClim-3. Se ha desarrollado una herramienta que permite visualizar los datos de irradiancia solar y el´ındice de cielo despejado cada 15 minutos durante todo el a˜no de medida estudiado. De esta manera se permite visualizar las variaciones intradiarias de la radiaci´on en la geograf´ıa insular. C´alculo de las correlaciones entre los´ındices de cielo despejado de los datos terrestres de las estaciones de medida y cada uno de los p´ıxeles de la cuadr´ıcula de datos satelitales obtenida. Las correlaciones se han calculado con los datos satelitales desfasados con respecto a los datos terrestres para poder estudiar la relaci´on entre las estaciones de la isla y su alrededor en distintos instantes. Desarrollo de una herramienta que permite visualizar los resultados de estas correlaciones para cada estaci´on y para cada desfase temporal entre ambas series. Estos resultados se pueden observar para todo el conjunto de datos anuales y seg´un la estaci´on del a˜no elegida (invierno, primavera, verano y oto˜no). La correlaci´on nos permite calcular la relaci´on entre cada p´ıxel y la estaci´on de medida, con lo que nos dar´a una pista de los datos m´as significativos para obtener mejores resultados de predicci´on. Estudio de los modelos lineales de predicci´on AR y ARMA. Se calcula la complejidad de ambos modelos a partir las series temporales del ´ındice de cielo despejado y se decide el modelo lineal ´optimo para realizar las predicciones de radiaci´on solar. Estudio del modelo ´optimo de Redes Neuronales Artificiales para realizar las predicciones de radiaci´on solar a partir de datos de radiaci´on terrestres y datos 11 ex´ogenos. Se decide el n´umero de neuronas ocultas de la capa intermedia y el n´umero de entradas relevantes para cada caso. Elecci´on de los p´ıxeles de la cuadr´ıcula de datos satelitales que consiguen mejores resultados de predicci´on en conjunto con los datos terrestres como entradas de las RNAs. Elecci´on de la altitud atmosf´erica a la que extraer los datos de velocidad del viento y humedad relativa del modelo num´erico de predicci´on. Se estudian diferentes altitudes y se comparan los resultados de predicci´on de radiaci´on solar obtenidos. Aplicaci´on de los diferentes modelos estad´ısticos para la predicci´on de la radiaci´on en las diferentes estaciones de la isla de Gran Canaria. Los modelos se han estudiado utilizando datos terrestres de radiaci´on solar y datos ex´ogenos como los datos satelitales o los datos del modelo num´erico de predicci´on elegido. Para realizar el tratamiento de datos de radiaci´on solar se ha partido del c´alculo de modelo de cielo despejado siguiendo los trabajos de Bird & Hulstrom [Bird81] utilizando los valores de aerosoles AOD500 nm, AOD380 nm y el contendio en vapor de agua de la columna vertical obtenidos de la red AERONET [AERONET14, Holben98]. Una vez calculado el ´ındice de cielo despejado se ha realizado un filtrado de los datos disponibles en cada estaci´on siguiendo el m´etodo SERI-QC [Maxwell93, Younes05] para bases de datos de radiaci´on global horizontal ´unicamente. Por otro lado, para trabajar en los modelos estad´ısticos de predicci´on con series temporales continuas se filtran las horas nocturnas. Este filtro se basa en el ´angulo cenital de cada medida. El criterio seguido en esta tesis establece la frontera de datos v´alidos en 80o, a partir del cual se considera que los datos son horas nocturnas. El an´alisis de la climatolog´ıa de la isla a partir de los datos horarios de radiaci´on solar de las seis estaciones de medida se realiza a partir de la caracterizaci´on del comportamiento meteorol´ogico de cada estaci´on. Para caracterizar cada estaci´on, la base de datos de radiaci´on solar horaria se ha clasificado seg´un una distribuci´on de los tipos de d´ıas en funci´on de la media diaria del ´ındice de cielo despejado y la desviaci´on est´andar de la variaci´on del mismo en los a˜nos estudiados [Dambreville14b, Dambreville14a]. Los d´ıas se dividen en nueve tipos, donde la media se divide en d´ıas que presentan medias diarias bajas, por lo que se consideran d´ıas nubosos, hasta los d´ıas que presentan media altas de radiaci´on, por lo que representan los d´ıas despejados. La variabilidad se divide en d´ıas que presentan una variabilidad bastante baja, por lo que se consideran d´ıas con los valores de radiaci´on estables, y d´ıas que presentan variabilidades altas. Los datos satelitales utilizados en esta tesis se obtuvieron de la base de datos del Helioclim-3, en particular, de la versi´on 5 (HC3v5). El an´alisis de los datos satelitales ofrecidos por la cuadr´ıcula obtenida alrededor de la isla nos permite confirmar los efectos de la meteorolog´ıa de la zona en la radiaci´on solar. Con la herramienta de visualizaci´on de los datos desarrollada se puede analizar la variaci´on de la radiaci´on cada 15 minutos. Tambi´en se ha realizado la correlaci´on entre las dos bases de datos (datos satelitales y datos terrestres de cada estaci´on) utilizando los ´ındices de cielo despejado para estudiar la relaci´on entre ambos [Dambreville14b]. Para poder evaluar la correlaci´on entre las series 12 CAP´ ITULO 1. INTRODUCCI ´ ON temporales en diferentes momentos se establecieron desfases temporales entre ambas series. En cada estaci´on se dispone de las correlaciones entre cada p´ıxel de los datos sat´elite desfasado hasta 3 horas y el dato terrestre correspondiente a la propia estaci´on sin desfasar. De esta manera se pretende conocer y valorar la relaci´on entre cada punto geogr´afico de los alrededores y la estaci´on en la que se desear realizar la predicci´on. Los modelos estad´ısticos utilizados en esta tesis establecen una relaci´on entre unos datos de entrada y una variable de salida esperada. Para obtener la relaci´on entre ambas series temporales es necesario un entrenamiento a partir de un conjunto de datos hist´oricos recogidos por las estaciones de medida. Durante este entrenamiento se debe estudiar adem´as la complejidad de dichos modelos estad´ısticos. En el caso de los modelos lineales, esto conlleva decidir el orden de los mismos, es decir, el n´umero de par´ametros relevantes para realizar la predicci´on. En esta tesis se ha estudiado la complejidad de los modelos lineales siguiendo el m´etodo descrito por J. Boland [Boland95, Boland08], utilizando las Funciones de Autocorrelaci´on y Autocorrelaci´on Parcial, ACF y PACF, y los criterios de decisi´on bayesiana, BIC. En cuanto a las RNAs, la complejidad del modelo cosiste en estudiar el n´umero de variables de entrada y n´umero de neuronas ocultas relevantes para la predicci´on. En esta tesis las teor´ıas que se van a utilizar est´an basadas en la interpretaci´on probabilista del algoritmo de retropropagaci´on para el perceptr´on multicapa realizado por Mackay [MacKay03, MacKay92]. Como ya se ha comentado, las RNAs se utilizan partiendo de diferentes conjuntos de datos de entradas. Para cada simulaci´on que se realice durante el trabajo se ha estudiado la complejidad del modelo seg´un las t´ecnicas probabilistas bayesianas. En los ´ultimos a˜nos se han realizado diversos trabajos de predicci´on de radiaci´on solar partiendo de datos terrestres en conjunto con datos satelitales o datos de un modelo num´erico de predicci´on [Marquez13, Dambreville14b, Diagne14, Zagouras15]. Uno de los aspectos m´as importantes destacados en estos trabajos es la elecci´on de los datos ex´ogenos relevantes para la predicci´on. En este caso, se han utilizado los datos satelitales con mayores valores de correlaci´on con respecto a la estaci´on de medida para mejorar los resultados obtenidos con los datos terrestres. Utilizando la herramienta desarrollada para visualizar las correlaciones calculadas se realizaron diversas simulaciones cambiando la elecci´on de p´ıxeles. En cuanto a los datos de viento y humedad relativa del modelo num´erico ECMWF se ha trabajado con el m´odulo y direcci´on de la velocidad y la humedad relativa a una altitud atmosf´erica determinada. Para determinar la altitud a la que extraer los datos referidos se utilizaron tres diferentes criterios en cada estaci´on [Lave13, DA96, Badosa15]. Por ´ultimo se ha realizado la aplicaci´on de los modelos estad´ısticos de predicci´on a los datos de Gran Canaria para establecer el modelo que mejor se ajusta a las condiciones clim´aticas. Los resultados obtenidos con todos los modelos se han estudiado teniendo en cuenta los valores anuales de radiaci´on, pero tambi´en se han estudiado los resultados seg´un la estaci´on del a˜no y el tipo de d´ıa (siguiendo la distribuci´on de nueve d´ıas obtenida). De esta manera se pretende conocer el comportamiento general de los modelos y su fiabilidad para distintas condiciones meteorol´ogicas. 13 Contenido de la tesis En el Cap´ıtulo 2 se presentan los datos terrestres de radiaci´on solar para cada una de las estaciones de la isla de Gran Canaria. Se exponen los resultados del tratamiento y an´alisis de los datos, estableciendo una primera relaci´on entre las estaciones de medida seg´un su ubicaci´on. Tambi´en se describen los datos de radiaci´on solar, viento y humedad relativa obtenidos del modelo num´erico de predicci´on ECMWF. Los datos de radiaci´on solar obtenidos a partir de im´agenes satelitales se presentan en el Cap´ıtulo 3 de esta tesis. Se describen las herramientas desarrolladas para poder visualizar los datos del´ındice de cielo despejado y los valores obtenidos en las correlaciones entre cada p´ıxel satelital y el dato terrestre. Adem´as se presenta el an´alisis de los datos observados, tanto anualmente como por cada estaci´on meteorol´ogica del a˜no. Los modelos estad´ısticos utilizados durante la tesis se exponen en el Cap´ıtulo 4. Es en este cap´ıtulo donde se explican los modelos na¨ıve (Persistence, Smart-Persistence), el modelo Climatol´ogico, los modelos lineales (AR y ARMA) y los modelos basados en Redes Neuronales Artificiales (RNAs). Adem´as se describe la metodolog´ıa seguida para establecer una relaci´on entre los datos de entrada y la salida deseada en cada caso. En el caso de los modelos lineales y las RNAs se exponen las t´ecnicas utilizadas para establecer la complejidad de los modelos. Por ´ultimo, en el Cap´ıtulo 5 se presentan los resultados obtenidos para las estaciones de medida de la isla con cada uno de los modelos. Las simulaciones se han dividido en dos grandes grupos seg´un el n´umero de a˜nos de medida con los que se pueda trabajar. En un primer momento, se ha trabajado con un a˜no de medida para el entrenamiento de los modelos y otro a˜no de medida para comprobar el ajuste de los modelos entrenados. En este grupo se han simulado los modelos estad´ısticos con datos terrestres y datos de radiaci´on solar del ECMWF como variables de entrada. Los datos de radiaci´on solar satelitales y los datos de viento y humedad relativa del ECMWF est´an ´unicamente disponibles para el a˜no 2005. As´ı, en un segundo grupo se han realizado las simulaciones de los modelos estad´ısticos utilizando ´unicamente un a˜no de medida. En este caso se debe dividir el a˜no 2005 en dos conjuntos de datos, uno para realizar el entrenamiento y otro para el test. Finalmente, en el Cap´ıtulo 6 se presentan las principales conclusiones y las l´ıneas futuras de investigaci´on del presente trabajo. Cap´ıtulo 2 Fundamentos y datos de radiaci´on solar En este cap´ıtulo se realizar´a una breve introducci´on sobre la naturaleza de la radiaci´on solar que llega a la atm´osfera terrestre, de manera que se pueda entender la forma de proceder con los datos disponibles. Se explicar´an los distintos factores que influyen en la radiaci´on solar disponible en la atm´osfera terrestre, la descomposici´on de la misma al atravesar la atm´osfera y los ´angulos de incidencia seg´un la localizaci´on y la ´epoca del a˜no. Para trabajar con los modelos de predicci´on de radiaci´on solar elegidos en esta tesis se utilizar´an las series temporales del ´ındice de cielo despejado, por lo que se deber´a calcular el modelo de cielo despejado para cada estaci´on de medida. Adem´as tambi´en se presentar´an los datos disponibles para la isla de Gran Canaria en las diferentes estaciones de medida del Instituto Tecnol´ogico de Canarias, as´ı como los aparatos de medida utilizados en dichas estaciones. Cuando se trabaja con una base de datos tan amplia, se considera b´asico contar con unos datos fiables. En este cap´ıtulo tambi´en se describir´a el proceso de tratamiento de datos seguido en esta tesis y el posterior an´alisis de los mismos para realizar una caracterizaci´on del clima local. Por ´ultimo, se presentar´an los datos obtenidos del modelo num´erico de predicci´on solar, en concreto del European Centre for Medium-Range Weather Forecasts, ECMWF. 2.1. Fundamentos de radiaci´on solar 2.1.1. Naturaleza de la radiaci´on solar El sol es una esfera de gas caliente con un di´ametro de 1,39xxx109m y est´a, de media, a 1.5x1011 m. de la Tierra. Visto desde la Tierra, el Sol gira sobre su eje una vez cada cuatro semanas, aunque en realidad no gira como una masa s´olida. El ecuador gira en 27 d´ıas aproximadamente y las regiones polares se toman unos 30 d´ıas. Tiene una temperatura efectiva de 5777 oK (su n´ucleo puede llegar a los 15 millones oK), y una densidad estimada de 100 veces la del agua. El Sol es de hecho un gran reactor de fusi´on constituido por gases retenidos por fuerzas gravitacionales. 15 22 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR Latitud Φ, si se tiene en cuenta que la vertical de un lugar espec´ıfico de la tierra intersecta a la esfera celeste en dos puntos llamados cenit (Norte) y nadir (Sur). El complementario del ´angulo que forma esta recta con el eje polar es la latitud (posici´on angular norte o sur respecto del ecuador). Se denomina meridiano al c´ırculo m´aximo de la esfera terrestre que contiene los polos, al cenit y al nadir; −90◦<Φ<90◦. Se considera positiva en el hemisferio Norte y negativa en el hemisferio Sur. Declinaci´on δ, la posici´on angular del sol cuando ´este se encuentra en el meridiano del lugar (mediod´ıa solar) con respecto al plano del ecuador y representa el ´angulo entre el ecuador terrestre y la l´ınea que une los centros de la Tierra y el Sol; oscila a lo largo del a˜no entre −23.45◦< δ < 23.45◦. Acimut solar γs, ´angulo desviaci´on de la proyecci´on en el plano horizontal del meridiano del Sol con el meridiano del lugar (el Sur en el hemisferio norte y el Norte en el hemisferio sur). En el Este es negativo y positivo en el Oeste; −180◦< γs<180◦. Para latitudes, norte y sur, entre 23.45oy 66.45o, el acimut solar va a estar entre -90oy 90opara d´ıas de menos de 12 horas de duraci´on; si los d´ıas duran m´as entre la salida y la puesta de sol el acimut solar tendr´a valores mayores de 90oy -90o. El ´angulo azimut solar marca el recorrido solar a lo largo del d´ıa. ´ Angulo de incidencia θ, es el ´angulo entre la radiaci´on directa sobre la superficie y la normal a esa superficie. ´ Angulo cenital θzs, es el ´angulo entre la vertical y la l´ınea hacia el Sol, es decir, el ´angulo de incidencia de la radiaci´on directa sobre una superficie horizontal. La relaci´on entre la cantidad de atm´osfera que tiene que atravesar la radiaci´on directa y la que tendr´ıa que atravesar si el Sol estuviera en su cenit se denomina cantidad de aire m(air mass). Al nivel del mar m= 1 si el Sol est´a en su c´enit, y m= 2 si tenemos un ´angulo cenital de 60o. La parte de la radiaci´on exoatmosf´erica que llega a la superficie terrestre depender´a de este valor. ´ Angulo de altitud solar o elevaci´on h0, ´angulo entre la horizontal y la l´ınea hacia el sol, complementario del anterior. La elevaci´on es el ´angulo con el que vemos al sol si miramos en su direcci´on, tomando como origen la horizontal (suelo). Definiendo la latitud como positiva en el hemisferio Norte y negativa en el hemisferio Sur, se define el ´angulo cenital para una localidad y una ´epoca del a˜no seg´un la expresi´on (2.6). cos(θzs) = sin(h0) = sin δsin Φ + cos δcos Φ cos ω(2.6) El ´angulo horario ωse define como el desplazamiento angular, partiendo de Este a Oeste, del meridiano local seg´un la rotaci´on terrestre. Este valor se considera cero al mediod´ıa, negativo durante la ma˜nana y positivo por la tarde. El movimiento de rotaci´on de la Tierra se considera a una velocidad de 15opor hora. La expresi´on utilizada para calcular el angulo horario es (2.7), 2.1. FUNDAMENTOS DE RADIACI ´ ON SOLAR 23 Figura 2.8: Sistema de referencia sobre la superficie terrestre. Fuente [DR13] ω(◦) = 15((Hr −12) + (Lo −15T z)4 60 +ET/60) (2.7) donde Lo es la longitud local, T z es el tiempo de adelanto respecto al meridiano de Greenwich (meridiano origen del ´angulo horario) del meridiano correspondiente al huso horario local y ET es la ecuaci´on del tiempo que var´ıa de un d´ıa a otro seg´un la ecuaci´on (2.8). ET(min) = (0.000075 + 0.001868 cos(τ)−0.032077 sin(τ)− −0.014615 cos(2τ)−0.040849 sin(2τ))229.18 (2.8) El ´angulo horario es el tiempo basado en el movimiento angular aparente del sol sobre el cielo, con el mediod´ıa solar el sol cruza el meridiano local y la ecuaci´on del tiempo tiene en cuenta las perturbaciones en la rotaci´on de la tierra que afecta a la hora en la que el sol cruza el meridiano. 2.1.3. Radiaci´on exoatmosf´erica sobre superficie horizontal La radiaci´on exoatmosf´erica sobre un plano normal a la radiaci´on directa proveniente del Sol para cualquier d´ıa del a˜no, se puede calcular a partir de una expresi´on sencilla teniendo en cuenta la ligeras variaciones de la distancia entre el sol y la tierra. I0n=Isc0(2.9) 24 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR En esta ecuaci´on, como ya se explicado anteriormente, Isc representa la energ´ıa procedente del Sol por unidad de tiempo, recibida sobre una superficie unitaria perpendicular a la direcci´on de propagaci´on de la radiaci´on fuera de la atm´osfera. Esta energ´ıa se conoce como constante solar y se toma un valor constante, Isc = 1367W/m2. Por otra parte, 0es la variaci´on de la distancia entre el Sol y la Tierra lo largo del a˜no y se calcula seg´un la ecuaci´on (2.3). Teniendo en cuenta la relaci´on entre la radiaci´on sobre una superficie inclinada y una horizontal al plano terrestre, se puede expresar la radiaci´on exoatmosf´erica sobre una superficie horizontal en W/m2como (2.10). I0=Isc0cos(θzs) (2.10) Sustituyendo la ecuaci´on del ´angulo cenital (2.6) en la ecuaci´on (2.10) se obtiene la expresi´on final de la radiaci´on exoatmosf´erica. I0=Isc0(sin δsin Φ + cos δcos Φ cos ω) (2.11) A menudo es necesario para los c´alculos de irradiaci´on solar diaria tener la integral diaria de la radiaci´on exoatmosf´erica sobre una superficie horizontal, H0.´ Esta se obtiene integrando la ecuaci´on anterior en el periodo desde la salida hasta la puesta de sol, donde ωses el ´angulo de la puesta de sol. H0=24 ·3600 πIsc0(πωs 180 sin δsin Φ + cos δcos Φ cos ωs) (2.12) En este caso la radiaci´on exoatmosf´erica diaria viene definida en J/m2dia. 2.2. Modelo de Cielo despejado Las series temporales de radiaci´on solar no se consideran series estacionarias, ya que est´an influenciadas por la variabilidad diaria y anual del movimiento de la Tierra. Cuando se trabaja con modelos estad´ısticos de predicci´on solar, se debe partir de series estacionarias. Las series estacionarias son aquellas que no muestran una tendencia estacional y tienen una varianza constante, presentando una auto correlaci´on constante en todo el intervalo temporal de la misma [Chatfield13]. La radiaci´on solar incidente en una superficie horizontal en una localidad concreta para un tiempo espec´ıfico depende del ´angulo cenital. Para trabajar con modelos estad´ısticos se considera adecuado tratar las influencias deterministas dependientes de la geometr´ıa solar por separado de las influencias no deterministas generadas por fen´omenos atmosf´ericos [Diagne13]. As´ı, para conseguir transformar las series temporales de radiaci´on en series estacionarias se han introducido dos variables, el ´ındice de claridad ky el ´ındice de cielo despejado k∗. El ´ındice de claridad se considera como el resultado de dividir la radiaci´on solar en una superficie horizontal Ientre la radiaci´on exoatmosf´erica en la misma localidad para una 2.2. MODELO DE CIELO DESPEJADO 25 superficie horizontal I0, seg´un la ecuaci´on (2.13). El ´ındice de claridad consigue extraer la influencia determinista de la geometr´ıa solar, ya que la influencia del ´angulo cenital se encuentra incluido en el c´alculo de la radiaci´on exoatmosf´erica. Este ´ındice se utiliza ampliamente para tratar la influencia determinista. k=I I0 (2.13) La segunda variable introducida, el ´ındice de cielo despejado, es ampliamente utilizada en la bibliograf´ıa como m´etodo para evitar el car´acter estacional de las series temporales de radiaci´on solar horaria. En esta tesis se ha optado por utilizar las series temporales del ´ındice de cielo despejado como variable en los modelos de predicci´on. De esta manera se consiguen series temporales estacionarias mediante la siguiente f´ormula (2.14), k∗=I Ics (2.14) donde Ies la radiaci´on solar horaria sobre una superficie horizontal medida, por ejemplo, en una estaci´on de medida. Mientras que Ics es la radiaci´on solar calculada para un d´ıa a cielo despejado. Los modelos de radiaci´on a cielo despejado, clear sky models, estiman la radiaci´on que incide en cualquier superficie en cualquier instante de tiempo considerando unas condiciones de cielo limpio. Algunos de los modelos de cielo despejado consultados en la bibliograf´ıa se basan en diferentes variables clim´aticas como datos de entrada [Reno12, Younes07]. En la isla de Gran Canaria se han probado diversos modelos de cielo despejado, mostrando resultados satisfactorios y similares entre s´ı [DR13, BM01, ˇ S´uri04, Hofierka02, Bird81]. El modelo que se utilizar´a en esta tesis, por su precisi´on y sencillez de aplicaci´on, es el modelo de Bird & Hulstrom [Bird81]. Este modelo es ampliamente conocido y utilizado porque proporciona excelentes resultados a partir ´unicamente de una serie de datos meteorol´ogicos [Badescu13]. Las variables meteorol´ogicas que utiliza el modelo de cielo despejado son los espesores ´opticos de aerosol para longitudes de onda de 500 nm y 380 nm (Aerosol Optical Depths AOD500 nm and AOD380 nm), y el contenido de vapor de agua y ozono en la columna vertical de la atm´osfera. AODs miden la atenuaci´on de la radiaci´on s´olar como resultado de la dispersi´on y absorci´on de la luz solar en la columna vertical de la atm´osfera. Los valores AOD500 nm, AOD380 nm y el contendio en vapor de agua de la columna vertical se han obtenido de la red AERONET [AERONET14, Holben98]. De esta fuente se han podido extraer para la Islas Canarias las medias mensuales de los datos desde a˜no 2008 hasta el a˜no 2014. A partir de este conjunto de datos se calcul´o la media climatol´ogica de la zona en cuesti´on y se mantuvieron los valores obtenidos para la estimaci´on del modelo de cielo despejado para todo el a˜no. La medici´on de AERONET para las Islas Canarias se realiza en la isla de Tenerife. El contenido en Ozono de la atm´osfera se ha obtenido del World Ozone Monitoring Mapping perteneciente al Gobierno Canadiense [Canada’s15]. Al igual que en el caso anterior existe una dato disponible para las Islas Canarias y se han resumido todos los datos obtenidos en una media climatol´ogica anual. 26 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR C0-Pozo Izquierdo C1-Las Palmas Proveedor datos ITC ITC Periodo datos 2003 - 2005 2002 - Jan-Jun 2003 - Jul-Dec 2004 - 2005 Intervalo tiempo toma datos 1 min 1 min Longitud (o) -15.4244 -154269 Coordenadas Latitud (o) 27.8175 28.1108 Geogr´aficas Zona horaria (h) -1 -1 Altitud (m) 47 17 Media presi´on atmosf´erica 101202 101202 (Pa) Columna Ozono (cm) Media 0.3 0.3 Columna vapor agua (cm) Media 1.815 1.815 AOD500nm Media 0.159 0.159 AOD380nm Media 0.184 0.184 Ba (Asymetric factor) 0.84 0.84 Nohoras de cielo despejado 3224 1660 [Ineichen06] rRMSE % 3.97 % 3.98 % Variabilidad del lugar 0.14 0.18 [Hoff12] Tabla 2.1: Variables meteorol´ogicas del modelo de cielo despejado y resultados de su estimaci´on. El modelo de Bird [Bird81] estima la radiaci´on solar global en una superficie horizontal para un d´ıa despejado. A partir de los datos anteriores, se calculan una serie de variables como la dispersi´on de Rayleigh, las absortancias del Ozono, gases como el ox´ıgeno y di´oxido de carbono, el vapor de agua y la absortancia y dispersi´on de los aerosoles. Mediante estas variables se estimar´a finalmente la radiaci´on global incluyendo la radiaci´on exoatmosf´erica y el ´angulo cenital en el momento y lugar especificados. En esta tesis se ha estimado el modelo de cielo despejado de Bird para todas las horas del a˜no en cada una de las estaciones de medida del Instituto Tecnol´ogico de Canarias. En la Tabla 2.1, se pueden observar los datos utilizados y algunos resultados obtenidos al estimar el modelo de cielo despejado para dos estaciones de medida en la isla de Gran Canaria. En este caso se ha optado por mostrar dos estaciones significativas de los dos climas claramente diferenciales en la isla, Pozo Izquierdo-C0 en el sureste de la isla y Las Palmas-C1 en el noreste de la isla. Para comprobar la precisi´on del modelo de cielo despejado de Bird no se pueden comparar los resultados obtenidos con todos los datos disponibles, ya que en muchos casos no corresponden a condiciones de cielo limpio. La comprobaci´on se realizar´a ´unicamente comparando las horas de radiaci´on correspondientes a cielo despejado. En esta tesis se ha utilizado el modelo de Ineichen [Ineichen06] para reconocer los datos correspondientes a condiciones de cielo limpio entre todo el conjunto de datos. Finalmente, se calcular´a el 2.3. DATOS TERRESTRES DE RADIACI ´ ON 27 C0 C1 C2 C4 C5 C6 Nohoras de la muestra 17520 26280 14592 17544 17520 17520 Nohoras de cielo despejado 3224 1660 3159 1675 1448 2765 [Ineichen06] Precisi´on modelo Bird rRMSE % 3.97 3.98 4.06 4.93 4.23 4.01 Variabilidad del lugar [Hoff12] 0.14 0.15 0.13 0.16 0.18 0.20 Tabla 2.2: Estimaci´on de la precisi´on del modelo de cielo despejado de Bird. error cuadr´atico medio relativo ( %rRMSE) resultante de comparar el modelo de cielo despejado estimado con el modelo de Bird y los datos reales correspondientes a cielo despejado. En la Tabla 2.2 se puede observar el error global obtenido al aplicar el modelo de cielo despejado de Bird para cada hora del a˜no en todas las estaciones de medida de la isla de Gran Canaria. El error obtenido para todas las estaciones se encuentra alrededor del 4 %, por lo que el modelo se considera adecuado para nuestro prop´osito. Estos errores %rRMSE se encuentran en mismo orden que el modelo de McClear, que usa los datos de AOD, Ozono y vapor de agua del proyecto MACC [Lefevre13]. Esto siginifica que la variabilidad de estos par´ametros se puede considerar no relevante y por lo tanto utilizar las medias es suficiente. En la Tabla 2.2 no se observa ninguna diferencia apreciable en el error obtenido entre las estaciones del Sur y del Norte. Por otro lado, el n´umero de horas de sol en condiciones de cielo despejado detectadas por el m´etodo de Ineichen es, en general, mayor en las estaciones del Sur. Las estaciones del Norte de la isla presentan mayor n´umero de d´ıas nublados por la presencia de nubes en verano debido a los vientos Alisios. En este apartado tambi´en se propone una variable para estimar la variabilidad de la radiaci´on solar en cada estaci´on. Para cada estaci´on se ha calculado la variabilidad en los datos de radiaci´on solar horarios por el m´etodo propuesto por Hoff & Perez [Hoff12]. La variable propuesta por Hoff & Perez es la desviaci´on t´ıpica de la variaci´on del ´ındice de cielo despejado k∗para dos horas consecutivas. Un lugar que presente un valor superior a 0.2 en esta variables se considera que tiene condiciones clim´aticas inestables. En el caso de Gran Canaria, como se observa, todas las estaciones presentan valores inferiores, aunque siempre menores en las estaciones del Sur cuyo clima es m´as estable. 2.3. Datos de radiaci´on solar disponibles en las estaciones de medida Los modelos de predicci´on que se utilizar´an en esta tesis necesitan un conjunto de datos de radiaci´on solar medidos para poder realizar el aprendizaje. Por lo tanto, es muy importante disponer de una base de datos suficientemente amplia y fiable para que el proceso pueda ofrecer buenos resultados. Antes de comenzar a utilizar los datos se deber´a realizar un tratamiento y an´alisis de los mismos, de manera que la base de datos finalmente obtenida constituya un conjunto lo m´as fiable posible. En esta tesis se llevar´an a cabo los siguientes pasos para conseguir este prop´osito: 28 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR Figura 2.9: Distribuci´on geogr´afica de la estaciones de medida del Instituto Tecnol´ogico de Canarias en la isla de Gran Canaria Estudio del conjunto de datos y del porcentaje de datos sin huecos de los que se dispone en cada estaci´on. Tratamiento de los datos disponibles para filtrar aquellos datos considerados incorrectos. Analizar los datos disponibles. 2.3.1. Estaciones de medida en la isla de Gran Canaria Los datos de radiaci´on solar horaria utilizados en esta tesis provienen de las estaciones de medida del Instituto Tecnol´ogico de Canarias (ITC) en la isla de Gran Canaria. En las Islas Canarias el ITC dispone de 23 estaciones de medida, con las que se ha realizado el mapa solar de la regi´on [D´ıaz12]. En la isla de Gran Canaria existen siete estaciones de medida repartidas por la geograf´ıa insular, ver Figura 2.9. Se dispone de los datos de radiaci´on de dichas estaciones desde a˜no 1998 hasta 2010, aunque, dependiendo de la estaci´on el per´ıodo puede ser menor. Aunque se disponen de siete estaciones, ´unicamente se va a trabajar en esta tesis con seis ya que la estaci´on de G´aldar no dispone de datos suficientes para conseguir un conjunto fiable. Para reconocer las estaciones ´estas se han nombrado con una ’C’ y un 2.3. DATOS TERRESTRES DE RADIACI ´ ON 29 C´odigo estaci´on Latitud Longitud Altitud C0 27.82 -15.42 47 C1 28.11 -15.42 17 C2 27.99 -15.79 197 C4 27.77 -15.58 265 C5 28.03 -15.49 525 C6 27.88 -15.72 300 Tabla 2.3: Coordenadas geogr´aficas de las estaciones de medida del ITC en la isla de Gran Canaria. c´odigo num´erico, de manera que en la isla de Gran Canaria se encuentran las siguientes estaciones: C0-Pozo Izquierdo C1-Las Palmas. C2-La Aldea C4-Maspalomas C5-Sta. Br´ıgida C6-Mog´an C7-G´aldar (no se utilizar´a en esta tesis) Se puede observar que las estaciones cubren la mayor´ıa del territorio insular. Las estaciones C1 y C5 (Las Palmas y Sta. Br´ıgida) en la zona Norte, mientras que las estaciones C0, C2, C4 y C6 (Pozo Izquierdo, La Aldea, Maspalomas y Mog´an) se encuentran en la vertiente sur de la isla. Por otro lado la mayor´ıa de las estaciones se encuentran pr´acticamente en la costa, excepto C6 y C5 que se encuentran en el interior y a mayor altitud sobre el nivel del mar. En la Tabla 2.3 se pueden observar la coordenadas geogr´aficas de latitud, longitud y altitud de cada estaci´on. Estas coordenadas se utilizar´an para realizar cualquier c´alculo o estimaci´on de radiaci´on solar. En las seis estaciones que finalmente se van a utilizar en esta tesis se recogen datos de irradiancia solar en potencia W/m2. Los datos de radiaci´on se recogen cada minuto o cada 5 minutos, dependiendo de la estaci´on y del a˜no en concreto. A partir de los datos almacenados en la base de datos se realizan los c´alculos de la radiaci´on en Wh/m2horaria y diaria integrando los datos de potencia. Aunque en esta tesis, ´unicamente se van a utilizar los datos de radiaci´on global horizontal, en las estaciones de medida se dispone de equipos para medir otras variables meteorol´ogicas. Dependiendo de la estaci´on y del a˜no, nos encontramos con datos de presi´on atmosf´erica, temperatura ambiente, humedad relativa, horas de sol o radiaci´on difusa. 30 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR A˜no ene feb mar abr may jun jul ago sep oct nov dic 2003 100 99.7 96.6 99.7 93.5 100 100 100 100 93.5 100 100 2005 100 96.4 93.5 96.7 100 96.7 96.8 96.8 83.3 100 90 90.3 Tabla 2.4: Porcentaje de datos correctos disponibles en los a˜nos elegidos para la estaci´on de C0-Pozo Izquierdo. Como ya se ha dicho, las estaciones disponen de medidas registradas desde 1998 hasta 2010, excepto C0-Pozo Izquierdo cuyas medidas comenzaron a tomarse a partir de 2001. En cualquier caso, la base de datos muestra lagunas de datos derivadas del mal funcionamiento de los equipos o que se encuentren en revisi´on. Los datos horarios cedidos por el ITC ya hab´ıan sido previamente analizados. El ITC asigna un c´odigo de error negativo a todos los datos recogidos seg´un si se observa un fallo en los instrumentos de medida, fallo en el sensor o en el almacenado de la informaci´on. Cuando se va a proceder a la integraci´on de los datos horarios, ´esta no se realiza para los datos marcados con un c´odigo de error negativo. Los datos horarios cedidos por el ITC exigen un 50 % de datos correctos para calcular las estad´ısticas, en caso contrario se considera que no existe dato de medida en esa hora. Los modelos de predicci´on que se utilizan en esta tesis se basan en el aprendizaje a partir de un conjunto de datos hist´oricos. La predicci´on de datos futuros se basa en la observaci´on de un conjunto de datos pasados. En cualquier proceso de aprendizaje se debe disponer de un conjunto de datos de entrenamiento, con el que se definen los modelos te´oricos, y un conjunto de test, con el que se validan los modelos obtenidos. Ambos conjuntos deben ser comparables y contener una distribuci´on similar de todas las posibles situaciones clim´aticas que se pueden observar en la isla. En esta tesis, se ha decidido trabajar con un a˜no de datos de entrenamiento y un a˜no de datos para la validaci´on, de manera que en ambos casos se encuentren representadas todas las estaciones del a˜no. Los modelos estad´ısticos de series temporales reproducen mejores resultados de predicci´on utilizando bases de datos con una mayor continuidad en los mismos. Es por ello, que una vez analizados los datos se decidi´o trabajar ´unicamente con a˜nos pr´acticamente completos, en la que todos los meses presentaran un porcentaje de datos perdidos los m´as bajo posible. En algunas estaciones no se encontraron a˜nos completos en los que todos los meses presentaran una continuidad suficiente, por lo que se procedi´o a dise˜nar a˜nos ficticios como mezcla de dos a˜nos distintos. Una vez revisados los datos de todas las estaciones se muestran en las tablas (2.4-2.9) los porcentajes de datos de radiaci´on solar IGH de los a˜nos finalmente elegidos para este estudio. Por ´ultimo en la tabla 2.10 se muestran los conjuntos de datos finalmente elegidos para el estudio. Se ha dividido la base de datos en un conjunto de entrenamiento y otro de test para cada estaci´on. 2.3. DATOS TERRESTRES DE RADIACI ´ ON 31 A˜no ene feb mar abr may jun jul ago sep oct nov dic 2002 99.9 100 100 96.5 100 100 100 77.1 80 79.6 100 97.8 2003 99.9 100 99.8 100 99.9 78.75 0 0 0 6.5 96.5 100 2004 100 100 96.6 83.3 79 67.4 98.8 95.2 100 100 100 93 2005 96.6 99.4 77.3 80 98.1 100 100 99.9 100 100 93.1 78.9 Tabla 2.5: Porcentaje de datos correctos disponibles en los a˜nos elegidos para la estaci´on de C1-Las Palmas. A˜no ene feb mar abr may jun jul ago sep oct nov dic 2001 100 99.3 96.8 83.2 93.7 100 100 100 99.9 100 62.8 87 2002 98.8 100 99.9 100 95 100 100 100 99.9 100 89.4 38.4 Tabla 2.6: Porcentaje de datos correctos disponibles en los a˜nos elegidos para la estaci´on de C2-La Aldea. A˜no ene feb mar abr may jun jul ago sep oct nov dic 2000 100 100 100 100 100 100 0 0 0 0 0 0 2003 99.7 89.1 98.9 99.9 100 100 99.1 91.7 96.3 99.1 93.3 91.1 2005 99.9 71.4 99.7 96.8 73.5 38.2 100 100 96.7 100 99.6 98.7 Tabla 2.7: Porcentaje de datos correctos disponibles en los a˜nos elegidos para la estaci´on de C4-Maspalomas. A˜no ene feb mar abr may jun jul ago sep oct nov dic 2001 63.2 42.9 0 0 92.1 100 100 99.3 0 96.8 99.4 98.8 2002 99.6 99.9 83.3 99.4 83.9 6 2 80.4 99.9 99.7 99.4 99.1 2003 99.6 100 99.1 100 100 100 100 99.3 96.7 98.9 98.5 98.9 Tabla 2.8: Porcentaje de datos correctos disponibles en los a˜nos elegidos para la estaci´on de C5-Sta. Br´ıgida. A˜no ene feb mar abr may jun jul ago sep oct nov dic 2002 100 100 100 100 100 100 100 100 100 63.4 100 96.8 2003 100 100 100 100 100 100 100 67.7 96.7 100 100 100 2005 100 100 100 100 35.5 96.7 100 100 96.7 100 99.6 96.1 Tabla 2.9: Porcentaje de datos correctos disponibles en los a˜nos elegidos para la estaci´on de C6-Mog´an. 38 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR (a) C0 (b) C1 (c) C2 (d) C4 (e) C5 (f) C6 Figura 2.15: Distribucion de los datos diarios del ´ındice de cielo despejado para las estaciones de Gran Canaria en funci´on de la media y la variabilidad diaria. 2.3. DATOS TERRESTRES DE RADIACI ´ ON 39 k∗alta y una variabilidad baja std(∆k∗) se representa dentro del conjunto de d´ıas CI. Por otro lado, un d´ıa CIII presentar´ıa una media diaria de radiaci´on solar alta pero con alta variabilidad, mostrando inestabilidad en los niveles de radiaci´on. Los d´ıas tipo A, por contra, ofrecen medias del ´ındice de cielo despejado bajas y alejadas de la curva diaria estimada con el modelo de Bird, por la presencia de nubosidad en la zona. En la Figuras 2.16 y 2.17, se observan, para las estaciones de C0-Pozo Izquierdo y C1Las Palmas, un ejemplo de nueve d´ıas distintos correspondientes a las nueve clasificaciones presentadas. En azul se representa la curva diaria de irradiaci´on estimada con el modelo de cielo despejado de Bird (IGHcs), y con l´ınea roja se observa la curva real de datos de irradiaci´on global horizontal (IGH) medidos en la estaci´on. As´ı, por ejemplo en un d´ıa CI ambas curvas son muy similares debido a en este conjunto se incluyen los d´ıas despejados. Los d´ıas CIII muestran medias altas de radiaci´on, por lo que la curva alcanza los niveles m´aximos de radiaci´on estimados por el modelo de Bird pero aparecen irregularidades en la curva debido a la alta variabilidad. Por otro lado, los d´ıas tipo A, independientemente de si tienen o no alta variabilidad, la curva de radiaci´on siempre se encuentra muy por debajo de la estimaci´on de cielo despejado. Una vez se ha realizado la divisi´on del n´umero de d´ıas en cada estaci´on, es interesante estudiar la distribuci´on en cada conjunto para establecer las condiciones clim´aticas de cada estaci´on y conocer de manera general el comportamiento de la isla en su conjunto. Los resultados obtenidos muestran consistencia con la observaci´on emp´ırica del clima en la isla y confirman que la elecci´on arbitraria de los l´ımites de cada tipo de d´ıa [Dambreville14a] se puede utilizar en Gran Canaria. En las Tablas (2.13-2.18), se muestra la divisi´on de d´ıas seg´un los distintos tipo elegidos en cada estaci´on de medida. En la isla de Gran Canaria existen un mayor n´umero de d´ıas despejados en la parte Sur lo que se traduce en un mayor n´umero de d´ıas tipo C en las estaciones C0 y C2, mientras que en el resto de estaciones prevalecen los d´ıas tipo B. En el norte de la isla, donde se observan una mayor n´umero de d´ıas nublados debido al efecto de los vientos Alisios en su encuentro con las monta˜nas, las estaciones C1 y C5 presentan una mayor´ıa de d´ıas tipo B y III. Como norma general, en la isla se producen un bajo porcentaje de d´ıas tipo I, con baja variabilidad, predominando los d´ıas tipo II y III. De la misma manera, en ninguna estaci´on de medida existe un porcentaje significativo de d´ıas tipo A con baja radiaci´on. Total d´ıas A: alta nubosidad B: nubosidad media C: baja nubosidad 707 6 % 44 % 50 % III: alta variabilidad 34 % 3 26 5 II: variabilidad media 40 % 2 17 21 I: baja variabilidad 26 % 1 1 24 Tabla 2.13: Distribuci´on porcentual del n´umero de d´ıas de cada tipo para la estaci´on de C0-Pozo Izquierdo. Esta distribuci´on de d´ıas en nueve tipos distintos seg´un sus caracter´ısticas clim´aticas se utilizar´a posteriormente para comparar los resultados de las predicciones. Los modelos de 40 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR Figura 2.16: Representaci´on de la curva diaria de radiaci´on IGH (rojo) frente a la curva diaria de cielo despejado (azul) para cada uno de los tipo de d´ıas en la estaci´on de C0-Pozo Izquierdo 2.3. DATOS TERRESTRES DE RADIACI ´ ON 41 Figura 2.17: Representaci´on de la curva diaria de radiaci´on IGH (rojo) frente a la curva diaria de cielo despejado (azul) para cada uno de los tipo de d´ıas en la estaci´on de C1-Las Palmas 42 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR Total d´ıas A: alta nubosidad B: nubosidad media C: baja nubosidad 1045 19 % 69 % 12 % III: alta variabilidad 56 % 9 46 2 II: variabilidad media 39 % 9 23 7 I: baja variabilidad 5 % 1 1 3 Tabla 2.14: Distribuci´on porcentual del n´umero de d´ıas de cada tipo para la estaci´on de C1-Las Palmas. Total d´ıas A: alta nubosidad B: nubosidad media C: baja nubosidad 600 3 % 40 % 57 % III: alta variabilidad 27 % 2 19 6 II: variabilidad media 44 % 1 19 24 I: baja variabilidad 29 % 0 2 27 Tabla 2.15: Distribuci´on porcentual del n´umero de d´ıas de cada tipo para la estaci´on de C2-La Aldea. Total d´ıas A: alta nubosidad B: nubosidad media C: baja nubosidad 721 7 % 77 % 16 % III: alta variabilidad 46 % 4 41 1 II: variabilidad media 44 % 2 36 7 I: baja variabilidad 10 % 1 1 8 Tabla 2.16: Distribuci´on porcentual del n´umero de d´ıas de cada tipo para la estaci´on de C4-Maspalomas. Total d´ıas A: alta nubosidad B: nubosidad media C: baja nubosidad 713 26 % 63 % 11 % III: alta variabilidad 58 % 12 42 4 II: variabilidad media 38 % 13 20 5 I: baja variabilidad 4 % 1 1 2 Tabla 2.17: Distribuci´on porcentual del n´umero de d´ıas de cada tipo para la estaci´on de C5-Sta. Br´ıgida. predicci´on estimar´an la radiaci´on para los distintos tipos de d´ıas presentados y as´ı poder estudiar los m´etodos seg´un las condiciones clim´aticas del a˜no. 2.4. ADQUISICI ´ ON DE DATOS DEL ECMWF 43 Total d´ıas A: alta nubosidad B: nubosidad media C: baja nubosidad 1049 11 % 76 % 13 % III: alta variabilidad 72 % 6 57 9 II: variabilidad media 26 % 4 18 4 I: baja variabilidad 2 % 1 1 0 Tabla 2.18: Distribuci´on porcentual del n´umero de d´ıas de cada tipo para la estaci´on de C6-Mog´an. 2.4. Adquisici´on de datos de un m´etodo num´erico de predici´on (NWP) Los modelos num´ericos de predicci´on ’NWP’ (Numerical Weather Prediction) se emplean para la predicci´on del estado de la atm´osfera desde horas hasta 15 d´ıas de antelaci´on. La predicci´on de los cambios en la atm´osfera, incluyendo la formaci´on y disoluci´on de las nubes, se basan en modelos f´ısicos. Estos modelos f´ısicos se describen mediante ecuaciones diferenciales b´asicas que se resuelven utilizando de m´etodos num´ericos [Lorenz12, Diagne13]. Los NWPs globales calculan el futuro estado de la atm´osfera para una determinada superficie, partiendo de unas condiciones iniciales conocidas mediante las agencias meteorol´ogicas. Los modelos NWP globales se encuentran actualmente operados por 15 diferentes agencias meteorol´ogicas mundiales. Como ejemplo tenemos el Global Forecaste System (GFS) utilizado por la US National Oceanic and Atmospheric Administration (NOAA) y el Integrated Forecast System (IFS) operado por la European Centre for Medium-Range Weather Forecasts(ECMWF). Los modelos globales de predicci´on tienen normalmente una resoluci´on reducida para determinar las variaciones en una distribuci´on espacial muy localizada. Los NWP distribuyen la superficie terrestre en una cuadr´ıcula cuya precisi´on var´ıa. Las predicciones de los NWP cubren la superficie con una resoluci´on aproximada desde 0.125 grados hasta 0.5 grados actualmente dependiendo del modelo. La resoluci´on temporal var´ıa de 1 hora para los modelos t´ıpicos regionales y de 3 a 6 horas para los modelos globales. Los c´alculos realizados por los modelos NWP globales comienzan bas´andose en unos datos conocidos del estado de la atm´osfera. Los modelos globales obtienen esta informaci´on de la red de medidas meteorol´ogicas realizadas por las distintas agencias internacionales. Las variables clave con las que llevan a cabo estos c´alculos son la temperatura, humedad, la presi´on atmosf´erica, las condiciones espaciales del viento en las tres dimensiones o la temperatura de la superficie oce´anica entre otras. Los modelos NWP globales cubren la superficie completa de la Tierra. Muchos procesos f´ısicos ocurren para escalas espaciales mucho menores que la precisi´on de los modelos globales, como podr´ıan ser la dispersi´on y absorci´on de la radiaci´on o los procesos de condensaci´on y convecci´on que influyen en la formaci´on de nubes. Para mejorar la predicci´on local, tanto espacial como temporalmente, se han desarrollado modelos locales y regionales que reducen la distribuci´on de la cuadr´ıcula obteniendo superficies de 3 a 10 44 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR km2de precisi´on para datos horarios, aunque se pueden conseguir resoluciones superiores. Estos modelos s´olo cubren una parte de la superficie terrestre y utilizan como datos de partida los datos obtenidos por los modelos NWP globales como condiciones de contorno. Estos modelos son normalmente desarrollados por las distintas agencias nacionales de meteorolog´ıa o por compa˜n´ıas privadas. La precisi´on de estos modelos se est´a mejorando actualmente utilizando modelos estad´ısticos de predicci´on. Estos modelos, que se basan en el aprendizaje del comportamiento de los datos hist´oricos registrados para una localidad, est´an consiguiendo buenos resultados para diferentes horizontes de predicci´on combinando los datos hist´oricos con datos provenientes de modelos num´ericos de predicci´on. En general, se est´an desarrollando mejoras en los modelos de predicci´on estad´ısticos combinando datos hist´oricos de medida, modelos de predicci´on num´ericos o datos provenientes de la predicci´on sat´elite. 2.4.1. Datos de radiaci´on solar del ECMWF Como se ha explicado, existen diferentes modelos num´ericos de escala global desarrollados por diferentes agencias meteorol´ogicas. En esta tesis se han empleado valores de predicci´on del Integrated Forecast System (IFS) operado por la European Centre for Medium-Range Weather Forecasts(ECMWF). Este modelo proporciona las predicciones, incluyendo la irradiancia en la superficie terrestre y diferentes par´ametros de las nubes, con una precisi´on espacial horizontal de alrededor de 16 km x 16 km y 137 niveles verticales de resoluci´on. Los datos de salida est´an disponibles con una resoluci´on temporal de 3 horas, quedando disponible una resoluci´on horaria para trabajos de investigaci´on. Las predicciones del ECMWF han demostrado su calidad tanto para datos de radiaci´on solar como de condiciones de viento. En esta tesis se han utilizado datos de predicciones realizadas por el ECMWF para la isla de Gran Canaria gracias al Laboratoire de Physique et Ing´enierie Math´ematique pour l’Energie et l’environnement (PIMENT) de la Universidad de la isla de La Reuni´on. En concreto se obtuvieron datos para una superficie que cubre la isla con un rango de 27.5o a 28.5ode latitud Norte y de 15oa 16ode longitud Oeste y una resoluci´on de 0.125 o tanto en longitud como en latitud. En la Figura 2.18 se puede observar la distribuci´on de la cuadr´ıcula sobre la isla y la posici´on en rojo de cada una de las estaciones del ITC. Para cada elemento de la cuadr´ıcula se disponen de todos los datos de la predicci´on solicitada para los a˜nos entre el 2000 y 2005. Durante estos a˜nos el ECMWF ´unicamente puede ofrecer una resoluci´on temporal de 3 horas. As´ı, se disponen de las predicciones realizadas por el modelo num´erico de predicci´on ECMWF para la distribuci´on espacial especificada cada tres horas para el d´ıa siguiente (horizonte de una d´ıa de adelanto). La gesti´on de la redes el´ectricas necesita habitualmente de datos de radiaci´on horaria y adem´as, en esta tesis en concreto, se ha decidido trabajar con datos de radiaci´on horarios. En la bibliograf´ıa se encuentran diversas t´ecnicas de interpolaci´on de los datos obtenidos por el ECMWF cada 3 horas para convertirlos en datos horarios. En este trabajo se ha utilizado la t´ecnica descrita por el Proyecto Endorse [Espinar11] en el que se tiene en cuenta la conservaci´on de la energ´ıa solar en la integraci´on. 2.4. ADQUISICI ´ ON DE DATOS DEL ECMWF 45 Figura 2.18: Distribuci´on geogr´afica de los datos obtenidos del modelo de predicci´on num´erica ECMWF para la isla de Gran Canaria 46 CAP´ ITULO 2. FUNDAMENTOS Y DATOS DE RADIACI ´ ON SOLAR Una vez descargados los datos predecidos por el ECMWF para cada uno de los elementos de la cuadr´ıcula en los que se encuentren las estaciones de medida del ITC se interpolar´an. De esta manera, se dispondr´a de una base de datos horaria para los mismos a˜nos en los que se dispon´ıa de datos terrestres en las estaciones de medida. Las variables que se descargaron para dichas fechas y con las que se trabajar´a en esta tesis son las siguientes: Latitude, ofrece la latitud central de cada elemento de la cuadr´ıcula. Longitude, ofrece la longitud central de cada elemento de la cuadr´ıcula. Time, marca el momento al que corresponde cada medida. TCC Total Cloud Cover, se obtiene el ´ındice de nubosidad entre 0 (sin nubes) y 1 (totalmente nuboso) derivado del estudio de todos los niveles de altitud. SSRD Surface Solar Radiation Downwards, en este dato se dispone de la radiaci´on solar horizontal acumulada entre dos instantes sucesivos en J/m2. 2.4.2. Datos de viento y humedad del ECMWF Adem´as de los datos de radiaci´on, en esta tesis se estudiar´an los efectos de incluir los datos de velocidad y direcci´on del viento y la humedad relativa de la atm´osfera. En las estaciones de medida no se dispone de estos datos por lo que se opt´o por incluir las predicciones que realiza el ECMWF de estas variables. As´ı, de la misma manera que en el caso anterior se obtuvieron las siguientes medidas para la superficie terrestre con id´entica resoluci´on espacial y temporal. En este caso los datos extra´ıdos del programa ´unicamente est´an disponibles para el a˜no 2005, por lo que su contribuci´on a la predicci´on de la radiaci´on solar se estudiar´a en las estaciones de C0-Pozo Izquierdo y C1-Las Palmas, como estaciones representativas del clima en el sur y norte de la isla respectivamente. Latitude, ofrece la latitud central de cada elemento de la cuadr´ıcula. Longitude, ofrece la longitud central de cada elemento de la cuadr´ıcula. Time, marca el momento al que corresponde cada medida. LCC Low Cloud Cover, se obtiene el ´ındice de nubosidad entre 0 (sin nubes) y 1 (totalmente nuboso) derivado del estudio de los niveles entre la superficies y el 0.8 de la presi´on en la superficie. MCC Medium Cloud Cover, se obtiene el ´ındice de nubosidad entre 0 (sin nubes) y 1 (totalmente nuboso) derivado del estudio de los niveles entre el 0.8 y el 0.45 de la presi´on en la superficie. HCC High Cloud Cover, se obtiene el ´ındice de nubosidad entre 0 (sin nubes) y 1 (totalmente nuboso) derivado del estudio de los niveles entre el 0.8 de la presi´on en la superficie y el nivel m´as alto estudiado por el modelo. 2.4. ADQUISICI ´ ON DE DATOS DEL ECMWF 47 En este caso, no interesaba ´unicamente los datos de viento y humedad en la superficie terrestre. Se busca incluir los datos de viento y humedad a distintas altitudes de la atm´osfera por lo que se extrajeron los datos del ECMWF para todos los niveles disponibles, obteni´endose diferentes medidas entre las altitudes de 1000-1 hPa de presi´on. R Relative Humidity, se define con respecto a la saturaci´on sobre el hielo por debajo de 23oC y con respecto a la saturaci´on sobre el agua. U-velocity, velocidad del viento para la coordenada U del eje de coordenadas. V-velocity, velocidad del viento para la coordenada V del eje de coordenadas. Aunque se realiz´o la extracci´on de los datos para todas las altitudes se utilizar´an los correspondientes a ciertas altitudes en los modelos estad´ısticos elegidos. En concreto en esta tesis se decidi´o incluir los datos de viento y humedad para las siguientes altitudes de la atm´osfera: Altura de la base de las nubes, como en el ECMWF para el a˜no 2005 no se dispon´ıa de este dato, se obtuvo esta altitud como aquella a la cual la humedad relativa es mayor o igual al 95 % [Lave13]. En este caso por lo tanto ´unicamente se dispone en cada instante de los datos de m´odulo y direcci´on del viento. Altura de 700 HPa, se obtuvieron los datos de velocidad y direcci´on del viento y de humedad relativa para esta altitud siguiendo un estudio en la Isla Reuni´on donde se indicaba esta altitud como la m´as influyente en la predicci´on de la radiaci´on solar [Badosa15]. Esta altitud se define como aquella a la cual la humedad relativa presenta una mayor desviaci´on t´ıpica. Altura de la inversi´on t´ermica en las Islas Canarias seg´un el estudio [DA96], se obtuvieron los datos de velocidad y direcci´on del viento y de humedad relativa para la altitud en que se produce la inversi´on t´ermica en cada mes del a˜no. El ECMWF ofrece los datos del viento seg´un sus coordenadas (u,v) pero en los modelos de predicci´on estos valores se utilizaron como el m´odulo y direcci´on del vector correspondiente. 54 CAP´ ITULO 3. RADIACI ´ ON SOLAR SATELITAL meteorol´ogicos conocidos como la presencia est´atica de nubosidad sobre la ciudad de Las Palmas, en el noreste de la isla, durante el verano. Este fen´omeno se representa en dos diferentes escenarios, en uno se observa la acumulaci´on de nubes aparentemente debidas al movimiento de los vientos Alisios predominantes, ver Figura 3.4, y en otro se observa la formaci´on local de las nubes sobre la ciudad sin que aparentemente ´estas hayan sido provocadas debidas al movimiento, ya que no se observa su presencia durante las primeras horas, ver Figura 3.5. En las Figuras 3.4 y 3.5, la secuencia superior muestra la evoluci´on durante un d´ıa cada 30 minutos del ´ındice de cielo despejado para toda la superficie de datos satelitales disponible. Mientras que la secuencia inferior muestra la secuencia de datos de radiaci´on solar para la misma superficie. La presencia de nubes se puede inferir estudiando la secuencia de datos del ´ındice de cielo despejado. En la Figura 3.4 se puede observar claramente la gran influencia de los vientos Alisios predominantes en la formaci´on de nubes sobre todo el frente norte de la isla. Esta nubosidad permanece estable durante la ma˜nana, incluso con la presencia de fuertes vientos, debido al fen´omeno de inversi´on t´ermica en conjunto con la orograf´ıa escarpada de la zona norte [D´ıaz12]. Se observa que las nubes siguen un movimiento de norte a sur principalmente debido a los vientos predominantes en la zona donde no se ven interrumpidas por la orograf´ıa de la isla. Por otro lado, en la Figura 3.5, se muestra la creaci´on de nubosidad est´atica de nubes local sin ninguna aparente contribuci´on de nubes provenientes de la direcci´on principal de los vientos Alisios. Esta formaci´on de nubes aparecen y crecen gracias a la acumulaci´on de part´ıculas de agua presentes en las corrientes de aire y que permanecen indetectables para las im´agenes satelitales. Esta observaci´on nos llevar a pensar que los m´etodos basados ´unicamente en el movimiento de las nubes no van a proporcionar toda la informaci´on necesaria para predecir la formaci´on de nubes en las zonas monta˜nosas. Aunque estos dos sucesos aparentemente muestran diferentes condiciones clim´aticas, ambos son provocados por el mismo fen´omeno principal, la inversi´on t´ermica en el norte de la isla y los vientos Alisios. Estos ejemplos son un simple vistazo a la complejidad en la predicci´on meteorol´ogica cuando la disipaci´on, creaci´on y evoluci´on de las masas de vapor de agua est´an sujetas a gran cantidad de variables. 3.3. An´alisis espacio temporal de los datos Los datos satelitales se ha estudiado para utilizarlos, junto con los datos terrestres, a la hora de mejorar las predicciones de radiaci´on solar. Para estudiar la relaci´on con el clima y establecer alg´un tipo de relaci´on ´util entre ambas bases de datos, se opt´o por utilizar la correlaci´on de Pearson. Con este estudio se pretende establecer qu´e elementos de la cuadr´ıcula de datos sat´elites tienen mayor relaci´on con la radiaci´on medida en las dos estaciones terrestres. Hay que recordar que se va a trabajar con las estaciones de C0-Pozo Izquierdo y C1-Las Palmas, por lo que siempre se estudiar´a por separado la relaci´on de todos los p´ıxeles satelitales con ambas estaciones. 3.3. AN ´ ALISIS ESPACIO TEMPORAL DE LOS DATOS 55 Figura 3.4: Evoluci´on intradiaria de la radiaci´on global horizontal y del ´ındice de cielo despejado de los datos satelitales para cada 30 minutos. D´ıa 07/07/2005. Figura 3.5: Evoluci´on intradiaria de la radiaci´on global horizontal y del ´ındice de cielo despejado de los datos satelitales para cada 30 minutos. D´ıa 12/08/2005. 56 CAP´ ITULO 3. RADIACI ´ ON SOLAR SATELITAL 3.3.1. Correlaci´on con el ´ındice de cielo despejado Primero se realiz´o la correlaci´on entre las dos bases de datos (datos satelitales y datos terrestres de cada estaci´on) utilizando los ´ındices de cielo despejado para estudiar la relaci´on entre ambos [Dambreville14b]. Para poder evaluar la correlaci´on entre las series temporales en diferentes momentos se establecieron desfases temporales entre ambas series. Los desfases temporales seleccionados se mueven desde h= 0, ambas series sin ninguna diferencia temporal, y h= 3, la serie de datos satelitales retrasada 3 horas respecto a los datos terrestres, ecuaci´on (3.1). De esta manera se puede establecer una relaci´on entre los datos sat´elites en los instantes pasados y el dato medido en la superficie terrestre de la estaci´on objeto de estudio. En cada estaci´on se dispondr´a de las correlaciones entre cada p´ıxel de los datos sat´elite desfasado hasta 3 horas y el dato terrestre correspondiente a la propia estaci´on sin desfasar. Ck∗(i, j)h=corr(k∗ ground(t), k∗ satellite(t−h)) para h = 0,1,2 & 3 (3.1) En las Figuras 3.6 y 3.7 se muestran los resultados de las correlaciones anuales para cada estaci´on. Luego con estas im´agenes se pueden deducir los p´ıxeles con una mayor relaci´on con nuestra estaci´on, y que por tanto aportar´an una informaci´on m´as relevante para mejorar la predicci´on. En ambas estaciones se han obtenido resultados que se parecen al comportamiento esperado. En la estaci´on C1-Las Palmas, Figura 3.7, las correlaciones m´as altas se obtienen con los p´ıxeles de la zona norte de la isla, mientras que en la estaci´on C0-Pozo Izquierdo la relaci´on es con la zona sur, Figura 3.6. En ambos casos, los valores obtenidos en las correlaciones decaen conforme aumenta el desfase temporal entre los datos sat´elites y terrestres. Un vistazo m´as exhaustivo muestra correlaciones inversamente proporcionales entre el norte y el sur, dividiendo la isla en dos sectores. Este factor es crucial para permitir una mejor compresi´on del microclima de la zona. Otro fen´omeno importante que queda reflejado en las im´agenes es que el norte de la isla est´a influenciado por los vientos Alisios, mientras que las monta˜nas protegen el sur del efecto de los mismos. Todos estas observaciones encajan perfectamente con el conocimiento que se tiene de que el sur de la isla tiene un mayor n´umero de d´ıas despejados que el norte. En las correlaciones anuales la correlaci´on se realiza incluyendo toda la base de datos satelitales y terrestres, con lo que quedan mezclados en los resultados todas las relaciones existentes en las diferentes estaciones del a˜no. As´ı, para poder observar las posibles diferencias en las relaciones entre los datos sat´elites en las distintas estaciones del a˜no, se calcularon correlaciones trimestrales independientes [Zagouras15]. Estas correlaciones trimestrales est´an divididas en los trimestres meteorol´ogicos para dividir los datos seg´un su parecido clim´atico. Los resultados obtenidos y representados en las im´agenes de las correlaciones trimestrales para todos los desfases temporales son coherentes con las condiciones clim´aticas observadas. La estaci´on m´as representativa es el verano por la presencia de los vientos Alisios, como se observa en las Figuras 3.8 y 3.9. Sobre todo en las costa este, se puede observar claramente como los vientos Alisios provocan anomal´ıas sobre la superficie ocupada por el oc´eano, mientras que la zona sur de la isla se mantiene protegida 3.3. AN ´ ALISIS ESPACIO TEMPORAL DE LOS DATOS 57 Figura 3.6: Mapa de intercorrelaci´on anual para el ´ındice de cielo despejado entre los datos de tierra y cada p´ıxel de datos satelitales con un retraso de h = 0, 1, 2 & 3 horas en la estaci´on C0-Pozo Izquierdo. 58 CAP´ ITULO 3. RADIACI ´ ON SOLAR SATELITAL Figura 3.7: Mapa de intercorrelaci´on anual para el ´ındice de cielo despejado entre los datos de tierra y cada p´ıxel de datos satelitales con un retraso de h = 0, 1, 2 & 3 horas en la estaci´on C1-Las Palmas. 3.3. AN ´ ALISIS ESPACIO TEMPORAL DE LOS DATOS 59 Figura 3.8: Mapa de intercorrelaci´on en verano para el ´ındice de cielo despejado entre los datos de tierra y cada p´ıxel de datos satelitales con un retraso de h = 0, 1, 2 & 3 horas en la estaci´on C0-Pozo Izquierdo. 60 CAP´ ITULO 3. RADIACI ´ ON SOLAR SATELITAL Figura 3.9: Mapa de intercorrelaci´on en verano para el ´ındice de cielo despejado entre los datos de tierra y cada p´ıxel de datos satelitales con un retraso de h = 0, 1, 2 & 3 horas en la estaci´on C1-Las Palmas. 3.3. AN ´ ALISIS ESPACIO TEMPORAL DE LOS DATOS 61 Figura 3.10: Mapa de intercorrelaci´on en oto˜no para el ´ındice de cielo despejado entre los datos de tierra y cada p´ıxel de datos satelitales con un retraso de h = 0, 1, 2 & 3 horas en la estaci´on C0-Pozo Izquierdo. 62 CAP´ ITULO 3. RADIACI ´ ON SOLAR SATELITAL Figura 3.11: Mapa de intercorrelaci´on en oto˜no para el ´ındice de cielo despejado entre los datos de tierra y cada p´ıxel de datos satelitales con un retraso de h = 0, 1, 2 & 3 horas en la estaci´on C1-Las Palmas. 3.3. AN ´ ALISIS ESPACIO TEMPORAL DE LOS DATOS 63 por las monta˜nas. En la Figura 3.8 incluso se puede observar como existe un ´area de la imagen en la derecha con una relaci´on mayor provocada por el mismo efecto en la isla de Fuerteventura. Por ´ultimo, se observa claramente, al igual que anualmente, la relaci´on inversamente proporcional entre la estaci´on del sur con los p´ıxeles del norte de la isla y viceversa. En contraste con las im´agenes del verano, se observa en las Figuras 3.10 y 3.11, que representan el oto˜no, como la superficie de p´ıxeles con mayor correlaci´on se concentra alrededor de las estaciones. Se conserva la divisi´on norte-sur de la isla pero no muestra una influencia de los vientos predominantes en los alrededores. Todos estos an´alisis conducen a la conclusi´on de que ambas correlaciones estudiadas, tanto la cuatrimestral como la anual, presentan informaci´on adicional al estudio de la radiaci´on. Adem´as del estudio espacial, nos permite tener informaci´on de la relaci´on en diferentes horizontes temporales. Al demostrar que existe una relaci´on fuerte entre los datos sat´elites con distintos desfases temporales y el dato terrestre, nos permite pensar que la inclusi´on de algunos p´ıxeles satelitales mejorar´an la predicci´on. Este estudio ser´a utilizado en el Cap´ıtulo 5 para realizar la elecci´on de los p´ıxeles m´as relevantes. A diferencia del estudio de Dambreville [Dambreville14b], que selecciona ´unicamente 9 p´ıxeles independientemente del horizonte temporal, se ha optado por no limitar el n´umero de p´ıxeles y variar el n´umero en cada desfase temporal dependiendo de las correlaciones obtenidas. 3.3.2. Correlaci´on con la variaci´on del ´ındice de cielo despejado De la misma manera que se procedi´o con las correlaciones entre los ´ındices de cielo despejado de las bases de datos terrestres y satelitales, se intent´o conseguir informaci´on adicional con una correlaci´on similar utilizando la variaci´on horaria entre los ´ındices de cielo despejado [Dambreville14b]. Esta informaci´on deber´ıa aportar la direcci´on de los sucesos climatol´ogicos que est´an llegando a la isla. Las f´ormulas que describen el proceso seguido se observan en las Ecuaciones (3.2) y (3.3). Adem´as, como en el caso descrito en el apartado 3.3.1, tambi´en se realiz´o un an´alisis anual y trimestral de las bases de datos disponibles para el a˜no 2005. ∆k∗=k∗(t+ 1) −k∗(t) (3.2) C∆k∗(i, j)h=corr(∆k∗ ground(t),∆k∗ satellite(t−h)) para h = 0,1,2 & 3 (3.3) En ambos casos, anual y trimestral, esta correlaci´on basada en la variaci´on del´ındice de cielo despejado se descart´o porque los resultados obtenidos aportaban informaci´on err´atica. Adem´as, esta informaci´on no se consider´o relevante debido a los bajos valores obtenidos en las correlaciones. Como se observa en la Figura 3.12, para cada desfase temporal, los m´aximos valores de correlaci´on obtenidos son menores de 0.1. Un an´alisis m´as profundo de la informaci´on, realizando una aproximaci´on a una distribuci´on normal de los valores para cada desfase temporal, Figura 3.13, se concluye que las curvas obtenidas no ofrecen informaci´on relevante ya que todos los valores se concentran muy cerca del cero, excepto 70 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR 4.3. Modelos lineales La radiaci´on solar est´a considerada un proceso estoc´astico, cuyos niveles esperados dependen de los valores anteriores de la serie temporal y de la variaciones meteorol´ogicas que se puedan producir en la zona de estudio, en un rango temporal que puede ir desde los segundos hasta los d´ıas. Las series temporales que describen fen´omenos f´ısicos se pueden definir como estacionarias o no estacionarias. Las series estacionarias se mantienen constantes respecto a su tendencia general y las fluctuaciones pueden aparecer de manera completamente aleatoria. Como ya se ha explicado, las series temporales de radiaci´on solar se transforman en una serie estacionaria dividiendo por el valor de radiaci´on estimado por un modelo de cielo despejado [Diagne13]. Los modelos estad´ısticos lineales han sido utilizados en el modelado de series temporales durante a˜nos en diversos campos, como econom´ıa, demograf´ıa, meteorolog´ıa o medio ambiente. Entre los modelos lineales m´as utilizados nos encontramos: Modelos Autorregresivos (AR): Se trata de un modelo de regresi´on lineal basado ´unicamente en datos pasados de la serie temporal que se estudia. El orden de modelo viene definido por la variable p, que representa el n´umero de par´ametros de los que consta la funci´on de transferencia. Modelos de medias m´oviles (MA): En este caso el orden de modelo viene definido por la variable q, que representa el n´umero de par´ametros de los que consta la funci´on de transferencia. Modelo Autorregresivos de medias m´oviles (ARMA): Se trata de un modelo formado por dos partes, una basada en el modelo Autorregresivo (AR) y otra basada en el modelo de medias m´oviles (MA). La complejidad del modelo vendr´a definida por los par´ametros pyq, que representan el orden de la parte (AR) y (MA) respectivamente. Modelo Autorregresivos Integrados de media m´oviles (ARIMA): Se trata de un modelo estad´ıstico utilizado para encontrar patrones para una predicci´on futura a partir de variaciones y regresiones de datos estad´ısticos. Las complejidad del modelo vendr´a definida por los par´ametros p,dyq, que representan el orden de la parte (AR), la parte de integraci´on (I) y (MA) respectivamente. Modelo lineales con datos ex´ogenos: En este caso los modelos utilizan, adem´as de los datos pasados de la serie temporal, otros provenientes de series temporales ex´ogenas a la que se pretende modelar. Se encuentran t´ecnicas autorregresivas (ARX), Autorregresivas de medias m´oviles (ARMAX) o Autorregresivas Integradas de medias m´oviles (ARIMAX). Estas t´ecnicas has sido ampliamente utilizadas desde los a˜no setenta en el campo de la radiaci´on solar. Los modelos de series lineales se han utilizado para generar series sint´eticas de valores o para definir A˜nos Meteorol´ogicos T´ıpicos [Mazorra10] y para realizar predicciones de radiaci´on solar. Aguiar & Collares-Pereira [CP89] realizaron trabajos con modelos ARMA a partir de datos del ´ındice de claridad horario y J. Boland 4.3. MODELOS LINEALES 71 [Boland95, Boland08] realiz´o una descripci´on de la radiaci´on solar diaria y horaria, un estudio del modelado de series temporales de radiaci´on solar y un m´etodo de estimaci´on de la radiaci´on solar difusa en Australia. En esta tesis, como se podr´a ver en el cap´ıtulo 5, se han evaluado dos modelos lineales basados ´unicamente en valores pasados, el modelo (AR) y el modelo (ARMA). 4.3.1. Modelo lineal Autorregresivo AR En un modelo Autorregresivo (AR) [Chatfield13], la predicci´on de un valor de la serie temporal en un horizonte temporal h, se asume vendr´a representada por una combinaci´on lineal de valores pasados de la propia serie, seg´un la ecuaci´on (4.6), b k∗(t+h) = p−1 X i=0 [Φi+1k∗(t−i)] + t+h(4.6) donde b k∗(t+h) representa el valor de radiaci´on en el horizonte temporal hytes un ruido blanco con varianza σ2. Los valores k∗(t−i) son los datos pasados de la serie temporal que se han elegido para establecer una relaci´on lineal con el valor a predecir, mientras que los par´ametros del modelo son p, que representa el orden de AR o, lo que es lo mismo, el n´umero de valores pasados utilizados, y {Φi}i=1,2,...,p, que muestra los par´ametros de autorregresi´on obtenidos a partir de los datos de la muestra durante el entrenamiento. Para optimizar el modelo AR, una de las claves del desarrollo del mismo es determinar el n´umero de orden del mismo, p. En esta tesis, se proponen m´etodos basados en el estudio de la muestra de Funciones de Autocorrelaci´on Parcial [Boland95, Boland08], PACF por sus siglas en ingl´es (Partial Autocorrelation Function), y el Criterio de Informaci´on Bayesiana, BIC por sus siglas en ingl´es (Bayesian Information Criterion). 4.3.2. Modelo lineal autorregresivo de medias m´oviles ARMA En los modelos Autorregresivos de medias m´oviles ARMA, los valores futuros de la serie temporal est´an basados en dos modelos b´asicos, un modelo Autorregresivo (AR) y un modelo de medias m´oviles (MA), como una combinaci´on lineal de una cantidad determinada de valores pasados de la serie y de errores, seg´un la ecuaci´on (4.7) b k∗(t+h) = p−1 X i=0 [Φi+1k∗(t−i)] + t+h+ q−1 X j=0 [Θj+1t−j] (4.7) donde b k∗(t+h) representa el valor del ´ındice de cielo despejado en el horizonte temporal h. Los valores k∗(t−i) son los datos pasados de la serie temporal que se han elegido para establecer una relaci´on lineal con el valor a predecir y tes una serie de ruido blanco con con media cero y varianza σ2. En este caso, los par´ametros del modelo son p, que representa el orden de AR o, lo que es lo mismo, el n´umero de valores pasados utilizados, {Φi}i=1,2,...,p, que muestra los par´ametros de autorregresi´on obtenidos a partir de los datos de la muestra durante el entrenamiento, q, que indica orden del modelo MA, 72 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR y{Θi}i=1,2,...,q que representan los par´ametros que acompa˜nan a la serie de errores t. Los par´ametros de los modelos, tanto AR como MA, se obtienen mediante una regresi´on por m´ınimos cuadrados resultante de la comparaci´on del conjunto de datos pasados utilizados como entrada, y de datos futuros que se desean obtener [Box98]. Al igual que para el modelo AR, existen diversas t´ecnicas para establecer los ´ordenes del modelo AR, p, y del modelo MA, q, con lo que el modelo quedar´a definido como ARMA(p,q). En esta tesis, se proponen m´etodos basados en el estudio de la muestra de Funciones de Autocorrelaci´on Parcial, PACF por sus siglas en ingl´es (Partial Autocorrelation Function), la muestra de Funciones de Autocorrelaci´on Simple, ACF por su siglas en ingl´es (Autocorrelation Function), y el Criterio de Informaci´on Bayesiana, BIC por sus siglas en ingl´es (Bayesian Information Criterion)[Boland95, Boland08]. La popularidad de las t´ecnicas ARMA radica es su flexibilidad para representar diferentes tipos de series temporales seg´un el orden de los modelos. Las series temporales deben ser series estacionarias para poder realizar ajustes adecuados con estos modelos [Hamilton94], habiendo demostrado ser adecuadas para la predicci´on de valores futuros. En la mayor´ıa de los casos, se ha demostrado que para series estacionarias los modelos ´optimos obtenidos tienen ´ordenes pyqno superiores a dos [Box98]. 4.3.3. Estudio de la complejidad de los modelos lineales En los modelos lineales utilizados en esta tesis una de las decisiones m´as importantes ser´a el orden de los modelos AR y MA respectivamente. Para ello, siguiendo el m´etodo descrito por J. Boland [Boland95, Boland08], utilizaremos las Funciones de Autocorrelaci´on y Autocorrelaci´on Parcial, ACF y PACF, y los criterios de decisi´on bayesiana, BIC. Para entender el c´alculo de estos par´ametros, primero realizaremos una serie de definiciones. Partiendo de dos series temporales XeYcon medias µXyµY respectivamente, se define la covarianza de XeYseg´un la ecuaci´on (4.8) Cov(X, Y ) = E(X−µX)(Y−µY) (4.8) donde E(x) representa el valor esperado de una variable. Si las series XeYno son independientes la covarianza toma valores positivos o negativos. En el caso de que para valores altos de Xtiendan a coincidir valores altos de Y, la covarianza ser´a positiva, mientras que si Ymuestra valores bajos en el mismo instante, la covarianza ser´a negativa. El coeficiente de correlaci´on se obtiene dividiendo la covarianza por el producto de la desviaciones est´andar de ambas series (4.9) ρ=Corr(X, Y ) = E[(X−µX)(Y−µY)] pE(X−µX)2(Y−µY)2(4.9) Si en lugar de dos series temporales aleatorias, se trabaja con dos series temporales que representan el mismo proceso estoc´astico en diferentes momentos, XtyXt+τ, la expresi´on anterior se conoce como el coeficiente de autocorrelaci´on para el instante τ(4.10) 4.3. MODELOS LINEALES 73 ρτ=Corr(Xt, Xt+τ) = E[(Xt−µt)(Xt+τ−µt+τ)] pE(Xt−µt)2(Xt+τ−µt+τ)2(4.10) En una serie temporal se podr´ıa calcular, variando el valor de τ, un coeficiente de autocorrelaci´on entre dos instantes temporales cualesquiera. Por otro lado, en el caso de que la serie temporal Xtse considere un proceso estacionario, las desviaciones t´ıpicas de todas las variables en el tiempo son id´enticas y por lo tanto la expresi´on de la Funci´on de Autocorrelaci´on (ACF) es seg´un (4.11) ρτ=Corr(Xt, Xt+τ) = E[(Xt−µt)(Xt+τ−µt+τ)] E(Xt−µ)2(4.11) La autocorrelaci´on se puede calcular para todos los valores de la serie temporal Xt=X1, X2, ...XN, y as´ı obtener la autocorrelaci´on muestral SACF para cada diferencia temporal τ(4.12) ˆρτ= N−τ P t=1 [(xt−µ)(xt+τ−µ)] N P t=1 (xt−µ)2 (4.12) donde µrepresenta la media de toda la muestra de la serie temporal Xt, compuesta por Nvalores, y ˆρτser´a el valor de cada una de las Autocorrelaciones de la muestra. De esta manera, la Autocorrelaci´on muestral SACF es una medida de la relaci´on lineal que existe entre dos instantes cualesquiera dentro de la muestra, separados por un tiempo τ. La autocorrelaci´on puede tomar valores entre -1 y +1, siendo la relaci´on lineal m´as fuerte cuanto m´as se acerque a ambos extremos. En el caso de valores positivos, la relaci´on entre los distintos valores de la serie temporal ser´a directamente proporcional. Adem´as, si el valor de correlaci´on es alto para una diferencia de tiempo τ= 1, se est´a indicando una fuerte relaci´on entre XtyXt−1,Xt−1yXt−2, y as´ı hasta N. Los valores de autocorrelaci´on nos indican la relaci´on lineal existente entre los instantes temporales de la serie. Por lo tanto, los valores significativos de la muestra de ACF representan el n´umero de valores pasados que pueden ser necesarios para generar un modelo lineal. En el caso que nos ocupa, considerando un proceso estacionario gaussiano y que para valores altos de τla autocorrelaci´on tiende a cero, se puede considerar que el 95 % se encuentran dentro del intervalo [−1.96std(ρτ),+1.96std(ρτ)]. La autocorrelaci´on parcial PACF establece la correlaci´on entre dos instantes de tiempo de la serie con un retardo τ, sin tener en cuenta la dependencia creada por los retardos intermedios existentes entre ambos instantes, es decir, muestra la correlaci´on pura. La muestra de Autocorrelaci´on Parcial SPACF para todos los retardos de tiempo τse puede obtener a partir de la SACF por medio de las ecuaciones de Yule-Walker. Los par´ametros ˆ φττ ofrecen la estimaci´on de SPACF para cada retardo a partir la ecuaci´on (4.13) 74 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR Figura 4.1: Ejemplo de SACF y SPACF para serie temporal estacionaria aleatoria. En rojo se representan los intervalos de elecci´on ˆ φττ = ρτ+ τ−1 P t=1 ˆ φτ−1,t ˆρτ−1 1 + τ−1 P t=1 ˆ φτ−1,t ˆρt (4.13) De la misma manera que en el caso de SACF, se busca identificar los valores de SPACF significativos para establecer el retardo a partir del cual la relaci´on entre las variables es pr´acticamente nula, ver Figura 4.1. Se puede considerar que el 95 % de los valores se encuentran dentro del intervalo (−1.96/√N, +1.96/√N). Una vez calculados las muestras de SACF y SPACF para la serie temporal estacionaria, se puede utilizar esta informaci´on para estimar el orden de los modelo lineales AR(p), MA(q) y ARMA(p,q). El criterio general, seg´un J. Boland [Boland95, Boland08] es: Cuando el SACF decae paulatinamente y el SPACF contiene picos significativos hasta un valor de retardo τ=p, se considera que la muestra se puede ajustar con modelo Autorregresivo AR(p). Cuando el SPACF decae paulatinamente y el SACF contiene picos significativos hasta un valor de retardo τ=q, se considera que la muestra se puede ajustar con modelo de medias m´oviles MA(q). Cuando ambas muestras, SACF y SPACF, decaen paulatinamente, se considera que la muestra se puede ajustar con modelo Autorregresivo de medias m´oviles 4.4. REDES NEURONALES ARTIFICALES 75 ARMA(p,q). En este caso, pyqse deber´an aumentar gradualmente hasta obtener un modelo ´optimo. En algunas ocasiones es muy dif´ıcil decidir entre los diferentes modelos, por ejemplo entre un AR(6) o un ARMA(2,1). Es por ello, que se han desarrollado unas t´ecnicas de decisi´on basada en la teor´ıa bayesiana [Lauret12, Lebarbier04, LR12], donde se tiene en cuenta el error cometido por el modelo y la complejidad del mismo a la hora de decidir el orden. En esta tesis, se ha optado por el Criterio de Informaci´on Bayesiana, BIC por sus siglas en ingl´es (Bayesian Information Criterion), ya que penaliza m´as el n´umero de par´ametros del modelo que otras t´ecnicas existentes como el Akaike Information Criterion AIC, BIC =−2ln(ˆ L) + mkln(N) (4.14) donde ˆ Les el valor de m´axima verosimilitud del modelo, mkes el n´umero de par´ametros a definir del modelo, en nuestro caso depende de los ´ordenes pyq, y Nes el n´umero de valores de la serie temporal en la ecuaci´on (4.14). Para realizar una simplicaci´on en el c´alculo de la m´axima verosimilitud, se asume que la distribuci´on es gaussiana y por lo tanto la ecuaci´on finalmente utilizada para definir el BIC es (4.15), BIC =ln(ˆσ2) + mk ln(N) N(4.15) donde σ2es la varianza de los errores. Para tomar una decisi´on del modelo ´optimo a utilizar, se utilizar´an las gr´aficas de SPACF y SACF para identificar el n´umero m´aximo de orden de los modelos AR y MA respectivamente. Estos valores pyqno ser´an los valores finalmente utilizados para ajustar la serie, sino que se realizar´an los c´alculos del par´ametro BIC para todas las combinaciones posibles de modelos variando dichos par´ametros {1,2, ...p}y {1,2, ...q}. El modelo ´optimo a utilizar teniendo en cuenta la complejidad del mismo ser´a el que muestre un valor BIC menor. El cap´ıtulo 5, donde se muestran los resultados obtenidos para la Isla de Gran Canaria, se observar´a que finalmente, adem´as del estudio del Bayesian Information Criterion BIC, se estudiaron los valores del error cuadr´atico medio relativo de modelo de ajuste. En muchas ocasiones, el modelo ´optimo se˜nalado por el BIC no ofrec´ıa una mejora sustancial del error con respecto a modelo m´as sencillo, en el que se utilizaban un menor n´umero de valores pasados de la serie temporal. 4.4. Redes Neuronales Artificales En la mayor´ıa de los casos, la respuesta inmediata a cu´al es la diferencia m´as significativa entre el ser humano y el resto de los animales, ser´ıa sin lugar a dudas nuestra capacidad de raciocinio. Esta caracter´ıstica nos ha permitido desarrollar una tecnolog´ıa capaz de emular nuestras capacidades. En los ´ultimos a˜nos, con el desarrollo de la ciencia y la tecnolog´ıa, se ha manifestado la necesidad de realizar el tratamiento de gran cantidad de informaci´on, lo que ha conducido a enfrentarnos al reto de dise˜nar sistemas inteligentes 76 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR capaces de emular el funcionamiento del cerebro humano. Este es el objetivo de una rama cient´ıfica conocida como Inteligencia Artificial (IA). La Inteligencia Artificial se suele dividir en dos grandes ´areas, la IA Simb´olica y la IA Subsimb´olica [Isasi04]. En la primera, los sistemas que se utilizan para resolver el problema se dise˜nan siguiendo esquemas previamente fijados. En la IA se dice que estos sistemas siguen un esquema de arriba hacia abajo, ya que necesitan disponer de una soluci´on aproximada. Por el contrario, en la IA Subsimb´olica no se utilizan esquemas dise˜nados previamente, sino que se parte de sistemas gen´ericos que ir´an modific´andose y gener´andose mediante mecanismos de aprendizaje hasta formar un sistema capaz de resolver el problema. Es decir, desde esta perspectiva se estudian los mecanismos de los sistemas nerviosos del cerebro que nos hacen inteligentes, para poder dise˜nar sistemas basados en su estructura, funcionamiento y caracter´ısticas que se adapten a los problemas que se deseen resolver. Las t´ecnicas de Machine Learning, o de aprendizaje autom´atico en espa˜nol, se encuentran dentro del ´ambito de la Inteligencia Artificial Subsimb´olica. Estas t´ecnicas trabajan en el desarrollo y estudio de algoritmos capaces de resolver un problema o tomar decisiones a partir de un conjunto de datos conocido. Por ejemplo, para conseguir crear una hip´otesis y ofrecer una respuesta a una situaci´on no conocida, los algoritmos generan un modelo capaz de aprender de los datos de partida. El modelo generado se ajustar´a a la realidad en la medida en que pueda realizar predicciones de datos nuevos, que no se hayan presentado durante el entrenamiento. Por lo tanto, el conjunto de datos observados de un fen´omeno concreto se deber´a dividir en un conjunto de entrenamiento y otro de test, de manera que el primero nos permita generar un modelo y el segundo se utilice para validar el ajuste del mismo. La validaci´on del modelo se calcular´a estableciendo una desviaci´on sobre la realidad. Entre las t´ecnicas de Machine Learning m´as conocidas se encuentran las Redes Neuronales Artificiales (RNA), un sistema de aprendizaje basado en la estructura de las neuronas biol´ogicas, las Redes Bayesianas, en las que los nodos de la red se encuentran asociados a una funci´on de probabilidad, los Support Vector Machines, utilizados para resolver problemas de clasificaci´on y regresi´on, los Gaussian Process, ´utiles para el estudio de an´alisis espacial, etc. Aunque todos estos m´etodos han sido utilizados satisfactoriamente en el campo de la radiaci´on solar [Diagne13, Mellit08], en esta tesis nos hemos centrado en el trabajo con Redes Neuronales Artificiales y Redes Bayesianas. 4.4.1. Introducci´on a las redes neuronales En general las redes neuronales tratan de emular el funcionamiento de una neurona como el elemento m´as simple en la estructura del cerebro. En la bibliograf´ıa especializada nos encontramos con diferentes definiciones, como la que realiza Mackay ”Muchos investigadores les gustar´ıa crear m´aquinas que puedan aprender, reproducir patrones o descubrir patrones en los datos” [MacKay03]. En la mayor´ıa de las definiciones se establece una analog´ıa con las neuronas del cerebro humano, aunque, como se ver´a m´as adelante, no todas las RNA emulan una 4.4. REDES NEURONALES ARTIFICALES 77 determinada estructura neuronal. Todas la redes neuronales realizan operaciones en una serie de elementos b´asicos que se conocen como neuronas por analog´ıa. Estas unidades o neuronas est´an interconectadas entre s´ı por unos pesos sin´apticos, lo cuales variar´an con el tiempo durante el periodo de aprendizaje. En general siempre que se describa una RNA se especificar´an los siguientes tres aspectos [MacKay03]: Arquitectura.- Se debe tener en cuenta el n´umero de neuronas involucradas y la forma de relaci´on entre cada una de ellas. Regla de activaci´on.- En esta parte se define en qu´e manera las neuronas responden a la interconexiones entre ellas. M´etodo de entrenamiento.- Los m´etodos de entrenamiento definen la manera en la que los pesos cambian durante el mismo. Aunque se pueden situar las primeras investigaciones sobre Redes Neuronales Artificiales en el siglo XIX, no es hasta la segunda mitad del siglo XX en que se comienza realmente a desarrollar esta herramienta, gracias entre otras cosas al desarrollo del hardware. El primero que intent´o realizar avances en computaci´on a partir del estudiar el cerebro fue Alan Turing en 1936. Warren McCulloch y Walter Pitts desarrollaron en 1943 una teor´ıa sobre c´omo trabajar con neuronas modelando una red simple a partir de circuitos el´ectricos [McCulloch43, Gonz´alez95]. Posteriormente, en 1957 Frank Rosenblatt comenz´o el desarrollo de uno de los tipos de redes que utilizamos en esta tesis, llamado Perceptr´on.´ Este era capaz de reconocer una serie de patrones que se le hubieran presentado anteriormente, pero presentaba una serie de limitaciones a la hora de clasificar clases no separables linealmente [Rosenblatt58]. En 1960 Widrow y Hoff presentaron el ADALINE, Adaptive Linear Element, el cual, a trav´es de un algoritmo de entrenamiento sencillo denominado LMS (Least Mean Square), consegu´ıa un sistema adaptativo que aprend´ıa de forma m´as precisa que el perceptr´on [Widrow60, Haykin96]. Los inconvenientes que presentaba el perceptr´on a la hora de resolver problemas no separables linealmente, llevaron a Minsky y Papert en 1969 a publicar su libro Perceptrons [Minsky88]. En ´este se expon´ıan con mucha claridad los problemas que era capaz de resolver el perceptr´on, pero se incid´ıa en los inconvenientes que presentaban en la mayor´ıa de los casos. Adem´as tambi´en expusieron sus opiniones en contra del uso de las extensiones del perceptr´on, como el multicapa, hecho por el que fueron criticados posteriormente ya que sus conjeturas parecen haber sido err´oneas. En cualquier caso, este trabajo detuvo el desarrollo de las Redes Neuronales Artificiales durante diez a˜nos. Durante estos a˜nos, se sucedieron algunas investigaciones por parte de cient´ıficos diversos como Kohonen y Aderson, que proponen independientemente un modelo similar de memoria asociativa, Asociador Lineal. En 1982, con la publicaci´on de John Hopfield [Hopfield82], se recuper´o el inter´es por este tipo de herramientas. En este trabajo, el autor describe de manera muy clara una variante del Asociador Lineal. As´ı mismo, Hopfield mostr´o c´omo esta variante es capaz de trabajar. El uso principal que se le ha dado a estas redes es como memorias y 78 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR para resolver problemas de optimizaci´on. En ese mismo a˜no, Teuvo Kohonen [Kohonen88] publica un trabajo sobre mapas autoorganizativos mediante reglas simples, cuyo sistema de aprendizaje es de tipo no supervisado, mientras que al a˜no siguiente, Fukushima, Miyake e Ito presentan un dispositivo capaz de realizar reconocimiento de patrones con ´exito, el Neocognitr´on [Fukushima83]. En 1986, se presenta el trabajo de Rumelhart, Hinton y Williams, en el que se desarrolla un algoritmo de aprendizaje conocido como retropropagaci´on (backpropagation) [Rumelhart86] para redes neuronales multicapa. A partir de esta publicaci´on, el n´umero de trabajos sobre RNA se han multiplicado apareciendo un gran n´umero de aportaciones en los m´etodos de aprendizaje y tipo de estructuras. 4.4.2. Aplicaciones Las RNA tienen m´ultiples aplicaciones que se pueden dividir seg´un el problema a resolver y el campo de conocimiento en el que se aplican. As´ı, para resolver problemas de Clasificaci´on de elementos nos encontramos los siguientes trabajos: Medicina.- En el campo del diagn´ostico m´edico, se han desarrollado trabajos en detecci´on de cardiopat´ıas o en detecci´on de tumores cancer´ıgenos, en la que una red neuronal analiza la posible presencia de tumores en una imagen. Farmacia.- Se utilizan para diagnosticar posibles efectos adversos al administrar un determinado f´armaco para el tratamiento de cardiopat´ıas o tratamientos oncol´ogicos. Procesado de se˜nal.- En este campo las redes neuronales se han desarrollado en un amplio espectro. Por ejemplo en la ecualizaci´on de canales de comunicaci´on se ha mostrado m´as efectiva que otros m´etodos. Se emplea tambi´en en el reconocimiento de patrones en im´agenes o en el reconocimiento de voz. Econom´ıa.- Debido a la necesidad de tomar decisiones entre una gran cantidad de opciones las RNA son bastante aplicables en la concesi´on de cr´editos, detecci´on de posibles fraudes o la posibilidad de la quiebra de un banco. Por otro lado tambi´en se pueden aplicar para resolver problemas de Modelizaci´on. Se pueden encontrar diversos trabajos con RNAs en los siguientes campos: Medicina.- Se relacionan con el modelado de diferentes tipo de se˜nales como el electrocardiograma (FEGG), el electromiograma (EMG), electroencefalograma (EEG), etc. Farmacia.- Determinar la concentraci´on de un determinado f´armaco en sangre. Procesado de se˜nal.- En este campo se han hecho trabajos en la eliminaci´on activa de ruido y en el control de sistemas. 4.4. REDES NEURONALES ARTIFICALES 79 Econom´ıa.- Las RNA has sido utilizadas para intentar predecir comportamientos futuros que permitan valorar el ´exito de una operaci´on. En concreto encontramos trabajos en predicci´on del gasto el´ectrico, cambio de moneda, tendencias a corto y medio plazo en la bolsa, predicci´on de stocks, etc. Medio Ambiente.- Las variaciones en el medio ambiente no son lineales y dependen de muchas variables por lo que lo hacen un campo muy apetecible para utilizar las RNA. En concreto, en la predicci´on de la radiaci´on solar que nos ocupa en esta tesis, se encuentran multitud de trabajos al respecto. Estos trabajos se han desarrollado partiendo de datos de radiaci´on como variable de entrada [Hontoria01, Lauret15] o de la mezcla ´estos con otros par´ametros meteorol´ogicos como la temperatura, humedad relativa, horas de sol, velocidad y direcci´on del viento [AA98, Rehman08, Mellit10, Ghanbarzadeh09]. Por otro lado, tambi´en se han aplicado en la predicci´on de niveles t´oxicos de ozono en zonas urbanas y rurales o en la predicci´on de temperatura [Almonacid13] As´ı, el objetivo de las Redes Neuronales Artificiales es crear un sistema con una serie de unidades, de manera que el comportamiento global trate de emular el funcionamiento del sistema neuronal humano. Esto hace imprescindible el estudio de las caracter´ısticas y funcionamiento de una neurona simple. 4.4.3. Fundamentos te´oricos de las redes neuronales El sistema nervioso humano es el encargado de realizar las comunicaciones entre los ´organos de los sentidos, que reciben los est´ımulos del exterior, y los ´organos diana, que realizan las acciones (m´usculos, gl´andulas). As´ı pues, se encarga de recoger la informaci´on del exterior, trasmitirla, procesarla y enviarla una vez elaboradas. En este proceso nos encontramos con los receptores, c´elulas encargadas de recoger la informaci´on del exterior o del propio interior del cuerpo, el sistema nervioso, que recoge, elabora y env´ıa esta informaci´on, y los ´organos diana, que reciben la informaci´on y la interpretan. La transmisi´on de toda esta informaci´on se realiza por medio de unas c´elulas espec´ıficas conocidas como neuronas. Estas c´elulas forman redes en las cuales se elabora y transmite la informaci´on mediante se˜nales electro-qu´ımicas. Por lo tanto, una parte de la red de c´elulas neuronales debe estar conectada a los ´organos receptores, donde reciben la informaci´on del exterior o del interior del cuerpo. La informaci´on recibida de los receptores se denomina est´ımulo y pone en funcionamiento la red. Las neuronas de la red conducen la informaci´on hasta las conexiones con los ´organos diana. La red neuronal no es m´as que multitud de neuronas simples interconectadas entre s´ı. Por lo tanto, una neurona puede recibir informaci´on de un ´organo receptor o de otra neurona, y la transmite a otra neurona o a un ´organo diana [Isasi04, Gonz´alez95]. Los componentes principales de una c´elula neuronal son los siguientes: Dendritas, la neurona recibe la informaci´on del exterior y la transmiten al interior de la c´elula a trav´es de estas ramificaciones de entrada. 86 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR Figura 4.6: Ejemplo de los elementos de una red multicapa general aj=fj[ Ne X i=1 (ω1 jixi)] (4.25) donde ω1representa una matriz con todos los pesos sin´apticos de conexi´on entre la capa de entrada y la primera capa oculta. Esta matriz contendr´a en cada fila todos los pesos de conexi´on entre cada una de la entradas y una neurona espec´ıfica de la capa oculta. Por contra, en cada columna nos encontramos la conexi´on entra una entrada cualquiera y cada una de las neuronas de la capa oculta (4.26). ω1=      ω1 11 ω1 12 ... ω1 1Ne ω1 21 ω1 22 ... ω1 2Ne ... ... ... ... ω1 j1... ω1 ji ω1 jNe ... ... ... ...      (4.26) La funci´on de activaci´on de cada neurona de la capa oculta est´a representada por fj, que se activa si la suma ponderada de las entradas supera el umbral. Las salidas intermedias ajse propagan hacia la siguiente capa oculta o hacia la capa de salida de la misma manera que entre las capas anteriores. Las conexiones entre ellas se realizar´a mediante los pesos ω2. As´ı en la capa de salida final, la neurona de sproduce una salida seg´un la ecuaci´on (4.27), 4.4. REDES NEURONALES ARTIFICALES 87 Figura 4.7: Arquitectura del Perceptr´on Multicapa (MLP) ys=fs[ H X j=1 ω2 sjfj[ Ne X i=1 (ω1 jixi)] (4.27) donde Hmuestra el n´umero de neuronas de la capa intermedia y Nemuestra el n´umero de entradas a la RNA. En este caso, ω2 sj representa el peso de conexi´on de la neurona de la capa intermedia jcon la neurona de la capa de salida s, donde el super´ındice dos indica que nos encontramos en el segundo conjunto de pesos. La funci´on fses la funci´on de activaci´on de la neurona de salida s. El Perceptr´on simple mostraba limitaciones para resolver problemas de separaciones no lineales. Por este motivo, se desarroll´o el Perceptr´on Multicapa, en ingl´es Multilayer Perceptron (MLP), como una combinaci´on de varias unidades simples. Algunos autores han demostrado que el MLP es un aproximador universal de cualquier funci´on continua en el espacio Rn[Hornik89]. Al igual en que se explic´o para las RNA multicapa, el MLP dispone las neuronas entre la capa de entrada, al menos una capa oculta o intermedia y la capa final de salida. En la Figura 4.7 se representa un MLP con una capa oculta y una salida con una ´unica neurona. Las entradas se representan por las variables xi, los pesos ω1son los pesos sin´apticos entre la capa de entrada y la capa intermedia, mientras que ω2 son las conexiones entre la capa intermedia y la ´unica neurona de salida, y por ´ultimo la salida final de la red se representa por y. Las conexiones del MLP siempre son hacia adelante, feed-forward, es decir que cada capa se conecta con la capa siguiente. Generalmente, la red est´a conectada totalmente, con lo que todas las neuronas de una capa est´an conectadas con la todas las neuronas de la capa siguiente. Por otro lado, cada neurona tiene asociado un umbral de activaci´on, 88 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR que en el caso del MLP se suele tratar como una conexi´on m´as asociada a una entrada de valor la unidad. Las funciones de activaci´on fde cada neurona m´as utilizadas son la funci´on sigmoidal (4.18), con un rango de valores continuo entre [0,1] y la funci´on tangente hiperb´olica (4.19), con un rango de valores continuo entre [−1,1]. El dise˜nador elegir´a el tipo de funci´on de activaci´on bas´andose en los valores que desee que tomen las salidas de las neuronas. Por contra, las neuronas de salida tienen una funci´on de activaci´on distinta y, es muy com´un que ´esta sea la funci´on lineal (4.17). Matem´aticamente, un MLP con una capa intermedia con funci´on de activaci´on fjen la neurona jy una ´unica neurona de salida con una funci´on de activaci´on lineal se puede representar con la siguiente formulaci´on. Las salidas de la capa intermedia se obtienen seg´un la expresi´on (4.28). aj=fj[ Ne X i=1 (ω1 jixi) + ω1 0] (4.28) Las salidas ajde la capa intermedia se propagar´an hacia la ´unica neurona de salida y producir´an la salida yseg´un la ecuaci´on (4.29). y= H X j=1 ω2 sjfj[ Ne X i=1 (ω1 jixi) + ω1 0] + ω2 0(4.29) Cuando se utiliza un MLP para resolver un problema, uno de los pasos m´as importantes a realizar es la elecci´on de la arquitectura del mismo. Este dise˜no implica la elecci´on del n´umero de neuronas de cada capa y el n´umero de capas. Algunos de estos par´ametros los debe seleccionar el dise˜nador y otros vienen impuestos por la propia naturaleza del problema. El n´umero de neuronas de la capa de salida viene definido por las variables que se est´an trabajando. En el caso de esta tesis, la salida que se pretende es ´unica porque se busca la predicci´on de la radiaci´on solar. Sin embargo, el n´umero de entradas no siempre viene definido por el problema. Se puede dar el caso, de que existan algunas entradas que no sean relevantes para la soluci´on del problema y la informaci´on que aportan es irrelevante, por lo que su uso complica la red neuronal sin obtener ning´un beneficio. Para evitar estas situaciones, es muy conveniente hacer un estudio de la importancia de cada una de las entradas para resolver el problema y descartar previamente las que no sean relevantes. Existen multitud de t´ecnicas para realizar este tipo de an´alisis (algoritmos gen´eticos, an´alisis de sensibilidad, etc). Para definir el n´umero de capas ocultas y el n´umero de neuronas de cada una de ellas no existe ninguna regla te´orica definida. As´ı, el dise˜nador realizar´a la elecci´on bas´andose en t´ecnicas de ensayo y error. El n´umero de neuronas ocultas debe ser suficiente para que exista el n´umero de pesos adecuado para resolver el problema. Sin embargo, un n´umero elevado de conexiones complicar´ıa la estructura de la red, elevar´ıa la carga computacional derivada del proceso de entrenamiento e impedir´ıa la generalizaci´on de la soluci´on, ya que se ajustar´ıa demasiado a los datos que se le muestran en el entrenamiento pero no a otro conjunto de datos diferentes. El n´umero de neuronas 4.4. REDES NEURONALES ARTIFICALES 89 intermedias puede influir en el comportamiento pero, generalmente, no es un par´ametro significativo, ya que el MLP puede mostrar muy diferentes arquitecturas para resolver un mismo problema. En general con una sola capa intermedia el MLP puede resolver la mayor´ıa de los problemas, por lo que en esta tesis ´unicamente se estudiar´a el n´umero de neuronas de una ´unica capa intermedia. 4.4.6. Regla de aprendizaje. Backpropagation Con todo lo explicado, la importancia de las RNA radica en su capacidad de aprendizaje conforme se le muestran las relaciones entre los datos de entrada y salida. As´ı, la parte m´as importante del proceso de resoluci´on del problema es el aprendizaje. ´ Este consiste en general en la modificaci´on paulatina de los pesos sin´apticos de las conexiones entre neuronas, de manera que se capacite a la red para resolver de la mejor manera posible el problema. Las RNA son sistemas que aprenden a partir de ejemplos, por lo que ´estos deben ser suficientes en n´umero, para darle tiempo a la red a adaptar los pesos, y representativos, ya que si los ejemplos mostrados son espec´ıficos de un tipo de evento la red se especializar´a en resolver este tipo de problema. A la hora de preparar el conjunto de datos de entrenamiento es importante que todas las regiones significativas del espacio est´en representadas en el mismo. El proceso de aprendizaje de una red neuronal presenta una serie de pasos, Figura 4.8. Primero se establece el conjunto de datos de entrenamiento, conformado por datos de entrada xiy datos de salidas esperadas t, posteriormente se realiza un inicio aleatorio de los pesos de la RNA, se calcula la salida y se compara con la salida deseada, modificando los pesos seg´un la diferencia que exista hasta conseguir un determinado criterio de convergencia. Mientras no se cumpla dicho criterio, se repetir´a el proceso introduciendo todos los ejemplos del conjunto de entrenamiento. Cuando se cumple el criterio se dice que ha concluido el proceso de aprendizaje. La modificaci´on de los pesos puede realizarse despu´es de introducir cada ejemplo del conjunto, aprendizaje on-line, o una vez introducidos todos los ejemplos, aprendizaje Batch. El criterio de convergencia depender´a del tipo de problema a resolver y del tipo de estructura elegido. La finalizaci´on del proceso de aprendizaje se puede determinar: Mediante un n´umero de ciclos fijos, se fija el n´umero de iteraciones que se van a realizar de antemano. Cuando se alcance se detiene el aprendizaje y se acepta la RNA obtenida. Cuando el error descienda por debajo de una cantidad, se define antes de comenzar un error m´ınimo que se desea alcanzar. Puede darse el caso de que la red no consiga reducir el error hasta dicho objetivo, en ese caso se debe introducir otro criterio de parada (como puede ser el n´umero de ciclos). Si ocurre esto se dice que no se ha conseguido encontrar la soluci´on y deber´an cambiarse algunos par´ametros para probar otra vez. Cuando la modificaci´on del error sea insignificante en cada iteraci´on, puede ocurrir que, independientemente del valor nominal del error, ´este no se disminuya de manera 90 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR Figura 4.8: Esquema de aprendizaje supervisado de una RNA significativa con cada iteraci´on. Se fija de antemano un valor m´ınimo de modificaci´on del error entre una iteraci´on y la siguiente. Cuando la modificaci´on de los pesos sea irrelevante, si en el proceso de entrenamiento llega el momento en que ya no se var´ıan los valores de los pesos en cada iteraci´on, se da por finalizado el entrenamiento. El m´etodo de aprendizaje es el proceso seguido para ir adaptando los par´ametros de la red. En esta tesis el criterio que se va a utilizar se basa en el estudio del error cometido por la red. Existen diversos m´etodos para minimizar la funci´on de error elegida. En el caso del Perceptr´on Multicapa (MLP), los m´as utilizados son los algoritmos por descenso de gradiente que se basan en la minimizaci´on de una determinada funci´on de error entre la salida deseada ty la obtenida en la predicci´on y. La funci´on de error EDampliamente utilizada en la bibliograf´ıa es el error cuadr´atico medio [Bishop95], seg´un la ecuaci´on (4.30), ED(ω) = 1 2 N X i=1 (ti−yi)2(4.30) donde Nes el n´umero total de patrones del conjunto de datos de entrenamiento. Cada patr´on se considera una pareja de valores entrada-salida relacionados entre s´ı {xi, ti}. El algoritmo de aprendizaje que utilizar´a el MLP es el de RetroPropagaci´on, conocido por su nombre en ingl´es Backpropagation. En su caso, es un algoritmo de descenso por gradiente que actualiza los pesos en cada iteraci´on retropropagando la se˜nal desde la capa salida hasta la capa de entrada. El algoritmo primero propaga la se˜nal hacia adelante, desde 4.4. REDES NEURONALES ARTIFICALES 91 Figura 4.9: Esquema del avance en la optimizaci´on. En la izquierda se toma un valor alto de coeficiente de aprendizaje y en la derecha valores peque˜nos. la entrada a la salida, calculando la salida ya partir de los conjuntos de pesos actuales (ω1, ω2) y estimando el error entre esta salida obtenida y la deseada t. Posteriormente, en funci´on del error obtenido en la salida, el algoritmo actualiza los valores de los pesos sin´apticos que determinan las conexiones entre las neuronas partiendo de la capa de salida a la capa de entrada. La minimizaci´on de la funci´on de error es un problema no lineal y, como consecuencia, tienen que utilizarse t´ecnicas de optimizaci´on no lineales para su resoluci´on. En el caso del Perceptr´on la regla de minimizaci´on es ajustar los par´ametros de la red siguiendo la direcci´on negativa del gradiente de la funci´on de error. Por tanto, aplicando el m´etodo de descenso de gradiente, cada peso ωde la RNA se modifica en cada iteraci´on con la siguiente regla de aprendizaje general (4.31), ωk+1 =ωk−η∂E ∂ωk (4.31) donde kyk+ 1 representan las sucesivas iteraciones durante el entrenamiento. Se debe tener en cuenta que en cada entrenamiento los pesos se modifican y se vuelve a calcular el error. El par´ametro ηse conoce como el coeficiente de aprendizaje e influir´a en la velocidad de convergencia durante el entrenamiento. Si se elije un valor de ηmuy alto, en cada iteraci´on se produce una modificaci´on grande de los pesos y se desplazar´a r´apidamente por la superficie de la funci´on de error. De esta manera se corre el riesgo de pasar por encima de un punto m´ınimo y oscilar alrededor de ´el sin alcanzarlo. Por el contrario, si se elije un valor muy peque˜no las variaciones de pesos ser´an muy peque˜nas y el descenso ser´a muy paulatino, por que ser´an necesarias muchas iteraciones, ver Figura 4.9. La actualizaci´on del los pesos en el algoritmo de backpropagation se inicia por la capa de salida y se va propagando hasta la capa de entrada. En esta tesis, se trabajar´a en todos los casos con redes de una sola capa intermedia con funci´on tangente hiperb´olica y una 92 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR salida ´unica con funci´on lineal. La regla general de aprendizaje para este tipo de redes conlleva pues, la actualizaci´on del conjunto de pesos que une la capa oculta con la de salida, ω2, y la los pesos que unen la capa de entrada con la capa oculta, ω1. Teniendo en cuenta estos aspectos, la actualizaci´on de los pesos de la capa de salida queda definida seg´un (4.32). ω1j(k+ 1) = ω1j(k)−η∂ED ∂ω1j(k)(4.32) Teniendo en cuenta la funci´on de error (4.30) y que la salida deseada tno depende de los pesos y por lo tanto es una constante, la derivada del error respecto de los pesos de salida queda seg´un la ecuaci´on (4.33). ∂ED ∂ω1j(k)=−(t(k)−y(k)) ∂y(k) ∂ω1j(k)(4.33) Llegados a este punto se deber´a calcular la derivada de la funci´on de salida yrespecto del peso ω1j. Como se explic´o, la funci´on de salida yqueda definida seg´un la ecuaci´on (4.29). Al aplicar la regla de la cadena se debe tener en cuenta que el ´unico t´ermino de la ecuaci´on (4.29) cuya derivada es distinta de cero es ω1jaj. El t´ermino ajest´a definido por la ecuaci´on (4.25), por lo que La derivada de la salida yrespecto de los pesos quedar´a seg´un la ecuaci´on (4.34). ∂y(k) ∂ω1j(k)=fj[ N X i=1 (ω1 jixi) + ω1 0] = aj(4.34) Siendo com´un definir δcomo el t´ermino de actualizaci´on de la salida en cuesti´on para la iteraci´on de entrenamiento kseg´un la ecuaci´on (4.35) δ(k) = −(t(k)−y(k)) (4.35) Finalmente se concluye que la actualizaci´on de los pesos del conjunto de salida se realiza seg´un la ecuaci´on (4.36). ω1j(k+ 1) = ω1j(k) + η(t(k)−y(k))aj=ω1j(k) + ηδ(k)aj(4.36) Se puede observar que la actualizaci´on del peso de conexi´on entre la neurona jde la capa oculta y la salida ´unicamente depender´a de la salida intermedia ajde dicha neurona y del error cometido por la salida obtenida yrespecto a la salida deseada t, expresado en la funci´on δ(k). Siguiendo el mismo m´etodo se puede obtener la expresi´on de actualizaci´on del umbral de la neurona de salida ω2 0seg´un la ecuaci´on (4.37). ω2 0(k+ 1) = ω2 0(k) + ηδ(k) (4.37) Una vez se han actualizado los pesos de la capa de salida, se continua con la retropropagaci´on y se realizar´a la actualizaci´on de los pesos que unen las capa de entrada 4.4. REDES NEURONALES ARTIFICALES 93 con la capa oculta, ω1. As´ı, la actualizaci´on de la conexi´on entre la entrada iy la neurona intermedia jqueda definida seg´un (4.38). ωji(k+ 1) = ωji(k)−η∂E ∂ωji(k)(4.38) Al igual que en el caso anterior se utiliza la regla de la cadena para hallar las expresiones generales de la derivada. Se debe tener en cuenta que el peso ωji influye en la activaci´on de la neurona j,aj, y que el resto de activaciones de la capa oculta no dependen de dicho peso, ωji(k+ 1) = ωji(k)−ηδj(k)xi(4.39) donde xies la entrada correspondiente al peso que se est´a actualizando y δj(k) viene dado por la expresi´on (4.40) para la neurona jen la iteraci´on k, δj(k) = f0 j[ N X i=1 (ω1 jixi) + ω1 0]δ(k)ω1j(4.40) donde fjes la funci´on de activaci´on de la neurona intermedia j, cuya expresi´on es la tangente hiperb´olica en el caso que estamos explicando. Por lo que su derivada se puede expresar de forma sencilla seg´un (4.41). f0 j= 1 −f2 j(4.41) Por ´ultimo, de la misma manera que en el caso anterior, se obtiene la regla de actualizaci´on de los umbrales de las neuronas de la capa oculta. ω1 0(k+ 1) = ω1 0(k) + ηδj(k) (4.42) El m´etodo del descenso de gradiente, para un MLP como el descrito, es un algoritmo simple de optimizaci´on pero presenta una serie de inconvenientes. Se puede dar el caso del m´etodo alcance un m´ınimo local y por lo tanto no consiga minimizar la funci´on hasta el valor global. Por otro lado, en las zonas donde la pendiente de la funci´on de error sea casi nula, el n´umero de iteraciones necesario para alcanzar el m´ınimo podr´ıa ser exageradamente grande. Para resolver estos problemas se han presentado numerosas variantes y modificaciones del algoritmo de entrenamiento. Una de la variantes m´as comunes es aplicar un factor de momentum o factor de inercia. En este caso, se le a˜nade a la funci´on de actualizaci´on un nuevo t´ermino seg´un la ecuaci´on (4.43). ωk+1 =ωk−η∂E ∂ωk +βm(ω(k)−ω(k+ 1)) (4.43) En general, la adici´on de este par´ametro disminuye el n´umero de iteraciones necesarias para alcanzar la soluci´on y tiende a evitar los m´ınimos locales. De esta manera, cuando 94 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR el incremento de pesos es alto, en la siguiente iteraci´on lo ser´a a´un m´as. Por otro lado, en caso de que los incrementos oscilen entre positivo y negativo, el incremento efectivo final se reduce. En los ´ultimos a˜nos [Bishop95, BS10], para evitar los problemas derivados de este tipo m´etodos, se han desarrollado diversos m´etodos de optimizaci´on de la funci´on de error. En esta tesis se utilizar´a el algoritmo de gradiente conjugado escalado, considerado como un algoritmo de segundo orden en cuya formulaci´on interviene la matriz Hessiana. Son m´etodos m´as robustos que aceleran la velocidad de convergencia, aunque aumentan la complejidad del modelo y la carga computacional. 4.4.7. Algoritmo de gradiente conjugado escalado El algoritmo de gradiente conjugado escalado, en ingl´es scale conjugate gradient, es tambi´en un algoritmo de aprendizaje supervisado cuya principal ventaja es que converge m´as r´apido que el m´etodo de descenso del gradiente [Bishop95]. La diferencia m´as notable con respecto a este ´ultimo es que no utiliza una ´unica direcci´on para minimizar la funci´on ni un valor constante del coeficiente de aprendizaje η. En cada iteraci´on se busca la direcci´on ´optima entre un conjunto de direcciones conjugadas y se calcula el valor ´optimo del coeficiente de aprendizaje en esa direcci´on [BS10]. En general, el proceso de entrenamiento para la iteraci´on kse inicia con la definici´on aleatoria del primer conjunto de pesos ω(k) de las conexiones de la RNA y se eval´ua el gradiente seg´un dichos par´ametros, g(k) = ~ ∇E(4.44) y se toma la direcci´on inicial d(k) de b´usqueda seg´un el mismo criterio que el algoritmo de descenso de gradiente, ~ d(k) = −~ g(k) (4.45) y el coeficiente de aprendizaje ηqueda definido seg´un la ecuaci´on (4.46), η(k) = −~ d(k)t~ g(k) ~ d(k)tH(k)~ d(k)(4.46) donde H(k) representa la matriz hessiana en la iteraci´on kdefinida seg´un la expresi´on (4.47). H(k) = ∂2E ∂ωi∂ωj (4.47) As´ı, en cada iteraci´on se debe calcular y almacenar los valores correspondientes a la matriz hessiana, lo que conlleva mayor complejidad de c´alculo. Una vez se dispone de estos par´ametros se procede a actualizar los pesos seg´un la ecuaci´on (4.31). Con los nuevos pesos se estudia el resultado de la RNA con el grupo de pesos ω(k+ 1) y se comprueba que se haya alcanzado el criterio de convergencia fijado. Si no se alcanza este criterio se deber´a 4.4. REDES NEURONALES ARTIFICALES 95 realizar el c´alculo del nuevo gradiente con los nuevos par´ametros ~g(k+ 1) y evaluar la nueva direcci´on de b´usqueda usando las siguientes ecuaciones. ~ d(k+ 1) = −~g(k+ 1) + βm(k)~ d(k) (4.48) βm(k) = −~gt(k+ 1)(~g(k+ 1) −~g(k)) ~ gt(k)~ g(k)(4.49) Este proceso se repetir´a en cada iteraci´on hasta alcanzar el criterio de convergencia. En esta tesis se ha utilizado este algoritmo seg´un se ha implementado en la herramienta NETLAB en el entorno MATLAB [Bishop95, Nabney02]. 4.4.8. T´ecnicas de optimizaci´on de la arquitetura de la red La elecci´on de la red m´as sencilla posible capaz de resolver el problema es una de las decisiones m´as importantes, ya que ayuda a simplificar los c´alculos y permite una mejor generalizaci´on del resultado. Algunos de los m´etodos m´as utilizados son los m´etodos de poda, que consisten en la eliminaci´on de pesos o neuronas innecesarios durante el aprendizaje de la red. El algoritmo consistir´a en realizar una elecci´on de una arquitectura m´as compleja de la necesaria y que sea el algoritmo el que se encargue de podar la conexiones innecesarias. Adem´as de permitir reducir la carga computacional, nos ayudar´a a encontrar posibles entradas que no aporten nada al c´alculo de la salida, ya que los pesos correspondientes a dichas entradas ser´an podados. Los m´etodos basados en aplicar t´erminos de penalizaci´on son los m´etodos m´as simples y utilizados en la pr´actica. Se basan en incluir un t´ermino adicional de regularizaci´on Eωen la ecuaci´on general del error (4.30), con el objetivo de eliminar los pesos innecesarios. Este t´ermino provoca un decaimiento en estos pesos de manera que se puede optar por eliminar aquellos que, al final del entrenamiento, queden por debajo de un umbral determinado. S(ω) = ED(ω) + µEω(ω) (4.50) Eω(ω) = 1 2 m X j=1 (ω2 j) (4.51) donde mrepresenta el n´umero de par´ametros de la red y µes el par´ametro de regularizaci´on o penalizaci´on, que determina la importancia de los t´erminos de regularizaci´on. Utilizar el valor correcto de este par´ametro es unos de los aspectos importantes a la hora de definir una RNA. Un valor peque˜no del mismo puede conducir a un sobreajuste, en ingl´es overfitting, con lo que la funci´on representar´ıa correctamente los datos del conjunto de entrenamiento pero no ser´ıa capaz de ajustar un conjunto no conocido de datos, es decir, no se conseguir´ıa la generalizaci´on de la soluci´on, Figura 4.10. Mientras que si el valor es muy grande puede conducir a no ajustar los datos, en ingl´es underfitting. 102 CAP´ ITULO 4. MODELOS DE PREDICCI ´ ON DE RADIACI ´ ON SOLAR Figura 4.12: Resultados de la t´ecnica de log of evidence para decidir el n´umero de neuronas ocultas ´optimo para resolver un problema. Ejemplo de la estaci´on C6 en Gran Canaria. la entrada correspondiente al conjunto de pesos regido por el hiperpar´ametro en cuesti´on se considera no relevante para la red y puede ser eliminada. En la pr´actica se muestran figuras de barras como la Figura 4.13, donde cada barra muestra el valor de la varianza de los hiperpar´ametros, es decir la inversa del valor del hiperpar´ametro αg. Las entradas que corresponden con una barra que muestre un valor bajo con respecto al resto de entradas se considera irrelevante y se puede probar a realizar una poda de dicha entrada. Esta t´ecnica muestra una estimaci´on de la probabilidad de la importancia de cada entrada. As´ı que se recomienda comprobar si el error disminuye al realizar una poda de las entradas se˜naladas como poco relevantes. Figura 4.13: Resultados de la t´ecnica de selecci´on autom´atica del n´umero de entradas necesarias para resolver un problema. Ejemplo de la estaci´on C5 en Gran Canaria. Cap´ıtulo 5 Aplicaci´on de los modelos de predicci´on de Radiaci´on Solar 5.1. Introducci´on En este cap´ıtulo se expondr´an los resultados obtenidos en la predicci´on de los valores de radiaci´on solar horario para la isla de Gran Canaria. El objetivo de esta tesis es el estudio de distintos m´etodos de predicci´on de la radiaci´on solar global horizontal con horizontes temporales desde h= 1 hasta h= 6 horas. Los modelos de predicci´on elegidos para realizar este estudio son, como ya se ha comentado en el Cap´ıtulo 4, distintos m´etodos estad´ısticos. En un primer momento se estudiar´a la complejidad de dichos modelos para obtener, en cada caso, una soluci´on ´optima adaptada a nuestros datos. Se ha realizado el estudio de la complejidad de los modelos lineales Autorregresivos (AR), Autorregresivos de Medias M´oviles (ARMA) y de las Redes Neuronales Artificiales (RNAs). Este estudio se realiz´o tomando como datos de partida los datos hist´oricos de radiaci´on solar recogidos por las estaciones de medida terrestres de la isla de Gran Canaria. Posteriormente, se estudiar´a si las predicciones de radiaci´on solar mejoran aportando diferentes datos de partida. De esta manera, adem´as de los datos hist´oricos de radiaci´on solar, se expondr´an simulaciones de predicci´on a partir de datos de radiaci´on satelitales, datos de radiaci´on del modelo de predicci´on num´erico ECMWF y datos de viento y humedad relativa del ECMWF. En este caso, el modelo estad´ıstico utilizado ser´an las RNAs ya que aporta los mejores resultados de predicci´on. Las predicciones de radiaci´on solar a partir de m´etodos estad´ısticos necesitan disponer de suficientes datos hist´oricos para realizar el entrenamiento de los modelos. En esta tesis se ha optado por trabajar con un a˜no de medidas para realizar el entrenamiento y otro para comprobar el ajuste del modelo obtenido, test. Los errores finales que se utilizar´an para comprobar la bondad de los modelos se obtendr´an realizando las predicciones con el a˜no del test. Con esto se pretende evitar el sobreajuste (overfitting) de los modelos a unos datos en concreto. Los a˜nos de medida de radiaci´on solar, tanto para el entrenamiento como para el test, deben tener el menor porcentaje de huecos de medidas posible, ya 105 106 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON que se pretende tener un conjunto de datos continuo que represente todas las diferencias clim´aticas anuales. En todas las estaciones de medida terrestres disponibles en la isla de Gran Canaria se pueden establecer dos a˜nos de media completos para trabajar con los modelos estad´ısticos. El modelo num´erico de predicci´on meteorol´ogico ECMWF tambi´en ofrece la posibilidad de contar con los dos mismos a˜nos de medida en cada estaci´on. Como se comenta en el Cap´ıtulo 2, entre los a˜nos 2000 y 2005 el ECMWF ´unicamente puede ofrecer una resoluci´on temporal de 3 horas. As´ı, se disponen de las predicciones realizadas por el modelo num´erico de predicci´on ECMWF para la isla de Gran Canaria cada tres horas para el d´ıa siguiente (horizonte de un d´ıa de adelanto). Por otro lado, los datos satelitales se obtuvieron de la base de datos del Helioclim3, en particular, de la versi´on 5 (HC3v5). La base de datos Helioclim-3 no dispone datos horarios para un a˜no completo anteriores a 2005. Los datos de radiaci´on terrestres disponibles en la isla de Gran Canaria pertenecen a a˜nos entre 2000 y 2005 dependiendo de la estaci´on, ver Tabla 2.10. De esta manera, cuando se trabaje con datos provenientes de im´agenes satelitales ´unicamente se puede utilizar el a˜no 2005. Los datos de velocidad y direcci´on del viento y la humedad relativa de la atm´osfera obtenidos del ECMWF tambi´en est´an disponibles ´unicamente para el a˜no 2005 en la isla de Gran Canaria. As´ı, se han realizado dos conjuntos de simulaciones con diferentes modelos estad´ısticos y diferentes variables de entrada. En un primer caso, se utilizar´an dos a˜nos de medida para realizar el entrenamiento y el test y las variables de entrada se reducen a los datos hist´oricos de radiaci´on terrestre y los datos de radiaci´on obtenidos de las predicciones del ECMWF. Los resultados de estas simulaciones se muestran en la Secci´on 5.2. Mientras que al incluir como variables de entrada los datos satelitales y los datos de viento y humedad del ECMWF se trabajar´a ´unicamente con el a˜no 2005. En este caso, el a˜no de medida se dividir´a en un conjunto para el entrenamiento y otro para el test. Los resultados de estas simulaciones se muestran en la Secci´on 5.3. Para evaluar el ajuste de las predicciones de cada m´etodo, se utilizar´an la medidas de error est´andar. En esta tesis se han elegido el Error Cuadr´atico Medio, en ingl´es Root Mean Square Error (RMSE), y el Error Medio Absoluto, en ingl´es Mean Absolute Error (MAE). Ambos errores han sido ampliamente referenciados en la bibliograf´ıa especializada [Lauret15, Willmott05], al igual que sus correspondiente medidas relativas (rRMSE y rMAE). Las medidas relativas se obtienen dividiendo el error absoluto entre la media del conjunto de datos de radiaci´on solar de la zona en cuesti´on. En esta tesis, para decidir la calidad de los diferentes m´etodos de predicci´on se utilizar´an los errores relativos. Adem´as, para comparar la mejora relativa de los diferentes modelos con respecto al modelo m´as sencillo, Persistence Model, se ha calculado un par´ametro llamado %SKILL para cada modelo [Coimbra13]. Un valor positivo y alto de este par´ametro supone que dicho modelo mejora sustancialmente el Persistence Model, mientras que un valor negativo significa que el modelo simulado no mejora los resultados ya obtenidos. RMSEmodelo =v u u t1 N N X i=1 (ˆ Ig,i −Imeasure,i)2(5.1) 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 107 MAEmodelo =1 N N X i=1 (ˆ Ig,i −Imeasure,i)(5.2) SKILL( %) = 1−RMSEmodelo RMSEpersistence x100 (5.3) Donde Nrepresenta el n´umero de datos horarios de la muestra sobre la que se est´a calculando el error, ˆ Ig,i es el valor de radiaci´on solar para la hora icalculado por el modelo de predicci´on, Imeasure,i es el valor de radiaci´on solar obtenido por las estaciones de medida terrestres y RMSEpersistence es el error obtenido por las predicciones realizadas con el modelo Persistence para dicha estaci´on. Las predicciones en los modelos estad´ısticos se realizar´an con el ´ındice de cielo despejado, ˆ (k∗). Cuando se hayan realizado los procesos de entrenamiento de cada uno de los modelos, el c´alculo y discusi´on de los errores se llevar´a a cabo en t´erminos de irradiancias horarias en W/m2, seg´un ˆ Ig=ˆ (k∗)·Ics. 5.2. Predicci´on de Radiaci´on Solar a partir de datos hist´oricos terrestres y del ECMWF En esta Secci´on se describir´an los c´alculos de predicci´on de la radiaci´on solar en la isla de Gran Canaria utilizando datos hist´oricos y datos de radiaci´on solar obtenidos por el modelo num´erico ECMWF. En la Secci´on 2.3 se describieron las estaciones de medida utilizadas en esta tesis. Se pueden observar los a˜nos de medida que se utilizar´an en cada caso, Tabla 2.10, y el porcentaje de datos mensuales disponibles. Todos los modelos estad´ısticos descritos en esta Secci´on se simular´an con los datos disponibles de cada una de las estaciones de la isla: C0-Pozo Izquierdo C1-Las Palmas C2-La Aldea C4-Maspalomas C5-Sta. Br´ıgida C6-Mog´an En un primer momento se estudiar´an las predicciones ´unicamente utilizando la base de datos hist´oricos de cada una de las estaciones. Como ya se ha comentado, para estudiar la complejidad del modelo y obtener el modelo ´optimo para realizar las predicciones se utilizar´a el a˜no definido como de entrenamiento. Mientras que la comprobaci´on de los errores obtenidos por dicho modelo se obtendr´an con el a˜no definido como de test. En esta tesis se estudiar´an los siguientes modelos estad´ısticos: 108 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON Modelo Persistence Modelo Smart-Persistence Modelo Climatol´ogico Modelo lineal Autorregresivo AR Modelo lineal Autorregresivo de Medias M´oviles ARMA Modelo de Redes Neuronales Artificiales RNAs 5.2.1. Modelos lineales de predicci´on Los modelos estad´ısticos lineales han sido utilizados en el modelado de series temporales durante a˜nos. En esta tesis se estudiar´a la capacidad de predicci´on de la radiaci´on solar de los dos siguientes modelos lineales: Modelos Autorregresivos (AR): Se trata de un modelo de regresi´on lineal basado ´unicamente en datos pasados de la serie temporal que se estudia. El orden de modelo viene definido por la variable p, que representa el n´umero de par´ametros de los que consta la funci´on de transferencia. Modelos Autorregresivos de medias m´oviles (ARMA): Se trata de un modelo formado por dos partes, una basada en el modelo autorregresivo (AR) y otra basada en el modelo de medias m´oviles (MA). Las complejidad del modelo vendr´a definida por los par´ametros pyq, que representan el orden de la parte (AR) y (MA) respectivamente. 5.2.1.1. Estudio de la complejidad de los modelos lineales A la vista de lo anterior, parece obvio que estudiar la complejidad del modelo lineal que mejor se ajusta a nuestros datos es un primer paso muy importante. La variable p representar´a en ambos casos el orden del modelo AR, definiendo el n´umero de par´ametros del mismo. Por tanto, la definici´on de esta variable se˜nalar´a el n´umero de datos pasados de la serie temporal necesarios para establecer una relaci´on lineal con el valor a predecir. Como se ha explicado, la variable utilizada para trabajar con los modelos de cielo despejado es el ´ındice de cielo despejado k∗. As´ı, los datos pasados de la serie temporal se refieren a los valores del ´ındice de cielo despejado para las horas anteriores al momento en el que se quiere realizar la predicci´on t. Modelo autorregresivo AR El modelo autorregresivo para predecir un ´ındice de cielo despejado en el horizonte temporal hrespecto al momento actual queda definido con la ecuaci´on 4.6. b k∗(t+h) = p−1 X i=0 [Φi+1k∗(t−i)] + t+h 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 109 Los valores k∗(t−i) son los datos pasados de la serie temporal, mientras prepresenta el orden de AR. Por otro lado, {Φi}i=1,2,...,p muestra los par´ametros de autorregresi´on obtenidos a partir de los datos de la muestra durante el entrenamiento. Para tomar una decisi´on del modelo ´optimo a utilizar, es decir el orden pdel mismo, se utilizar´an las gr´aficas de SPACF para identificar el m´aximo orden de los modelos AR. En cada estaci´on se realizar´a la muestra de Autocorrelaci´on Parcial SPACF para un conjunto de retardos de tiempo. La autocorrelaci´on parcial PACF establece la correlaci´on entre dos instantes de tiempo de la serie con un retardo determinado. De esta manera, para identificar los valores de SPACF significativos y establecer el retardo a partir del cual la relaci´on entre las variables es pr´acticamente nula se ha considerado que el 95 % de los valores se encuentran dentro del intervalo (−1.96/√N, +1.96/√N). Utilizando este l´ımite se establecer´a para cada estaci´on el valor m´aximo de retardo significativo, que representar´a el n´umero de datos pasados necesarios para realizar una predicci´on ´optima y, por lo tanto, el orden del modelo p. Figura 5.1: SPACF de la serie temporal del ´ındice de claridad horario para la estaci´on de C0-Pozo Izquierdo. En rojo se representan los intervalos de elecci´on Como se observa para las estaciones de C0 y C5, el primer valor tiene una importancia muy superior a los siguientes, cuyo valor de autocorrelaci´on parcial va decayendo hasta superar el l´ımite marcado, Figuras 5.1 y 5.2. En la estaci´on de C0 el valor m´aximo se encuentra entorno al 10, mientras que en C5 se puede entender que el valor m´aximo est´a en 14. Para realizar unas simulaciones homog´eneas en todas las estaciones se ha elegido como valor m´aximo del orden del modelo p= 14. Estos valores de pno ser´an los valores finalmente utilizados para ajustar la serie, sino que se realizar´an los c´alculos del par´ametro BIC, ecuaci´on 4.15, para todos los 110 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON Figura 5.2: SPACF de la serie temporal del ´ındice de claridad horario para la estaci´on de C5-Sta. Br´ıgida. En rojo se representan los intervalos de elecci´on ´ordenes, p= 1...14, hasta el m´aximo seleccionado para todos los horizontes temporales de predicci´on, h= 1...6h. El par´ametro BIC (Bayesian Information Criterion) tiene en cuenta la complejidad del modelo a la hora de decidir el ´optimo. El modelo ´optimo a utilizar ser´a el que muestre un valor BIC menor. BIC =ln(ˆσ2) + mk ln(N) N Adem´as del estudio del Bayesian Information Criterion BIC, se estudiaron los valores del error cuadr´atico medio relativo de modelo de ajuste %rRMSE. En muchas ocasiones, el modelo ´optimo se˜nalado por el BIC no ofrec´ıa una mejora sustancial del error con respecto a un modelo m´as sencillo, en el que se utilizaban un menor n´umero de valores pasados de la serie temporal. Se puede observar en las Figuras 5.3 y 5.4, que el valor del par´ametro BIC desciende significativamente conforme se aumenta el orden del modelo AR para todos los horizontes temporales. Este descenso se detiene entorno al orden p= 10, por lo que el modelo ´optimo seg´un este criterio se encontrar´a entorno a este orden. Dependiendo de la estaci´on y del horizonte temporal se calcul´o el n´umero de orden que proporcionaba el BIC m´as bajo. Por ejemplo, para estaci´on C0-Pozo Izquierdo se obtuvieron los ´ordenes p={9,9,9,10,10,10} para los horizontes temporales de h={1,2,3,4,5,6}respectivamente. Por otro lado, en las Figuras 5.5 y 5.6, se muestran los errores obtenidos en la predicciones realizadas por los distintos modelos AR. En este caso, se observa en las estaciones del sur, como C0, que la diferencia de error comienza a ser significativa a partir 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 111 Figura 5.3: Valores del par´ametro BIC obtenidos para la estaci´on de C0-Pozo Izquierdo. Se muestran la simulaciones de un modelo AR con ´ordenes desde 1 hasta 14 para todos los horizontes temporales de predicci´on h= 1...6h Figura 5.4: Valores del par´ametro BIC obtenidos para la estaci´on de C5-Sta. Brigida. Se muestran la simulaciones de un modelo AR con ´ordenes desde 1 hasta 14 para todos los horizontes temporales de predicci´on h= 1...6h 118 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON (a) (b) Figura 5.11: Error cuadr´atico medio relativo %rRMSE y BIC para la estaci´on de C5-Sta. Br´ıgida. Se muestran la simulaciones de un modelo ARMA para todas las combinaciones posibles de pyq. En cada parte se observa a) BIC 1h, b) %rRMSE 1h 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 119 (a) (b) Figura 5.12: Error cuadr´atico medio relativo %rRMSE y BIC para la estaci´on de C5-Sta. Br´ıgida. Se muestran la simulaciones de un modelo ARMA para todas las combinaciones posibles de pyq. En cada parte se observa a) BIC 6h, b) %rRMSE 6h 120 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON Estaciones 1 h 2 h 3 h 4 h 5 h 6 h C0 9 9 9 10 10 10 C1 11 11 12 14 14 14 C2 8 8 9 11 11 11 C4 10 10 11 11 11 11 C5 9 9 10 10 10 10 C6 7 8 9 10 11 11 Tabla 5.1: ´ Ordenes del modelo AR elegido para cada estaci´on y horizonte temporal h Otro dato importante y que podemos ir observando desde este momento, es que en las estaciones C0, C2 y C6 los resultados de predicci´on son claramente mejores que en el resto de las estaciones. Estas estaciones se encuentran en el sur de la isla y poseen, como ya se describi´o en el Cap´ıtulo 2, un clima m´as estable y con un mayor n´umero de d´ıas despejados. A la vista de estos resultados, a partir de este momento se utilizar´a el modelo lineal ARMA para comparar con el resto de modelos de predicci´on. Los resultados obtenidos son ligeramente mejores que con el modelo AR utilizando un n´umero menor de datos pasados de radiaci´on solar. 5.2.2. Modelos de predicci´on basados en Redes Neuronales Artificiales La Redes Neuronales Artificiales (RNAs) [Bishop95] son modelos matem´aticos de predicci´on que encuentran una relaci´on ´optima entre un conjunto de datos de entrada y una salida deseada. La relaci´on se establece mediante un proceso de entrenamiento. Las RNA est´an compuestas por unidades, llamadas neuronas, que reciben una entrada de otra neurona o de una fuente externa. Todas las conexiones entre neuronas se representan con un peso asociado que se ir´a modificando durante el entrenamiento. Cada neurona recibe una entrada afectada por esos pesos y realiza la suma de todas las entradas para producir una nueva salida hacia la neurona siguiente. La suma se ve afectada por una funci´on de activaci´on o funci´on de transferencia para limitar la salida de cada neurona. En esta tesis, todas las neuronas se calcular´an utilizando una funci´on tangente hiperb´olica. El Perceptr´on Multicapa (MLP) es una de las arquitecturas de RNAs m´as utilizada por la bibliograf´ıa. Las RNAs utilizadas durante este trabajo tienen la estructura de un MLP. En concreto, se optado por una capa de entrada, una capa intermedia o capa oculta y una capa de salida, sin conexiones de retroalimentaci´on ni conexiones laterales entre las citadas capas. La capa oculta est´a formada por varias neuronas con una funci´on de activaci´on tangente hiperb´olica. La capa de entrada representa los datos de la serie temporal de ´ındices de cielo despejado, siendo Tel n´umero de datos pasados elegidos para predecir la salida, mientras que la capa de salida representa el ´ındice de cielo despejado predicho para el horizonte temporal h,ˆ k∗(t+h). 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 121 (a) C0 (b) C1 (c) C2 (d) C4 (e) C5 (f) C6 Figura 5.13: Error cuadr´atico medio relativo %rRMSE para todas las estaciones de medida y todos los horizontes temporales de predicci´on h= 1...6h. El modelo AR utilizado corresponde al ´optimo se˜nalado por el BIC, mientras que el modelo ARMA corresponde a un orden p= 2 y q= 1 122 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON ˆ k∗(t+h) = H X j=1 ω2 sjfj[ T−1 X i=0 (ω1 jik∗(t−i) + ω1 0] + ω2 0(5.4) Los pesos de la RNA se inician aleatoriamente y son luego optimizados durante el per´ıodo de entrenamiento. Este proceso se realiza minimizando una funci´on de error (funci´on de coste) mediante el algoritmo de retroalimentaci´on (backpropagation). La funci´on de error m´as utilizada es funci´on del error cuadr´atico medio entre la salida obtenida por la red, ˆ k∗(t+h), y la salida deseada k∗(t+h) (en nuestro caso el´ındice de cielo despejado medido en las estaciones), target set. La optimizaci´on de la RNA se realiza utilizando el conjunto de datos de entrenamiento. La precisi´on de las RNAs para aproximar las series temporales depende en gran medida de la estructura de la red. En muchas ocasiones se consigue realizar una buena aproximaci´on de los datos de entrenamiento, pero peores resultados cuando se utilizan datos nuevos (este problema se conoce como sobreajuste). Para evitar este problema se debe optimizar la estructura de la red utilizando tambi´en el conjunto de datos de test. En la bibliograf´ıa existen diversas t´ecnicas de regularizaci´on para estudiar la complejidad del modelo [Bishop95, Lauret08, Lauret06b]. En esta tesis se ha optado por las t´ecnicas Bayesianas de regularizaci´on de RNAs [MacKay03]. La aproximaci´on Bayesiana considera una funci´on de densidad de probabilidad en el espacio de los pesos. Los valores ´optimos del conjunto de pesos corresponden a aquellos que hacen m´aximo el valor de la funci´on de probabilidad. Las t´ecnicas Bayesianas introducen dos nuevos hiperpar´ametros αyβa la funci´on de error para controlar la complejidad del modelo. S(ω) = β 2 N X i=1 [ˆ k∗ i(t+h)−k∗ i(t+h)] + α 2 m X k=1 ω2 k=β 2ED+α 2Eω(5.5) Para realizar todos los c´alculos de entrenamiento y optimizaci´on de la complejidad del modelo se ha utilizado la herramienta NETLAB, desarrollada en el programa MATLAB para desarrollo de RNAs [Nabney02]. 5.2.2.1. Estudio de la complejidad de las RNAs Como se ha mencionado, uno de los aspectos m´as importantes a la hora de realizar una predicci´on es decidir la complejidad de la RNA. En el caso del MLP se controlar´a el n´umero de entradas y el n´umero de neuronas de la capa oculta. El n´umero de entradas representa, en el caso de las predicciones de radiaci´on solar, el n´umero de datos pasados de radiaci´on solar terrestre. Las t´ecnicas Bayesianas ofrecen diferentes posibilidades para estudiar la complejidad del modelo. Para valorar la importancia de cada entrada, se dividen los pesos en un grupo por cada entrada, un grupo para los pesos de la segunda capa y un grupo para los umbrales de la primera y segunda capa por separado. Cada grupo es controlado por un hiperpar´ametro independiente αg, as´ı que todos los pesos relacionados al mismo grupo est´an afectados por el mismo hiperpar´ametro. Esta t´ecnica, conocida como Automatic 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 123 Relevance Determination (ARD), asigna un coeficiente de regularizaci´on diferente a cada entrada y por lo tanto determina aqu´ellas que son m´as relevantes durante el proceso de entrenamiento. Al final del entrenamiento, los pesos con un gran valor de αgson cercanos a cero. Las entradas, cuyos pesos correspondientes son cercanos a cero, se consideran no relevantes para la RNA y podr´ıan ser eliminadas. Adem´as, las t´ecnicas Bayesianas tambi´en nos permiten elegir el n´umero de neuronas ´optimo. En este caso se calcula la probabilidad del modelo (evidence of the model) y se utiliza esta medida para elegir el mejor de los modelos propuestos. La bibliograf´ıa ofrece diferentes f´ormulas para calcular el log of evidence del modelo [Bishop95, MacKay03, Penny99]. En esta tesis se utiliz´o finalmente la ecuaci´on 4.64. logP(Mi|D) = −αMP EMP ω−βMP EMP D−1 2log|A|+m 2logαMP + +N 2logβMP +1 2log 2 γ+1 2log 2 N−γ Primero, se estudi´o la informaci´on ofrecida por la t´ecnica ARD para elegir el n´umero de entradas ´optimo para la RNA. El estudio se realiz´o con el n´umero m´aximo de 6 entradas, igual al horizonte temporal m´aximo. Como se explic´o, las t´ecnicas Bayesianas asignan durante el entrenamiento, valores de cercanos a cero a los pesos que se consideran no relevantes para la red. En las Figuras 5.14 y 5.15 se pueden observar los resultados de la t´ecnica ARD para todos los horizontes temporales en las estaciones C0 y C5. La informaci´on ofrecida no es constante para todos los horizontes temporales ni para los distintos n´umeros de neuronas ocultas utilizadas. En general, se puede considerar que las seis entradas elegidas aportan un grado de relevancia en algunos de los casos estudiados. En cualquier caso, se opt´o por podar aquellas entradas que la t´ecnica ARD consideraba no relevantes en cada simulaci´on y estudiar el error obtenido con el conjunto del test. Los resultados obtenidos indicaban que siempre era mejor utilizar las seis entradas y, durante el entrenamiento, permitir a las t´ecnicas Bayesianas establecer la importancia de los pesos. Es decir, finalmente se utilizar´an seis entradas de datos pasados del ´ındice de cielo despejado y la RNA establecer´a los valores de los pesos seg´un la relevancia de las entradas mediante las t´ecnicas Bayesianas. Una vez qued´o establecido el n´umero de datos pasados que se utilizar´an como entradas, se estudi´o el n´umero de neuronas ocultas. Para cada estaci´on y cada horizonte temporal de predicci´on se simularon varios modelos con distinto n´umero de neuronas ocultas y la decisi´on se tom´o en base al valor del log of evidence de cada modelo simulado. La mayor´ıa de los resultados obtenidos conducen a RNAs con n´umero de neuronas ocultas reducido. En la Figura 5.16 se observan los resultados para la estaci´on C0 y los horizontes temporales h= 1 y h= 3. En el primer caso se puede observar que el m´aximo se sit´ua en 4 neuronas ocultas, mientras que en el segundo el m´aximo se sit´ua en 3 neuronas, aunque con una diferencia muy peque˜na con respecto a utilizar 2 ´o 4 neuronas. Como ya se indic´o en la Secci´on 4.5.3, la informaci´on ofrecida por el log of evidence se complementa con un estudio del %rRMSE de las predicciones con respecto al conjunto de datos del test. As´ı, a partir del resultado ´optimo observado en la Figura 5.16, se eligen diversos modelos con un n´umero de 124 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON (a) 1h (b) 2h (c) 3h (d) 4h (e) 5h (f) 6h Figura 5.14: Resultado del ARD utilizando el n´umero ´optimo de neuronas para cada horizonte temporal en la estaci´on de C0-Pozo Izquierdo. Cada entrada representa los datos terrestres pasados de la serie temporal 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 125 (a) 1h (b) 2h (c) 3h (d) 4h (e) 5h (f) 6h Figura 5.15: Resultado del ARD utilizando el n´umero ´optimo de neuronas para cada horizonte temporal en la estaci´on de C5-Sta. Br´ıgida. Cada entrada representa los datos terrestres pasados de la serie temporal 126 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON neuronas ocultas %rrmse entrenamiento %rmae entrenamiento %rrmse test %rmae test 2 14.96 9.74 17.37 11.91 3 14.86 9.65 17.45 12.00 4 14.74 9.47 17.39 11.62 5 14.56 9.19 17.36 11.40 6 14.55 9.23 17.58 11.73 7 14.60 9.35 17.67 11.75 10 14.42 9.09 17.95 11.73 15 14.03 8.71 18.75 11.94 20 13.60 8.46 20.07 12.15 Tabla 5.2: Errores del conjunto de datos de radiaci´on terrestres de entrenamiento y de test utilizando distintos n´umeros de neuronas ocultas en la estaci´on de C0-Pozo Izquierdo para el horizonte temporal 1 h. neuronas alrededor del mismo. Con estos modelos se realiza el estudio del error %rRMSE para decidir el n´umero final de neuronas ocultas. Para la estaci´on C0 y los horizontes temporales h= 1 y h= 3 se muestran los resultados obtenidos en las Tablas 5.2 y 5.3. Aunque los errores %rRMSE y %rMAE observados con el conjunto de datos del test no difieren significativamente, se eligi´o el n´umero de neuronas cuya predicci´on ofrec´ıa un menor error en cada caso. De esta manera, para un horizonte temporal h= 1 se obtuvo una RNA con 5 neuronas y para h= 3 se obtuvo una RNA con 4 neuronas. 5.2.3. Resultados obtenidos y comparaci´on entre los modelos En este apartado se estudian los resultados obtenidos en las predicciones para los distintos horizontes temporales en cada estaci´on de medida. El objetivo es encontrar una relaci´on entre las variables de entrada, datos hist´oricos del ´ındice de cielo despejado, y el valor futuro a predecir. La funci´on general que relaciona los par´ametros de entrada y salida es la ecuaci´on (4.2), b k∗(t+h) = Fk∗ g(t), . . . , k∗ g(t−i), k∗ e1(t), . . . , k∗ e1(t−j), . . . , k∗ en(t), . . . , k∗ en(t−j) donde b k∗(t+h) es el ´ındice de cielo despejado calculado para un horizonte de predicci´on h, en nuestro caso h= 1,2, ..., 6horas,k∗ g(t−i) es el ´ındice de cielo despejado para las i horas pasadas obtenido a partir de los datos terrestres medidos en la estaci´on y k∗ en(t−j) corresponde ´ındice de cielo despejado para un las jhoras pasadas obtenido a partir del dato ex´ogeno en. Si los modelos na¨ıve o modelos de referencia predicen con mayor precisi´on que los modelos propuestos no resulta interesante desarrollar metodolog´ıas sofisticadas. De esta manera, en este trabajo, se propone comparar las predicciones obtenidas con las RNAs 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 127 (a) 1h (b) 3h Figura 5.16: Resultado del Log of evidence para decidir el n´umero de neuronas ´optimo en la capa oculta para la estaci´on de C0-Pozo Izquierdo 134 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON estaci´on m´etodos 1 h 2 h 3 h 4 h 5 h 6 h C0 Persistence - - - - - - Smart Persistence 0.0000 -0.52 2.67 9.50 16.85 21.16 Climatol´ogico -58.79 -11.5864 7.35 18.34 23.67 25.72 ARMA 4.10 11.55 18.03 23.48 26.65 27.87 NN 2.21 10.77 18.45 24.10 28.03 28.44 NN+ECMWF 4.22 11.1 18.61 26.317 30.81 32.20 media IGH = 531.74 Wm−2 C1 Persistence - - - - - - Smart Persistence 0.00 -1.47 1.75 7.07 14.10 19.46 Climatol´ogico -44.70 -3.65 11.39 19.00 23.50 24.95 ARMA 6.21 13.16 18.47 22.37 25.25 26.34 NN 6.88 13.34 18.66 23.49 27.15 28.61 NN+ECMWF 7.75 16.52 23.64 28.45 31.49 32.77 media IGH = 441.98 Wm−2 C2 Persistence - - - - - - Smart Persistence 0.00 -1.13 -0.03 5.69 11.94 16.23 Climatol´ogico -63.50 -17.96 0.74 12.61 19.23 21.96 ARMA 3.35 9.19 14.48 19.63 23.17 24.58 NN 2.06 7.56 14.57 19.90 21.74 26.70 NN+ECMWF 5.01 10.79 16.22 24.74 28.83 30.86 media IGH = 551.93 Wm−2 C4 Persistence - - - - - - Smart Persistence 0.00 -4.87 5.34 14.65 24.38 31.34 Climatol´ogico -23.93 20.77 36.04 40.62 41.86 40.30 ARMA 8.36 20.55 29.66 34.46 36.70 36.25 NN 12.32 28.25 39.51 43.27 44.51 42.48 NN+ECMWF 17.93 32.74 39.66 45.57 45.89 44.85 media IGH = 441.49 Wm−2 C5 Persistence - - - - - - Smart Persistence 0.00 -2.87 0.25 7.70 15.42 20.31 Climatol´ogico -74.64 -21.35 0.15 11.90 17.72 19.60 ARMA 3.83 10.16 16.09 21.18 24.51 25.82 NN 2.24 8.95 16.47 21.94 26.54 28.12 NN+ECMWF 7.30 15.41 22.94 28.54 32.47 33.24 media IGH = 440.93 Wm−2 C6 Persistence - - - - - - Smart Persistence 0.00 3.00 10.44 17.87 22.54 25.50 Climatol´ogico -37.77 2.28 16.29 23.20 24.78 24.89 ARMA 8.01 18.73 25.78 30.29 31.36 31.25 NN 12.15 21.59 27.31 31.79 33.25 32.21 NN+ECMWF 18.90 30.12 35.78 39.81 40.04 38.77 media IGH = 491.36 Wm−2 Tabla 5.6: Comparaci´on con el modelo Persistence mediante el par´ametro %SKILL para todos los horizontes temporales h= 1...6 en todas las estaciones de medida. 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 135 el horizonte temporal. Con el Smart-Persistence se consigue mejorar los resultados al aumentar el horizonte de predicci´on porque cada vez utiliza mayor cantidad de datos pasados para estimar el valor futuro de radiaci´on. Este modelo tiende a igualar el modelo Climatol´ogico de predicci´on, como es obvio por su formulaci´on. Los modelos ARMA y los basados en las redes neuronales en las primeras horas muestran una ligera mejora con respecto a los modelos na¨ıve, sin embargo, conforme aumentan las horas esta mejor´ıa muestra valores relevantes. Adem´as, tanto el ARMA como las NN tienden a acercarse al modelo Climatol´ogico conforme aumenta el horizonte de predicci´on, aunque siempre mostrando una mejor´ıa. Al a˜nadir los datos del ECMWF a las redes neuronales se consigue establecer una mejor´ıa respecto a todos los modelos en todo el horizonte de predicci´on, incluyendo una mejor´ıa respecto al modelo Climatol´ogico para las ´ultimas horas. En la Tabla 5.6 se puede estudiar la mejor´ıa de cada uno de los modelos con respecto al Persistence ( %SKILL). En ella se observa que los modelos mejoran las predicciones respecto al Persistence (cuanto mayor sea el par´ametro %SKILL) conforme aumenta el horizonte temporal de predicci´on, lo cu´al quiere decir que la idoneidad de estos m´etodos respecto a los m´etodos sencillos es mayor conforme nos alejamos en el tiempo. A partir de la tercera hora de predicci´on ya no se observa ning´un valor negativo, con lo que todos los modelos muestran mejores resultados que el Persistence. Se concluye que tanto en las estaciones del sur como en el norte, los mejores resultados de predicci´on se consiguen con el modelo NN+ECMWF para todas las estaciones y todos los horizontes temporales. Por otro lado, aunque los resultados en cuesti´on de errores de predicci´on absolutos y relativos son mejores en las estaciones del sur, la mejor´ıa respecto al resto de modelos es independiente de la estaci´on de la isla. En general, en todas las estaciones de la isla el modelo NN+ECMWF ofrece una mejor´ıa entorno al 30 %SKILL. En el an´alisis de los datos de radiaci´on se explic´o la diferencia entre las estaciones de la isla de Gran Canaria seg´un las estaciones del a˜no. Adem´as de los resultados de predicci´on anuales se decidi´o estudiar el ajuste de los modelos seg´un la ´epoca del a˜no. En las Tablas 5.7 y 5.8 se muestran los resultados de predicci´on para un horizonte temporal de h= 1 yh= 6 respectivamente. Los errores se presentan en t´erminos del %rRMSE para el a˜no completo (Anual) y para cada una de las estaciones del a˜no (Invierno, Primavera, Verano y Oto˜no). En general se puede observar que en las estaciones del sur del isla, C0 y C2, existe una diferencia en los resultados de la predicci´on entre los meses de verano y los meses de invierno. As´ı por ejemplo, utilizando el mejor m´etodo de predicci´on (NN+ECMWF), los valores de %rRMSE obtenidos en C0 para el verano oscilan entre 10-19 % para todo el horizonte temporal de predicci´on, mientras que en invierno esta oscilaci´on se produce entre 25-38 %. En las estaciones del norte de la isla la diferencia entre el verano y el invierno no es tan acusada, debido a la formaci´on de nubosidad durante los meses de verano por el efecto continuo de los vientos Alisios. En la estaci´on de C1 por ejemplo, los errores para el verano oscilan entre 22-33 % para todo el horizonte temporal y en invierno oscilan entre 27-37 %. En los meses de verano, para la primera hora de predicci´on, el m´etodo de redes neuronales con datos del ECMWF consigue unos errores incluso menores del 10 % en algunas estaciones del sur, mientras que en C1 el error no baja del 20 %. Entre las estaciones del norte de la isla, C1-Las Palmas y C5-Sta. Br´ıgida, los errores 136 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON (a) C0 (b) C1 (c) C2 (d) C4 (e) C5 (f) C6 Figura 5.17: Evoluci´on del %rRMSE en funci´on del tiempo horizonte h= 1...6 para cada modelo de predicci´on 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 137 (a) C0 (b) C1 (c) C2 (d) C4 (e) C5 (f) C6 Figura 5.18: Evoluci´on del %rMAE en funci´on del tiempo horizonte h= 1...6 para cada modelo de predicci´on 138 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON estaci´on m´etodos Anual Invierno Primavera Verano Oto˜no C0 Persistence 17.75 25.37 17.60 10.93 22.44 Smart Persistence 17.75 25.37 17.60 10.93 22.44 Climatol´ogico 28.18 45.39 24.34 20.52 33.20 ARMA 17.02 24.47 16.67 10.97 21.12 NN 17.36 26.42 16.77 10.51 21.53 NN+ECMWF 17.00 25.30 16.62 10.09 21.41 C1 Persistence 26.00 29.20 23.04 24.20 29.61 Smart Persistence 26.00 29.20 23.04 24.20 29.61 Climatol´ogico 37.62 42.11 36.74 35.98 37.39 ARMA 24.39 27.36 21.76 22.76 27.23 NN 24.22 27.88 21.53 22.43 27.06 NN+ECMWF 23.99 26.95 21.26 22.42 27.10 C2 Persistence 15.35 19.45 18.80 8.95 16.16 Smart Persistence 15.35 19.45 18.80 8.95 16.16 Climatol´ogico 25.09 30.76 26.66 20.47 26.92 ARMA 14.84 18.60 17.76 9.28 15.85 NN 15.02 18.70 18.04 9.49 15.64 NN+ECMWF 14.56 17.76 17.48 9.30 15.31 C4 Persistence 27.93 37.81 18.61 22.87 30.81 Smart Persistence 27.93 37.81 18.61 22.87 30.81 Climatol´ogico 34.60 44.56 26.76 28.00 39.03 ARMA 25.60 34.48 17.85 20.71 28.17 NN 24.50 30.97 18.88 20.88 26.93 NN+ECMWF 22.93 30.64 - 17.93 26.37 C5 Persistence 24.10 25.68 30.84 17.53 25.28 Smart Persistence 24.10 25.68 30.84 17.53 25.28 Climatol´ogico 42.08 46.15 46.36 37.45 40.39 ARMA 23.18 25.20 29.54 16.52 24.55 NN 23.55 25.37 29.65 18.02 23.72 NN+ECMWF 22.33 24.05 28.97 15.84 23.50 C6 Persistence 23.06 38.61 21.21 13.77 28.48 Smart Persistence 23.06 38.61 21.21 13.77 28.48 Climatol´ogico 31.76 50.04 28.76 23.55 36.13 ARMA 21.21 35.20 19.55 12.60 26.23 NN 20.25 31.88 21.19 11.62 24.25 NN+ECMWF 18.70 30.96 18.94 10.10 22.66 Tabla 5.7: Errores %rRMSE anuales y por trimestres para el horizonte temporal h= 1 en todas las estaciones de medida. 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 139 estaci´on m´etodos Anual Invierno Primavera Verano Oto˜no C0 Persistence 37.96 44.59 34.69 32.35 46.90 Smart Persistence 29.93 42.03 29.15 24.26 32.07 Climatol´ogico 28.20 45.59 24.34 20.47 33.37 ARMA 27.38 39.63 25.24 20.57 32.40 NN 27.16 42.83 24.44 19.65 31.40 NN+ECMWF 25.73 38.17 23.84 19.06 29.78 C1 Persistence 50.12 53.22 45.80 50.19 51.49 Smart Persistence 40.37 45.01 39.17 39.41 38.70 Climatol´ogico 37.63 42.24 36.78 35.80 37.52 ARMA 36.92 41.06 34.60 35.31 37.25 NN 35.79 40.23 34.72 33.88 35.67 NN+ECMWF 33.71 36.57 31.86 33.02 34.23 C2 Persistence 32.07 36.44 38.52 23.27 31.42 Smart Persistence 26.87 34.86 29.61 20.51 27.50 Climatol´ogico 25.00 30.10 26.74 20.51 26.40 ARMA 24.19 29.88 26.84 18.84 24.78 NN 23.47 28.27 26.41 18.34 23.52 NN+ECMWF 22.14 26.63 25.05 16.25 23.43 C4 Persistence 57.90 57.93 36.60 59.71 54.12 Smart Persistence 39.76 47.74 31.41 38.36 39.71 Climatol´ogico 34.58 44.83 26.81 27.91 39.25 ARMA 36.91 44.93 26.32 32.72 38.66 NN 33.38 47.62 27.90 26.03 37.83 NN+ECMWF 32.01 40.27 26.98 35.56 C5 Persistence 52.38 53.38 60.90 45.28 52.11 Smart Persistence 41.74 45.37 49.95 33.39 43.088 Climatol´ogico 42.09 46.23 46.50 37.52 40.44 ARMA 38.86 44.31 46.65 29.80 40.21 NN 37.62 42.10 43.71 30.64 38.78 NN+ECMWF 34.94 37.07 41.78 27.67 37.58 C6 Persistence 42.30 58.17 40.06 33.08 49.22 Smart Persistence 31.51 48.45 29.49 22.22 37.80 Climatol´ogico 31.76 50.05 28.77 23.49 36.32 ARMA 29.08 46.59 26.86 18.28 35.73 NN 28.66 45.23 27.07 19.07 34.72 NN+ECMWF 25.89 43.26 24.16 15.82 31.57 Tabla 5.8: Errores %rRMSE anuales y por trimestres para el horizonte temporal h= 6 en todas las estaciones de medida. 140 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON en verano tampoco siguen una misma tendencia, encontr´andose unos errores claramente superiores en la estaci´on C1 para todos los horizontes de predicci´on. Esto se debe a que la estaci´on C5 se encuentra a mayor altitud sobre el nivel del mar (C1 est´a en la costa) y no se ve afectada en gran medida por la formaci´on de nubes t´ıpica de esta ´epoca. Por otro lado, entre las estaciones del sur, C2 muestra unos errores similares en las estaciones de invierno, primavera y oto˜no, con errores entre el 15-18 % para h= 1 y 23-26 % para h= 6, mientras que en la estaci´on de C0 en invierno y oto˜no los valores del %rRMSE son superiores, superando incluso el 30 % para h= 6. Aunque ambas estaciones muestran un clima m´as estable que en el norte y con mayor n´umero de d´ıas despejados, C2-La Aldea en el suroeste de la isla muestra un clima m´as estable debido a que se encuentra completamente opuesta a la direcci´on de los vientos predominantes de la isla. La diferencia del mejor modelo de predicci´on (NN+ECMWF) respecto al modelo Persistence aumenta con el horizonte de predicci´on en las estaciones del sur. Para las primeras horas de predicci´on pr´acticamente no se consigue una mejora relevante, incluso en la estaci´on C2 el modelo Persistence consigue resultados inferiores al 9 %. En el horizonte temporal h= 6 para dichas estaciones se consiguen mejoras que oscilan entre 8-15 %. Por otro lado, en las estaciones del norte, las mejoras se consiguen desde la primera hora de predicci´on, donde se observa entorno a un 2 %, hasta la ´ultima h= 6, donde oscila entre un 15-20 %. Al a˜nadir datos ex´ogenos del ECMWF (NN+ECMWF) las predicciones alcanzan un %rRMSE menor para todas las estaciones de a˜no independientemente de la ubicaci´on. Los resultados obtenidos, para un horizonte temporal h= 1, muestran una diferencia de error que no supera el 1 % en comparaci´on con los modelos ARMA y NN. En cambio, para h= 6, en las estaciones del sur el m´etodo NN+ECMWF logra unas mejoras entorno al 1-2 % dependiendo de la estaci´on del a˜no, y en las estaciones del norte estas mejoras alcanzan hasta un 4-5 % en los meses de invierno. Entre C1 y C5, en el norte ambas, se vuelve a observar una diferencia en verano, ya que en la primera la mejora se sit´ua en torno al 1 % debido a la nubosidad. 5.2.3.2. Resultados seg´un el tipo de d´ıa En la Secci´on 2.3.4 se realiz´o el an´alisis de los datos de radiaci´on solar en cada estaci´on de medida en funci´on de la media y la variabilidad diaria del ´ındice de cielo despejado. Los d´ıas se dividieron en nueve tipos de d´ıas seg´un la Tabla 5.9, donde por ejemplo un d´ıa CI representa un d´ıa completamente despajado con una alta radiaci´on y un d´ıa tipo AIII ser´a un d´ıa con alta variabilidad y una media de radiaci´on baja. En las Tablas 5.10-5.13 se muestran los resultados obtenidos para todos los tipos de d´ıas con cada uno de los modelos de predicci´on en las estaciones de C0, C1, C2 y C5. Se debe tener en cuenta la cantidad de datos disponible para cada tipo de d´ıa y estaci´on de medida. Los tipos de d´ıas cuyo n´umero en el conjunto de datos del test no superen el 5 % no se mostrar´an en las tablas ya que los resultados no se consideran concluyentes. As´ı, por ejemplo en las estaciones del sur de la isla no se muestran en las tablas d´ıas tipo A debido al escaso n´umero de los mismos. Los resultados para la predicci´on de d´ıas tipo I presentan los mejores valores, ya que son 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 141 Tipo de datos Condici´on Ak∗<0.5 B 0.5< k∗<0.9 Ck∗>0.9 Istd(∆k∗)<0.05 II 0.05 < std(∆k∗)<0.15 III std(∆k∗)>0.15 Tabla 5.9: Condiciones a cumplir por cada uno de los subconjuntos de datos. los d´ıas con menor variabilidad. En general, para todas las estaciones de medida, los d´ıas tipo C, cuya media de radiaci´on es la m´as alta, presentan los mejores resultados, mientras que los peores d´ıas son los tipo A, con medias de radiaci´on muy baja. En cuanto a los diferentes m´etodos de predicci´on, en los d´ıas tipo B y C los modelos de redes neuronales con datos de radiaci´on del ECMWF (NN+ECMWF) han obtenido los mejores resultados para todos los horizontes temporales en todas las estaciones de la isla. La diferencia de los resultados de estos modelos no lineales respecto a los modelos Persistence o ARMA es significativamente menor en los d´ıas tipo A. En concreto, en los d´ıas tipo AII los mejores modelos para las primeras horas de predicci´on son los modelos Persistence y en los d´ıas tipo AIII el modelo ARMA presenta los mejores resultados. En la estaci´on C0-Pozo Izquierdo, Tabla 5.10, el 94 % de los d´ıas del conjunto del test son tipo B y C. En ambos tipos de d´ıas el mejor modelo de predicci´on es el NN+ECMWF para todos los horizontes temporales de predicci´on. Se puede observar claramente que las mejoras m´as importantes se establecen para los d´ıas tipo C. En estos d´ıas, desde la segunda hora de predicci´on se consiguen mejoras del orden del 2 % respecto al modelo ARMA. Conforme se aumenta el horizonte de predicci´on se alcanzan mejoras m´as relevantes, llegando al 7 % respecto al modelo ARMA y al 2-3 % respecto a las redes neuronales sin datos ex´ogenos (NN). La diferencia entre el modelo NN y NN+ECMWF en los d´ıas BIII es pr´acticamente nula para todas las horas, mientras que para los d´ıas BII a partir de la tercera hora aumenta paulatinamente hasta el 2 %. En ambos tipos de d´ıas se mejora el modelo ARMA para todas las horas. En la estaci´on C2-La Aldea Tabla 5.12, con un 97 % de d´ıas tipo B y C, los resultados son muy similares a los observados para la estaci´on C0. Los d´ıas tipo C muestran una clara mejor´ıa de los resultados con el modelo NN+ECMWF desde la segunda hora de predicci´on, sobre todo con respecto al modelo ARMA (mejoras de hasta el 8 % para h= 6). En los d´ıas BII no se advierte una diferencia notable entre los modelos NN y NN+ECMWF, pero ambos modelos mejoran el ARMA. En cambio, en los d´ıas tipo BIII los tres modelos pr´acticamente muestran resultados similares. En las estaciones C1-Las Palmas y C5-Sta. Br´ıgida, Tablas 5.11 y 5.13 respectivamente, la mayor´ıa del conjunto de d´ıas presentan un perfil tipo BII o BIII, hasta un 62 % en la estaci´on C1 y un 68 % en la C5. Al igual que en las estaciones del sur de la isla, los d´ıas tipo C presentan los menores errores y la mejora obtenida por el modelo NN+ECMWF es la m´as significativa. En los d´ıas tipo BII y BIII este modelo sigue presentando mejoras notables sobre todo para las ´ultimas horas de predicci´on, ya que la diferencia con los 142 CAP´ ITULO 5. APLICACI ´ ON DE LOS MODELOS DE PREDICCI ´ ON tipo d´ıa modelo 1 h 2 h 3 h 4 h 5 h 6 h BIII Persistence 32.93 44.60 49.88 53.89 55.86 55.66 SmartPersistence 32.93 43.29 46.94 47.65 45.42 43.75 ARMA 29.87 36.61 38.23 38.99 39.06 38.77 NN 30.08 35.71 37.03 36.42 37.73 37.48 NN+ECMWF 30.33 36.54 37.48 37.74 37.85 37.87 CIII Persistence 17.51 23.57 31.34 37.15 41.32 44.40 SmartPersistence 17.51 24.19 31.29 35.66 37.76 38.26 ARMA 16.15 20.39 24.20 26.13 27.03 27.38 NN 15.69 18.62 23.21 23.64 21.83 22.61 NN+ECMWF 15.01 17.47 19.94 20.22 20.39 20.57 BII Persistence 14.46 22.63 29.46 34.82 37.85 38.67 SmartPersistence 14.46 24.21 30.23 31.90 30.57 28.18 ARMA 14.62 20.22 23.00 24.19 24.33 23.97 NN 14.47 19.53 21.40 22.59 21.19 22.53 NN+ECMWF 14.30 19.95 21.41 21.88 20.94 20.15 CII Persistence 8.53 13.58 18.46 22.04 24.71 26.17 SmartPersistence 8.53 14.07 17.72 19.35 19.30 18.48 ARMA 8.16 11.69 14.21 15.68 16.60 17.08 NN 8.28 12.25 13.59 15.64 15.49 14.77 NN+ECMWF 7.68 10.47 12.59 12.53 12.37 12.49 CI Persistence 6.12 11.02 15.53 19.12 21.52 22.58 SmartPersistence 6.12 11.68 15.41 16.70 16.10 14.99 ARMA 6.41 10.03 12.43 13.89 14.70 15.07 NN 5.72 9.51 11.40 13.79 13.48 13.70 NN+ECMWF 5.55 8.72 10.43 10.60 10.47 10.61 Tabla 5.10: Errores %rRMSE en los distintos tipos de d´ıas para todos los horizontes temporales en la estaci´on C0-Pozo Izquierdo 5.2. PREDICCI ´ ON A PARTIR DE DATOS TERRESTRES Y DEL ECMWF 143 tipo d´ıa modelo 1 h 2 h 3 h 4 h 5 h 6 h AIII Persistence 47.81 64.67 73.16 78.94 82.67 86.42 SmartPersistence 47.81 64.08 70.60 74.59 76.60 76.57 ARMA 44.50 56.34 61.46 64.76 66.94 68.66 NN 44.22 57.15 62.36 65.51 68.06 68.14 NN+ECMWF 44.56 56.97 62.63 65.57 67.96 67.64 BIII Persistence 32.81 45.20 51.53 54.80 56.80 56.99 SmartPersistence 32.81 45.31 49.46 49.86 47.74 44.47 ARMA 29.98 37.61 39.94 40.41 40.37 39.94 NN 29.57 36.91 38.90 38.69 38.04 37.61 NN+ECMWF 29.63 36.46 37.76 37.57 37.44 37.16 AII Persistence 28.79 40.36 49.53 57.71 64.40 68.81 SmartPersistence 28.79 42.04 52.25 58.54 62.17 64.75 ARMA 29.56 40.29 47.10 51.88 55.21 57.46 NN 31.08 42.69 50.03 55.53 59.89 62.47 NN+ECMWF 30.55 44.14 52.59 58.59 60.86 62.49 BII Persistence 15.37 24.03 30.90 36.40 40.07 41.08 SmartPersistence 15.37 25.99 32.90 35.19 34.14 31.89 ARMA 15.40 22.01 25.69 27.67 28.49 28.46 NN 14.82 21.51 25.29 26.66 26.97 26.73 NN+ECMWF 14.70 19.95 22.20 23.47 23.93 24.36 CII Persistence 10.82 15.86 20.84 24.20 27.12 28.55 SmartPersistence 10.82 16.68 20.42 21.77 21.27 20.26 ARMA 11.55 16.52 19.99 21.92 23.20 23.87 NN 12.67 18.71 22.67 24.10 24.67 24.73 NN+ECMWF 10.37 13.12 15.11 15.73 16.89 16.06 CI Persistence 5.07 8.90 12.70 15.74 18.05 20.11 SmartPersistence 5.07 9.30 12.26 14.62 17.71 22.11 ARMA 8.28 13.75 17.62 20.20 21.94 23.16 NN 9.36 16.48 21.08 23.63 24.67 25.58 NN+ECMWF 6.65 10.58 12.67 13.85 15.82 15.30 Tabla 5.11: Errores %rRMSE en los distintos tipos de d´ıas para todos los horizontes temporales en la estaci´on C1-Las Palmas