Data Mining Models for Short Term Solar Radiation Prediction and Forecast-Based Assessment of Photovoltaic Facilities

Dpto. Lenguajes y Ciencias de la Computaci´ on Tesis Doctoral Data mining models for short-term solar radiation prediction and forecast-based assessment of photovoltaic facilities Autor: Pedro Francisco Jim´ enez P´ erez Directora: Llanos Mora L´ opez 2016 AUTOR: Pedro Francisco Jiménez Pérez http://orcid.org/0000-0001-7858-4278 EDITA: Publicaciones y Divulgación Científica. Universidad de Málaga Esta obra está bajo una licencia de Creative Commons Reconocimiento-NoComercialSinObraDerivada 4.0 Internacional: http://creativecommons.org/licenses/by-nc-nd/4.0/legalcode Cualquier parte de esta obra se puede reproducir sin autorización pero con el reconocimiento y atribución de los autores. No se puede hacer uso comercial de la obra y no se puede alterar, transformar o hacer obras derivadas. Esta Tesis Doctoral está depositada en el Repositorio Institucional de la Universidad de Málaga (RIUMA): riuma.uma.es La Dra. Llanos Mora L´ opez, Titular de Universidad en el Departamento de Lenguajes y Ciencias de la Computaci´on de la Universidad de M´alaga Certifica Que Pedro Francisco Jim´ enez P´ erez, Ingeniero en Inform´atica, ha realizado bajo su direcci´on la tesis doctoral titulada Data mining models for short-term solar radiation prediction and forecast-based assessment of photovoltaic facilities, que se recoge en la presente memoria, cumpliendo todos los requisitos legales para optar al grado de Doctor, por lo que autoriza su lectura y defensa p´ublica. Y para que as´ı conste y tenga los efectos oportunos, firmo este certificado en M´alaga, a 9 de junio de 2016 Dra. Llanos Mora L´opez A mis padres y a Priscill Este trabajo de investigaci´on ha sido financiado por la Junta de Andaluc´ıa (proyecto No. P10-TIC-6441). Agradecimientos Obtener el t´ıtulo de doctor es bastante gratificante aunque tambi´en es un proceso largo y tedioso que no ser´ıa posible sobrellevar sin la ayuda y compresi´on de toda la gente que te rodea. Es por eso que las primeras palabras de esta tesis doctoral son de agradecimiento a todas esas personas que han tenido algo que ver en el desarrollo de esta larga aventura... Mi primer agradecimiento es para Llanos, de la que destacan su gran conocimiento, su rigor cient´ıfico, sus ganas de mejorar cada aspecto de la investigaci´on, su gran disponibilidad y dedicaci´on. Sin ella esta tesis no hubiera sido posible y estoy seguro de que esta tesis es la mejor posible gracias a ella. Tambi´en agradecer a Mariano su dedicaci´on, su gran conocimiento de la materia y sus ideas que seguro hacen de esta tesis un orgullo para mi. xvi LIST OF FIGURES 5.12 Mean hourly clearness index and hourly standard deviation for clusters 3 (left) and 4 (right) . . . . . . . . . . . . . . . . . . . . . . . . 99 5.13 Mean hourly clearness index and hourly standard deviation for clusters 5 (left) and 6 (right) . . . . . . . . . . . . . . . . . . . . . . . . 100 6.1 Daily profiles of hourly values of clearness index (left) and hourly values of k∗ h,d (right) for two different days . . . . . . . . . . . . . . 107 6.2 Scheme of the proposed procedures to obtain models to characterize hourly global solar radiation. In procedure 1, only one model is selected to estimate kd, regardless of cluster. In procedure 2, one model for each cluster is selected to estimate kd.∗training set . . . 109 6.3 Proposed procedure for forecasting next-day values of hourly solar radiation. ................................111 6.4 Values of index k∗ h,d (centroids) for each cluster. . . . . . . . . . . . 114 6.5 Values of index k∗ h,d for the days included in each cluster. Each line corresponds to the values of this index for a day. . . . . . . . . . . . 115 7.1 Main system architecture . . . . . . . . . . . . . . . . . . . . . . . . 126 7.2 Instantaneous values of solar radiation and photovoltaic power . . . 128 7.3 Flowchart for assessing solar energy plants. . . . . . . . . . . . . . . 131 7.4 Dailyparameters ............................132 7.5 Detecting problems in the inverter . . . . . . . . . . . . . . . . . . . 133 7.6 Detecting problems in one plant with two inverters . . . . . . . . . 134 A.1 Junta de Andaluc´ıa Web Site for meteorological data . . . . . . . . 144 A.2 Web application, step 1 . . . . . . . . . . . . . . . . . . . . . . . . . 145 A.3 Web application, step 2 . . . . . . . . . . . . . . . . . . . . . . . . . 146 xvii Table Index 1.1 Renewable power capacity until 2014 . . . . . . . . . . . . . . . . . 7 1.2 Energy demand evolution 2011-2015 . . . . . . . . . . . . . . . . . . 8 1.3 Renewable installed capacity and power production in 2015 . . . . . 8 1.4 Renewable installed capacity and power production in 2015 . . . . . 9 2.1 Results for different models tested in southern Spain and Italy . . . 29 3.1 Sky cover codes and descriptions from AEMET . . . . . . . . . . . 68 4.1 Measure types and units . . . . . . . . . . . . . . . . . . . . . . . . 74 5.1 Results of Kolmogorov-Smirnov two sample test for different numberofclusters.............................. 95 5.2 Energy error (EE)(%) for each cluster when forecasting hourly solar global radiation with proposed model and with a persistent naive model and percentage of energy received in all days included in each cluster respect to total energy received.(∗clusters are randomly built)100 6.1 Independent variables used to build the models (cols 1 and 2) and forecasting input variables. Note: d−1 refers to the previous day and dcurrent day, 12−15 refers to the mean value estimated using values from 12 to 15 hours, and 9−12 is for the mean value estimated using values from 9 to 12 hours, ˆ Xrefers to the forecasted value of X,Tis temperature, His humidity and Patmospheric pressure. . 112 6.2 Number of observations in each cluster and percentage of radiation received in days included in cluster . . . . . . . . . . . . . . . . . . 113 xviii TABLE INDEX 6.3 Training set: errors for each analyzed combination of methods for Experiment 1 (values of meteorological variables for the previous day) and Experiment 2 (values of meteorological variables for the same day to be modeled are the independent variables except the value of daily clearness index). (1) The same model is fitted for all observations (Procedure 1). (2) One model is fitted for the observations of each cluster (Procedure 2). . . . . . . . . . . . . . . . . . . 116 6.4 Test sets: errors for each analyzed combination of methods for Experiment 1 (values of meteorological variables for the previous day) and Experiment 2 (forecasts of meteorological variables for the same day to be modeled are the independent variables except the value of daily clearness index). (1) The same model is fitted for all observations (Procedure 1). (2) One model is fitted for the observations of each cluster (Procedure 2). . . . . . . . . . . . . . . . . . . . . . 117 6.5 Training set: errors by cluster for each analyzed combination of methods for Experiment 1 (values of meteorological variables for the previous day) and Experiment 2 (values of meteorological variables for the same day to be modeled are the independent variables except the value of daily clearness index). (1) The same model is fitted for all observations (Procedure 1). (2) One model is fitted for the observations of each cluster (Procedure 2). . . . . . . . . . . . . . . 117 6.6 Test sets: errors by cluster for each analyzed combination of methods for Experiment 1 (values of meteorological variables for the previous day) and Experiment 2 (forecasts of meteorological variables for the same day to be modeled are the independent variables except the value of daily clearness index). (1) The same model is fitted for all observations (Procedure 1). (2) One model is fitted for the observations of each cluster (Procedure 2). . . . . . . . . . . . . 118 6.7 Comparison of the performance of the SVM-C+SVM-R models. . . 119 1 Nomenclature αSolar elevation δDeclination γparameter for the kernel function in SVM ˆ AEstimated value of A λLongitude ω0Constant that determines the activation threshold ωhHour angle ωiWeight of input i ωsr Sunrise angle φLatitude dday E0Eccentricity factor GtAverage values of measured global horizontal irradiation for period t G0,d Daily extraterrestrial global radiation on a horizontal surface G0,h Hourly extraterrestrial global radiation on a horizontal surface G0,t Extraterrestrial global horizontal radiation on period t HHumidity hHour Isc Solar constant k∗ h,d Daily-detrended hourly clearness index kdDaily clearness index khHourly clearness index ktClearness index for period t 2TABLE INDEX MAE Mean absolute error MSE Mean square error PPresure rMAE Relative mean absolute error RMSE Root mean square error sForecast skill over 24 hour persistence forecast Ttemperature xValue of an input or of a neuron of previous layer in ANN ANN Artificial neural network DT Decision tree SVM Support vector machine SVM-C Support vector machine for classification SVM-R Support vector machine for regression 3 Chapter 1 Introduction As society progresses, the energy needs grow, we have more electronic devices that consume electricity, we travel to more places and more distant, more sophisticated goods are manufactured and factories that produce them require more energy. General energy consumption increases, especially the consumption of electric energy, devices that have traditionally worked with other energy sources now are electrifying itself, like vehicles (that worked with oil), bicycles (which traditionally had pedals), kitchens (traditionally powered with gas), etc. Energy is a strategic resource and countries seek to be energy independent since an energy dependent country has to pay a large bill for the energy consumed and therefore its development is slowed considerably. Integrating energy generation in factories and urban environments is interesting because that is where the most energy is consumed and avoid wasting energy in transport infrastructure as well as saving. Renewable energies are often the easiest to integrate into these environments because they require less infrastructure and cause fewer problems with noise, dirt, pollution, etc. Traditional energy sources have serious problems: they are limited, produce harmful pollution on health and the environment and someday (fossil fuels) they will run-out, so mankind is forced to find clean and renewable energy sources. At the same time scientific and technical developments advance, society is moving towards a more acute awareness of his responsibility to the environment. This is one fundamental reason that today renewable energy is undergoing a huge 41. INTRODUCTION expansion, which raises numerous challenges to the scientific and technical community. Renewable energy is clean, inexhaustible, and its technology is developing quickly thanks to initiatives as the Kyoto protocol, global commitment against global climate change. Renewable energy sources include: wind power, hydropower, solar energy, biomass, biofuel, etc. The motivation for searching and developing clean energy sources comes from several facts: the greenhouse gases effect over environment (temperature raising, rain and snow decrease), (Core Writing Team and , eds.), depletion of some energy sources like gas and oil, (Bentley, 2002), the oil price behaviour (Gori et al., 2007), the oil resources, (Kj¨arstad and Johnsson, 2009), (Owen et al., 2010) and (Hughes and Rudolph, 2011). In this scene, solar energy, and more specifically the photovoltaic solar energy, is beginning to be increasingly important in the energy mix. 1.1 Photovoltaic solar energy current situation The number of solar energy plants has increased significantly in recent years, mainly due to the following factors: the need to use energy sources that contribute to reducing carbon emissions, establishing support policies to introduce this type of systems, improving the efficiency of these systems and the significant reduction in the price of all the components that make up those systems. As the number of these systems rises, there is an increasingly greater need to develop systems that enable these energy sources to be integrated with the traditional generation system. Renewable energy continued growing in 2014 at the same time of global energy consumption and despite the dramatic decline in oil prices during the second half of the year (REN21, 2015). The biggest growth took place in the power sector and was dominated by three technologies: wind, solar photovoltaic and hydropower. Solar PV and wind power have suffered significant cost reductions that have played an important role in the increasing electrification of transportation and heating applications and highlighting the potential for these two technologies in the future. In many countries renewable energy is broadly competitive with conventional fuels, particularly in the power sector. In developing countries distributed renewable energy systems offer a great opportunity to speed up the transition to moderm 1.1. PHOTOVOLTAIC SOLAR ENERGY CURRENT SITUATION 5 energy services and to grant energy access, although prices and financing are still being one of the major barriers. Renewables have become a mainstream energy resource, in some cases with a rapid expansion, contributing to diversification of the energy mix, however, growth in renewables capacity and improvements in energy efficiency are below the necessary rates to achieve the Sustainable Energy for All (SE4ALL) goals of doubling the renewable energy levels and energy efficiency improvements in order to provide universal energy access by 2030. It’s not clear what a distributed energy generation system is. Different definitions for Distributed Generation (DG) can be found in literature (R´ujula et al., 2005), and some of then are not consistent, probably because they derive from different industries, experiences and objectives, but it’s a widely used term. Despite rising energy use, global carbon emissions associated with energy consumption remainded stable in 2014 due to the increased penetration of renewable energy and to improvements in energy efficiency. There is a global growing awareness that deployment of renewable energy is critical for addressing climate change, creating new economic opportunities and providing energy access to billions of people that still living without modern energy services. Renewable energy provided an estimated 19.1% of the total energy consumed in 2013 and continued growing in 2014. Renewable energy support policies and decreasing costs are the main factors that have encouraged growth of renewable energy sources. At the same time, subsidies to fossil fuels and nuclear power are an important obstacle in developing countries. China was the first country in new renewable power capacity installations in 2014 while Brazil, India and South Africa led for the capacity added in their respective regions. Renewables accounted for 58.5% of net additions to global power capacity in 2014 and comprised an estimated 27.7% of the world’s power generating capacity, enough to produce about 22.8% of global electricity, see figure 1.1. Falling costs have made unsubsidised solar PV-generated electricity cost-competitive with fossil fuels. In 2014, solar PV reached another important milestone with an estimated 40 GW installed for a total global capacity of about 177 GW. Figure 1.1 shows the evolution of the different renewable technologies over the last 10+ years, including solar PV. The major part of this new capacity was installed in China, Japan and United States. Latin America and African countries had a significant new capacity added while most EU markets continued to decline for the third consecutive year. Also, 2014 saw significants improvements in energy storage systems accross all 61. INTRODUCTION 78% 22% Fossils fuels and nuclear Renewable electricity Figure 1.1: Global electricity supply in 2013, Source: GSR2015 sectors and several regions have seen significant growth in numbers of residential “prosumers” (electricity customers who produce their own power). In general, the solar PV industry recovery that began in 2013 continued in 2014, thanks to a strong global market. According to data from (for Solar Energy Systems, 2015) photovoltaics is a fast growing market as the compound annual growth rate (CAGR) of PV installations was 44% from 2000 through 2014. China and Taiwan hold the lead of PV module production in 2014 with a share of 69% while Europe contributed with a 6% and USA and Canada each contributed with a 4%. In 2014 Europe contributed with 48% of the total cumulative PV installations (in 2013 it was 58%) while China in Taiwan accounted for 17% (compared to 13% in 2013). 92% of the total production in 2014 was for Si-wafer PV technology. The share of multi-crystalline technology is now about 56% of total production and the share of all thin film technologies amounted to about 9%. The best cell efficiency is 25.6% for mono-crystalline and 20.8% for multi-crystalline silicon wafer-based technology while the best efficiency in thin film technology is 21.0% for CdTe and 20.5% for CIGS solar cells. In the last 10 years, the efficiency of average commercial wafer-based silicon modules increased from about 12% to 16% and CdTe module efficiency increased from 9% to 13%. The Energy Payback Time of PV systems is dependent on the geographical location: while PV systems in northern Europe need around 2.5 years to balance 1.1. PHOTOVOLTAIC SOLAR ENERGY CURRENT SITUATION 7 START 2004 2013 2014 INVESTMET New investment in renewable power and fuels Billion USD 45 232 270 POWER Renewable power capacity (total, not including hydro) GW 85 560 657 Renewable power capacity (total, including hydro) GW 800 1578 1712 Hydropower capacity (total) GW 715 1018 1055 Bio-power capacity GW ¡36 88 93 Bio-power generation TWh 227 396 433 Geothermal power capacity GW 8.9 12.1 12.8 Solar PV capacity (total) GW 2.6 138 177 Concentrating solar thermal power (total) GW 0.4 3.4 4.4 Wind power capacity (total) GW 48 319 370 HEAT Solar hot water capacity (total) GWth 86 373 406 TRANSPORT Ethanol production (annual) billion litres 28.5 87.8 94 Biodiesel production (annual) billion litres 2.4 26.3 29.7 Table 1.1: Renewable power capacity until 2014 the input energy, PV systems in the south equal their energy input after 1.5 years and less, depending on the technology installed. For example, a PV system located in Sicily made with multi-Si modules has an Energy Payback Time of around one year, so, assuming 20 years of useful life, this system can produce twenty times the energy needed to produce it. Inverter efficiency actually stands at 98% and above. The market share of inverters is estimated to be 50% for string inverters (residential, small and medium commercial applications), 48% for central inverters (with applications mostly in large commercial and utility-scale systems) and about 1.5% belongs to microinverters (used on the module level). PV systems prices have reduced dramatically in last 20-25 years, for example, in Germany prices for a typical 10 to 100 kWp PV rooftop-system were around 14,000 e/kWp in 1990 while at the end of 2014, such systems cost about 1,300 e/kWp, showing a net-price regression of about 90% in that period of time (equivalent to an annual compound average price reduction rate of 9%). Nowadays China is the most important PV cell and module manufacturer, with about 60% of the total annual production in 2014. Europe and USA were big manufacturers 20 years ago but now they have a little share due to the extraordinary growth of China production. Europe and Asia are the regions with most cumulative installed PV systems followed by America in third position. Germany accounts for the 20% of the world total, Italy accounts for the 10% and the rest of Europe for the 18%. Multi-Si is 15 Bibliography Bentley, R., 2002. Global oil and gas depletion: an overview. Energy Policy 30 (3), 189 – 205. BOE, 2007. Bolet´ın oficial de estado. real decreto 661/2007, de 25 de mayo, por el que se regula la producci´on de energ´ıa el´ectrica en r´egimen especial. BOE, 2010. Bolet’´ın oficial del estado. real decreto 1565/2010, de 19 de noviembre, por el que se regulan y modifican determinados aspectos relativos a la actividad de producci´on de energ´ıa el´ectrica. Core Writing Team, R. P., (eds.), L. M., 2014. Ipcc, 2014: Climate change 2014: Synthesis report. contribution of working groups i, ii and iii to the fifth assessment report of the intergovernmental panel on climate change. Tech. rep., IPCC, Geneva, Switzerland. for Solar Energy Systems, F. I., 2015. Photovoltaics report. Tech. rep., Fraunhofer Institute for Solar Energy Systems, ISE. Gori, F., Ludovisi, D., Cerritelli, P., 2007. Forecast of oil price and consumption in the short term under three scenarios: Parabolic, linear and chaotic behaviour. Energy 32 (7), 1291 – 1296. URL http://www.sciencedirect.com/science/article/pii/ S0360544206001873 Hughes, L., Rudolph, J., 2011. Future world oil production: growth, plateau, or peak? Current Opinion in Environmental Sustainability 3 (4), 225 – 234, energy Systems. URL http://www.sciencedirect.com/science/article/pii/ S1877343511000509 Kj¨arstad, J., Johnsson, F., 2009. Resources and future supply of oil. Energy Policy 37 (2), 441 – 464. 16 BIBLIOGRAPHY URL http://www.sciencedirect.com/science/article/pii/ S0301421508005259 Owen, N. A., Inderwildi, O. R., King, D. A., 2010. The status of conventional world oil reserves—hype or cause for concern? Energy Policy 38 (8), 4743 – 4749. URL http://www.sciencedirect.com/science/article/pii/ S0301421510001072 REN21, 2015. Renewables 2015 global status report key findings. Tech. rep. R´ujula, A. B., Amada, J. M., Bernal-Agust´ın, J., Loyo, J. Y., Navarro, J. D., 2005. Definitions for distributed generation: a revision. Proceedings of the International Conference on Renewable Energy and Power Quality 05. 17 Chapter 2 State of the Art 2.1 Introduction This chapter reviews other works about solar radiation forecasting, including some located in the mediterranean zone. Performance forecasting of systems that use solar radiation as energy resource requires, on the one hand, using actual forecasted running conditions (meteorological parameters) and, on the other hand, using models that permit estimating the operation performance of these systems taking into account these conditions. Estimating the energy produced by solar plants is difficult mainly due to its dependence on meteorological variables, such as solar radiation and temperature, (Luque and Hegedus, 2002), (Chang, 2009). In fact, photovoltaic production prediction is mainly based on global solar irradiation forecasts. The main problem of having predictions of solar radiation that these systems will receive is in the nature of the solar resource, since it is an intermittent resource (due to day-night succession), it has a seasonal component, due to changes of relative position between Sun and Earth and has a certain stochastic behavior. The behavior of this variable can change dramatically on different days, even on the same day, due to the stochastic nature. Prediction models of solar radiation must be able to collect these trends and to reproduce the non-deterministic component of it. In the following sections some of the main approaches used to predict solar radiation are analyzed and described. 18 2. STATE OF THE ART It also includes a specific section in which the various errors of the main models are shown to have a reference to compare the proposals made in this work with. 2.2 Solar Radiation Forecasting Models Multiple methods like statistical and data mining techniques can be used to address the need to forecast solar radiation. First attempts in irradiance forecasting were made more than twenty years ago (Jensenius and Cotton, 1981), when daily solar radiation forecasts for one to two days in advance have been produced with the Model Output Statistics (MOS) technique (Glahn and Lowry, 1972). Subsequent years showed only minor attemps or progress with respect to the development of solar irradiance forecasting methods. (Heck and Takle, 1987) and (Jensenius, 1989) both presented variations of the MOS approach without introducing new elements. After these first approaches, numerous works and approaches have been proposed to characterize and predict solar radiation. Some of these approaches are based on using physical models, others, conversely, are based on the assumption that the prediction of solar radiation can be tackled as others forecasting processes. Here is a review of studies about solar irradiance forecasting using different methods. 2.3 Statistical and data mining models Statistical methods based on historical data can be divided in two categories: statistical and data mining models. Examples of statistical methods are: seasonality analysis, Auto Regressive Moving Average (ARMA). Examples of data mining models are fuzzy inference, genetic algorithm and neural networks. 2.3.1 Statistical models Statistical time series models are based on the assumption that the data have an internal structure and it can be identified by using simple and partial autocor- 2.3. STATISTICAL AND DATA MINING MODELS 19 relation, (Box and Jenkins, 1976), (Gooijer and Hyndman, 2005), (Brockwell and Davis, 2002). Time series forecasting methods detect and explore such a structure. In particular, ARMA (autoregressive moving average), ARIMA (autoregressive integrated moving average) models have been widely used; for instance, (Brinkworth, 1997), (Bartoli et al., 1983), (Aguiar et al., 1988), (Graham et al., 1988), (Aguiar and Collares-Pereira, 1992) and (Mora-L´opez and de Cardona, 1998) propose different methods for modeling hourly and daily series of clearness index (parameter related to solar global radiation). These models are particularly useful for long-term characterization and prediction of the clearness index as they pick up the statistical and sequential properties of series. However these methods have not been used for short-term prediction of clearness index as the error in the prediction of isolated values (next value in a series) is too large. (Mora-Lopez and Sidrach-de Cardona, 1998) proposes a methodology to generate hourly series of global irradiation. The only input parameter required is the monthly mean value of daily global irradiation, which is not difficult to be available for many locations. The procedure is based on multiplicative autoregressive moving-average (ARMA) statistical model for time series with regular and seasonal components and is able to capture the two relationships observed in recorded hourly series of global irradiation: the relationship between one value and the value of the previus hour and the relationship of the value and the value of the previous day at the same hour. Data from several Spanish cities and ranging from year 1976 to 1986 are used. Daily and seasonal trends are removed using the maximum hourly global irradiation value and a difference operator with hourly values from a given day and the prevoius day. In order to check the validity of the proposed methods to estimate the unknown parameters of the multiplicative ARMA models, synthetic hourly series of clearness index and global irradiation were obtained and it was verified that the generated series (for both parameters) have the same statistical characteristic as the real series, that is, same mean, variance and cumulative probability distribution function as the real series (using 0.05 as significance level in Kolmogorov-Smirnov test) In (Reikard, 2009) six data sets (three from Kansas City (Missouri), Denver (Colorado), and Phoenix (Arizona), and the others from the Measurements and Instrumentation Data Center baseline measurement system database: Clark power station in Nevada, Solar Radiation Research Laboratory (SRRL) and The National Wind Technology Laboratory) are used to run forecasting experiments at resolutions of 5, 15, 30, and 60 min, using the global horizontal component. Forecasting tests are run using regressions in logs, Autoregressive Integrated Moving Average (ARIMA), and Unobserved Components models but also transfer functions, neural networks, and hybrid models are evaluated. The models are estimated over 20 2. STATE OF THE ART history prior to the start of the forecast horizon, the data is forecasted and then the predicted values are compared with the actuals. The best results are obtained using the ARIMA in logs in nearly all the tests, with time-varying coefficients, but with some exceptions: at high resolutions, a transfer function using cloud cover is found to improve over the ARIMA and in a few cases, the neural net or hybrid models can improve at very high resolutions (5 min). The ability to capture the diurnal cycle more effectively than other methods explains the success of the ARIMA method. A new approach divided into two phases is used to predict the hourly solar radiation series in (Ji and Chee, 2011). The Autoregressive and Moving Average (ARMA) model is used to predict the stationary residual series, previously detrended selecting the best model based on the Augmented Dickey–Fuller method to test the stationarity of the residual. Furthermore, a hybrid model that combines both the ARMA and Time Delay Neural Network (TDNN), is applied to produce better prediction, where ARMA model is used to predict the linear component of the series and the TDNN model is used to predict the nonlinear component. The simulation shows that this hybrid model can take the advantages of both ARMA and TDNN and give better results than applying only the ARMA model. As explained in (Diagne et al., 2013), statistical models based on online irradiance measurements can be used to forecast at very short term, from 5 minutes up to 6 hours. Some examples are Auto Regressive (AR) and Auto Regressive Moving Average (ARMA) models. 2.3.2 Data mining models Data mining techniques have been also applied for process forecasting. These approaches have been proposed to overcome the limitations of statistical methods because they do not require any assumptions to be made, particularly with respect to the linearity of the series. In (Sfetsos and Coonick, 2000) a new approach for the forecasting of mean hourly global solar radiation on a horizontal surface is introduced. Arfificial Neural Networks (ANN) and the Adaptative Neuro-fuzzy Inference Scheme (ANFIS) are used. Initially one variable is used and then experiment is extended to include aditional meteorological parameters. Results indicate better performance for this artifical intelligence model over conventional procedures based on clearness index and are able to capture de periodic nature of this series. Levenberg- 2.3. STATISTICAL AND DATA MINING MODELS 21 Marquardt (LM) network was found to be the optimal prediction model (against back-propagation). Use of aditional meteorological variables such as temperature, pressure, wind direction and speed as potential input for the forecasting process is researched; some models can be further enhanced using aditional meteorological parameters but not all. LM and ANFIS performance are enhanced when using wind direction parameter as input. The best prediction is found to be that from the multivariate LM case, with an RMS error improvement about 74% compared with that of the bench-mark persistent aproach. In this approach no transformations of the solar radiation values are required. In (Perez et al., 2007) a simple solar forecast model using sky cover predictions is developed and tested against both ground-measured and satellite-derived irradiances data. Data is taken from the National Digital Forecast Database [NDFD] of United States National Weather Service, providing gridded forecasted parameters for the entire country. Three methods are tested for different forecasting periods, resulting that ’Best fit formula’ is the best method with a Relative Mean Bias Error of -2% and a Relative Root Mean Square Error of 35% for less than 4 hours ahead forecasting. Artificial neural network and ARIMA models are proposed in (Reikard, 2009); the errors range from 30 to 40% in energy terms. Similar models are used in (Voyant et al., 2013); the obtained errors for predicting hourly values for a day range from 23% to 28%. (Mellit and Pavan, 2010) proposes a practical method for solar irradiance forecast using artificial neural network (ANN), a Multilayer Perceptron MLP-model that makes it possible to forecast the solar irradiance on a base of 24 hours using the present values of the mean daily solar irradiance and air temperature. A database of solar irradiance and air temperature data (from July 1st 2008 to May 23rd 2009 and from November 23rd 2009 to January 24th 2010) collected in Trieste, Italy, is used. A K-fold cross-validation was carried out in order to check the generalization capability of the MLP-forecaster. The results indicate a good performance, as the correlation coefficient is in the range 98–99% for sunny days and 94–96% for cloudy days. A comparison between the forecasted energy and the one produced by the GCPV plant installed on the rooftop of the municipality of Trieste shows the goodness of the proposed model. The MLP input layer accepts the mean daily solar irradiance, the mean daily air temperature and the day of the month (at the time t) as parameters, while the output layer gives the 24 h of solar irradiance at the next day (time t + 1). The results show that the developed MLP-forcaster is suitable for the prediction of sunny days with rvalues between 98% and 99%), and it also provides acceptable results for cloudy days where ris 22 2. STATE OF THE ART between 92% and 95%. For the whole dataset the MBE and the RMSE are 32% and 67% respectively. In (Mora-L´opez et al., 2011) a model for short-term forecasting of continuous time series of solar radiation that binds the use of both statistical (regression) and machine learning methods has been performed. The mean square errors of the proposed models range from 0.04 to 0.4 depending on the value of clearness index. An artificial neural network (ANN) model was used to estimate the solar radiation parameters for seven cities from the mediterranean region of Anatolia in Turkey in (Koca et al., 2011). The maximum RMSE was found to be 6.9% for Mersin citation and the best value was obtained to be 3.6% for Isparta. A medium-term solar irradiance forecasting model was developed in (Marquez and Coimbra, 2011) adopting predicted meteorological variables from the US National Weather Service’s (NWS) forecasting database as inputs to an Artificial Neural Network (ANN) model. The inputs involved are the same from a validated forecasting model so mean bias error (MBE), root mean square error (RMSE) and correlation coefficient (R2) comparisons between the more established forecasting model and the proposed one are included. A set of criteria for selecting relevant inputs was developed, input variables were selected using a version of the Gamma test combined with a genetic algorithm. The solar geotemporal variables were found to be critically important, while the most relevant meteorological variables included sky cover, probability of precipitation, and maximum and minimum temperatures. Using the relevant input sets identified by the Gamma test, the developed forecasting models improve RMSE for GHI by 10–15% over the reference model. rRMSE range from 15% to 22% for different models constructed on 13 month data set for same-day forecasts of GHI. In (Wang et al., 2012) an ANN model using statistical feature parameters (ANN-SFP) for short-term solar irradiance forecasting is proposed where the input vector is reconstructed with several statistical feature parameters of irradiance and ambient temperature that include (i) daily maximum value of the third order derivative of the difference between surface global irradiance and extraterrestrial global irradiance and (ii) normalized discrete difference of solar surface and extraterrestrial irradiances. These statitistical features help in reducing the model complexity because it uses five inputs. The network training is done using the Levenberg-Marquardt algorithm (LMA). Network output are 24, one for each hour of the day with the predicted irradiance value. After simulations are carried out, the proposed model is validated and compared with the conventional ANN model using historical data series (ANN-HDS) dummy model, and the results showed 2.4. CLOUD IMAGINERY 23 an inproved accuracy under variable weather conditions. The proposed model is validated using solar irradiance recorded data at Yundian Science and Technology Park grid-connected PV plant from March 2011 to December 2011. 2.4 Cloud imaginery Using information about temporal developments of clouds may help in solar irradiance forecasting. Good performance is achieved from 30 minutes up to 6 hours using cloud motion vectors from satellite images. Cloud information from ground based sky images may be used to get irradiance forecasts with much higher spatial and temporal resolution. Besides the deterministic factors (Earth movement relative to Sun position) clouds are the main influence factor over the changing values of solar irradiance. Clouds show a high variability in time so determining clouds position at future time is essential to forecasting accuracy. Satellite and groud-based sky images are available techniques to forecast and model solar irradiance. The basis of this method relies upon determination of clouds structure at previous time and determine future position using extrapolation. With it’s high spatial and temporal resolution satellite and ground-based sky images offer the potential to derive the required information on cloud motion using motion vector fields. With this method, clouds position and structure can be determined up to 6 hours in advance. Using this information, irradiance for all sky conditions including cloudy skies maybe derived using radiative transfer models (RTM) (Diagne et al., 2013). 2.4.1 Satellite images In (Hammer et al., 1999) satellite images are used to forecast solar surface irradiance. The satellite data provides information about cloudiness which is used to detect motion of cloud structures. Extrapolating the temporal development of the clouds, solar irradiance can be predicted from 30 minutes up to 2 hours. Images used are taken from METEOSAT satellite and an enhanced version of semiempirical HELIOSAT method is used to derive surface irradiance. Test shown a certain improvement in forecast accuracy over persistence forecast. A statistical fit of the relationship between a normalised parameter of the solar 30 2. STATE OF THE ART key issue. 31 Bibliography Aguiar, L. M., Pereira, B., David, M., D´ıaz, F., Lauret, P., 2015. Use of satellite data to improve solar radiation forecasting with bayesian artificial neural networks. Solar Energy 122, 1309 – 1324. URL http://www.sciencedirect.com/science/article/pii/ S0038092X15005927 Aguiar, R., Collares-Pereira, M., 1992. T.a.g: A time dependent autoregressive gaussian model for generating synthetic hourly radiation. Solar Energy 49(3), 167–174. Aguiar, R., Collares-Pereira, M., Conde, J., 1988. Simple procedure for generating sequences of daily radiation values using a library of markov transition matrices. Solar Energy 4. (3), 269–279. Bartoli, B., Coluaai, B., Cuomo, V., Francesca, M., Serio, C., 1983. Autocorrelation of daily global solar radiation. Il nuovo cimento 40, 113–122. Box, G., Jenkins, G., 1976. Time Series Analysis forecasting and control. Prentice Hall. Brinkworth, B., 1997. Autocorrelation and stochastic modelling of insolation sequences. Solar Energy 19, 343–347. Brockwell, P., Davis, R. A., 2002. Introduction to Time Series and Forecasting. Springer Texts in Statistics. Cazorla, A., Olmo, F. J., Alados-Arboledas, L., Jan 2008. Development of a sky imager for cloud cover assessment. J. Opt. Soc. Am. A 25 (1), 29–39. URL http://josaa.osa.org/abstract.cfm?URI=josaa-25-1-29 Chang, T., 2009. Output energy of a photovoltaic module mounted on a single-axis tracking system. Applied Energy 86, 2071–2078. 32 BIBLIOGRAPHY Chow, C. W., Urquhart, B., Lave, M., Dominguez, A., Kleissl, J., Shields, J., Washom, B., 2011. Intra-hour forecasting with a total sky imager at the {UC} san diego solar energy testbed. Solar Energy 85 (11), 2881 – 2893. URL http://www.sciencedirect.com/science/article/pii/ S0038092X11002982 Diagne, M., David, M., Lauret, P., Boland, J., Schmutz, N., 2013. Review of solar irradiance forecasting methods and a proposition for small-scale insular grids. Renewable and Sustainable Energy Reviews 27, 65 – 76. URL http://www.sciencedirect.com/science/article/pii/ S1364032113004334 Glahn, H. R., Lowry, D. A., 1972. The use of model output statistics (mos) in objective weather forecasting. Journal of applied meteorology 11 (8), 1203–1211. Gooijer, J. G. D., Hyndman, R. J., 2005. 25 years of iif time series forecasting: A selective review. Monash Econometrics and Business Statistics Working Papers 12/05, Monash University, Department of Econometrics and Business Statistics. URL http://econpapers.repec.org/RePEc:msh:ebswps:2005-12 Graham, V., Hollands, K., Unny, T., 1988. A time series model for kt with application to global synthetic weather generation. Solar Energy 40, 83–92. Hammer, A., Heinemann, D., Lorenz, E., L¨uckehe, B., 1999. Short-term forecasting of solar radiation: a statistical approach using satellite data. Solar Energy 67 (1–3), 139 – 150. URL http://www.sciencedirect.com/science/article/pii/ S0038092X00000384 Heck, P., Takle, E., 1987. Objective forecasts of solar radiation and temperature. Iowa State Journal of Research 62, 29–42. HIRLAM, 2016. http://hirlam.org. Jensenius, J., 1989. Insolation forecasting. Solar Resources, MIT Press, Cambridge, 335–349. Jensenius, J., Cotton, G., 1981. The development and testing of automated solar energy forecasts based on the model output statistics (mos) technique. In: 1st Workshop on terrestrial solar resource forecasting and on use of satellites for terrestrial solar resource assessment, Washington, DC. Ji, W., Chee, K. C., 2011. Prediction of hourly solar radiation using a novel hybrid model of {ARMA}and {TDNN}. Solar Energy 85 (5), 808 – 817. BIBLIOGRAPHY 33 URL http://www.sciencedirect.com/science/article/pii/ S0038092X11000259 Koca, A., Oztop, H. F., Varol, Y., Koca, G. O., 2011. Estimation of solar radiation using artificial neural networks with different input parameters for mediterranean region of anatolia in turkey. Expert Systems with Applications 38 (7), 8756 – 8762. URL http://www.sciencedirect.com/science/article/pii/ S0957417411001059 Kostylev, V., Pavlovski, A., et al., 2011. Solar power forecasting performance– towards industry standards. In: 1st International Workshop on the Integration of Solar Power into Power Systems Aarhus, Denmark. Lara-Fanego, V., Ruiz-Arias, J., Pozo-V´azquez, D., Santos-Alamillos, F., TovarPescador, J., 2012. Evaluation of the {WRF}model solar irradiance forecasts in andalusia (southern spain). Solar Energy 86 (8), 2200 – 2217, progress in Solar Energy 3. URL http://www.sciencedirect.com/science/article/pii/ S0038092X11000582 Lorenz, E., Remund, J., M¨uller, S. C., Traunm¨uller, W., Steinmaurer, G., Pozo, D., Ruiz-Arias, J. A., Fanego, V. L., Ramirez, L., Romeo, M. G., et al., 2009. Benchmarking of different approaches to forecast solar irradiance. In: 24th European photovoltaic solar energy conference, Hamburg, Germany. Vol. 21. p. 25. Luque, A., Hegedus, S., 2002. Handbook of photovoltaic science and engineering. John Wiley & Sons Ltd., Berlin. Marquez, R., Coimbra, C. F., 2011. Forecasting of global and direct solar irradiance using stochastic learning methods, ground experiments and the {NWS}database. Solar Energy 85 (5), 746 – 756. URL http://www.sciencedirect.com/science/article/pii/ S0038092X11000193 Mellit, A., Pavan, A. M., 2010. A 24-h forecast of solar irradiance using artificial neural network: Application for performance prediction of a grid-connected {PV}plant at trieste, italy. Solar Energy 84 (5), 807 – 821. URL http://www.sciencedirect.com/science/article/pii/ S0038092X10000782 Mora-L´opez, L., de Cardona, M. S., 1998. Multiplicative arma models to generate hourly series of global irradiation. Solar Energy 63, 283–291. 34 BIBLIOGRAPHY Mora-Lopez, L., Sidrach-de Cardona, M., 1998. Multiplicative arma models to generate hourly series of global irradiation. Solar Energy 63 (5), 283–291. Mora-L´opez, L., Mart´ınez-Marchena, I., Piliougine, M., Sidrach-de Cardona, M., 2011. Binding statistical and machine learning models for short-term forecasting of global solar radiation. In: Proceedings of the 10th international conference on Advances in intelligent data analysis X. IDA’11. Springer-Verlag, Berlin, Heidelberg, pp. 294–305. URL http://dl.acm.org/citation.cfm?id=2075337.2075367 Perez, R., Lorenz, E., Pelland, S., Beauharnois, M., Knowe, G. V., Jr., K. H., Heinemann, D., Remund, J., M¨uller, S. C., Traunm¨uller, W., Steinmauer, G., Pozo, D., Ruiz-Arias, J. A., Lara-Fanego, V., Ramirez-Santigosa, L., Gaston-Romero, M., Pomares, L. M., 2013. Comparison of numerical weather prediction solar irradiance forecasts in the us, canada and europe. Solar Energy 94, 305 – 326. URL http://www.sciencedirect.com/science/article/pii/ S0038092X13001886 Perez, R., Moore, K., Stackhouse, P., 2007. Forecasting solar radiation preliminary evaluation of an approach based upon the national forecast database. Solar Energy 81(6), 809–812. Pierro, M., Bucci, F., Cornaro, C., Maggioni, E., Perotto, A., Pravettoni, M., Spada, F., 2015. Model output statistics cascade to improve day ahead solar irradiance forecast. Solar Energy 117, 99 – 113. URL http://www.sciencedirect.com/science/article/pii/ S0038092X15002212 Reikard, G., 2009. Predicting solar radiation at high resolutions: A comparison of time series forecasts. Solar Energy 83 (3), 342 – 349. URL http://www.sciencedirect.com/science/article/pii/ S0038092X08002107 Sfetsos, A., Coonick, A., 2000. Univariate and multivariate forecasting of hourly solar radiation with artificial intelligence techniques. Solar Energy 68 (2), 169 – 178. URL http://www.sciencedirect.com/science/article/pii/ S0038092X9900064X Traunm¨uller, W., Steinmaurer, G., 2010. Solar irradiance forecasting, benchmarking of different techniques and applications of energy meteorology. In: Proceedings of the EuroSun 2010 conference. BIBLIOGRAPHY 35 Voyant, C., Paoli, C., Muselli, M., Nivet, M.-L., 2013. Multi-horizon solar radiation forecasting for mediterranean locations using time series models. Renewable and Sustainable Energy Reviews 28 (0), 44 – 52. URL http://www.sciencedirect.com/science/article/pii/ S1364032113005030 Wang, F., Mi, Z., Su, S., Zhao, H., 2012. Short-term solar irradiance forecasting model based on artificial neural network using statistical feature parameters. Energies 5 (5), 1355. URL http://www.mdpi.com/1996-1073/5/5/1355 Zarzalejo, L. F., Polo, J., Mart´ın, L., Ram´ırez, L., Espinar, B., 2009. A new statistical approach for deriving global solar radiation from satellite images. Solar Energy 83 (4), 480 – 484. URL http://www.sciencedirect.com/science/article/pii/ S0038092X08002223 37 Chapter 3 Materials and methods 3.1 Introduction This chapter describes the different statistical and data mining models proposed for analyzing and forecasting solar global radiation. As mentioned in Chapter 2, new approaches based on data mining techniques have begun to be used in recent years for predicting solar radiation. The proposal of this thesis is to use both statistical and data mining techniques to these tasks. From the Statistics, the characteristics of solar goblal for a day can be analyzed using the cumulative probability distribution function. The function is used to know how the values are distributed along a day. The analysis of how many different types of day there are can also be done using the Kolmogorov-Smirnov as it can be used to test the similarity between two samples. From Data Mining, several clustering, classification and regression techniques will be tested for predicting solar radiation. The class prediction methods or classification methods allow to classify a sample in predetermined classes based on the data observed in the past. For example, it’s possible to classify the days into several types based on the observation of diffetent meteorological parameters. Classification methods presented here include Decision Trees (DT) and Support Vector Machines for classificatoin (SVC). Regression methods will be also tested. Finally some methods used to evaluate the error in predicting variables are 38 3. MATERIALS AND METHODS presented. These methods allow to evaluate and compare the performance of different modeling techniques based on the error committed. In the prediction of continuous values the error can be quantified based on the difference between the observed values and the predicted values, however, in classification models, simply determining the success rate of the technique used to classify each sample in the correct class is possible. 3.2 Statistical Methods In this section the statistical techniques used in this work are presented. The proposal is to review and analyze the possible utilization of the distribution function of values in a day and the usage of statistical tests to determine the similarity between two samples. Among those consulted in the bibliography, the proposal made in this work is to use the Kolmogorov-Smirnov two-sample test. This test allows the comparison of cumulative probability distribution functions of two different samples without the need to assume any distribution function on the underlying samples. 3.2.1 Cumulative probability distribution functions The cumulative probability distribution function (CPDF) describes the probability of a random variable Xwith a given probability distribution to have a value less than or equal to x: FX(x) = Pr(X≤x) (3.1) The right hand side represents the probability that the random variable X takes on a value less than or equal to x. Therefore, the probability that Xlies in the semi-closed interval (a, b] where a < b is: Pr(a < X ≤b) = Fx(b)−Fx(a) (3.2) Normally a capital Ffor a cumulative probability distribution function is used, while the lower-case fis used for probability distribution functions. For a continuos random variable Xthe CPDF can be expressed as the integral 3.2. STATISTICAL METHODS 39 of its probability density function Rxas follows: FX(x) = Zx −∞ fX(t)dt (3.3) 3.2.2 Kolmogorov-Smirnov two sample test The Kolmogorov-Smirnov test is based on Cumulative Probability Distribution Function (CPDF) and can be used to compare a sample with a reference probability distribution (one sample K-S test), or to compare two samples (two sample K-S test). The K-S test is a nonparametric test of the equality of continuous, one dimensional probability distributions. The K-S statistic quantifies a distance between the empirical distribution function of the sample and the cumulative distribution function of the reference distribution, or between the empirical distribution functions of two samples. Let the cumulative probability distribution function (CPDF) of Xas FX(·) and the CPDF of Yas FY(·),i.e. FX(t) = Pr(X≤t), FY(t) = Pr(Y≤t) (3.4) according to Section 3.2.1. Both FX(·) and FY(·) are assumed to be continuous. Suppose that we want to test the null hypothesis H0:FX(·) = FY(·), versus the general alternative hypothesis Ha:FX(·)6=FY(·), making no parametric assumption about the shape of these CPDF’s. The test can be performed using the Kolmogorov-Smirnov statistic that compares the empirical CPDF’s obtained with each sample. Specifically, if for any real number tit’s defined ˆ FX(t)≡n−1Pn i=1 I(Xi≤t) and ˆ FY(t)≡m−1Pm i=1 I(Yi≤t),where I(A) is the indicator function of event A, 46 3. MATERIALS AND METHODS •Supervised learning. A set of example pairs is given, (x, y), x ∈X, y ∈Yand the objetive is to find a function f:X→Yin the allowed class of functions that matches the examples. A commonly used cost is the mean-squared error, which tries to minimize the average squared error between the network’s output, f(x) and the target value yover all the example pairs. When gradient descent is used to minimize the cost for the class of neural networks called multilayer perceptrons, the well-known back-propagation algoritm is being used. •Unsupervised learning. The given cost function can be any function of the given data xand the network’s output f. The cost function depends on what is being modeled and the basic assumptions. •Reinforcement learning. Usually, the data is not given, but, is generated by the agent’s interaction with the environment. At each point in time t, the agent performs an action ytand the environment generates an observation xtand an instantaneous cost ct, according to some (usually unknown) dynamics. The aim is to discover a policy for selecting actions that minimizes some measure of a long-term cost; i.e., the expected cumulative cost. Training a neural network model means selecting one model from the set of allowed models that minimizes the cost criterion. Most of the models can be viewed as a straightforward application of optimization theory and statistical estimation. Most of the algorithms used in training artificial neural networks employ some form of gradient descent, taking the derivative of the cost function with respect to the network parameters and then changing those parameters in a gradient-related direction. The training data set consists of N training patterns (xp, tp), where p is the pattern number. xpis the input vector with dimension N and tpis the desired output vector with dimension M. ypis the network output vector for the pth pattern. The Levenberg-Marquardt (LM) method can be used for training the ANN and is one of the most popular algorithms. LM updates the neurons weight as follows: ∆w=−"µI+ P X p=1 Jp(w)TJp(w)#−1 ∇E(w) (3.10) where Jp(w) is the Jacobian matrix of the error vector Ep(w) evaluated in w, 3.3. DATA MINING TECHNIQUES 47 and I is the identity matrix. The vector error Ep(w) is the error of the network for pattern p:Ep(w) = Tp−Op(w). Some commonly used methods for training neural networks include evolutionary methods, gene expression programming, simulated annealing, expectation maximization, non-parametric methods and particle swarm optimization. 3.3.3 Support Vector Machines Support Vector Machines (SVM) is a supervised learning model (introduced by Vapnik et al. in 1992 (Berthold and Hand, 2003)), aimed to be an efficient way of learning ’good’ separating hyperplanes in a high dimensional feature space, capable of analyzing and recognizing patterns, and can be used for classification and non-linear regression analysis. A SVM is a linear machine based on Structural Risk minimization (SRM) method and on statistical learning theory and can provide good performance in pattern recognition problems without problem domain knowledge (Gorunescu, 2011). In the simplest case, given a set of training samples, each belonging to one of two categories, an SVM builds a model that assigns new samples to one of the two categories. The fundamental concept of a SVM, based on Cover’s Theorem (Cover, 1965), is that, while samples in the input (low-dimensional) space could not be linearly separable, mapping them into a sufficiently high-dimensional space (feature space) using a nonlinear mapping function ϕ(x) could more probably be linearly separable. A SVM is based on the implementation of the following two steps: •Mapping the training points by a nonlinear function ϕ(x) to a high-dimensional space in which the training points are linearly separable. •Determining the optimal separation hyperplane that maximizes margin (the distance between points of each category). In this work, SVM implementation from Chih-Chung Chang and Chih-Jen Lin is used (Chang and Lin, 2011), which is available at http://www.csie.ntu.edu.tw/ cjlin/libsvm/. This implementation supports the following types of tasks: 48 3. MATERIALS AND METHODS •SVC: Support Vector Classification, two-class and multi-class. •SVR: Support Vector Regression. •One-class SVM Four types of kernels are supported by LibSVM: linear, polynomial, radial basis function and sigmoid function. Using LIBSVM requires two steps: first, training a data set to construct a model and second, using the model to predict information from a testing data set. Getting output probabily estimates is also posible with this library. In the following section details abour ν-SVC and ν-SVR are provided as these are the variants of LibSVM that are used in this thesis. ν-Support Vector Classification Given training vectors xi∈Rn, i = 1..., l in two clases, and given y∈Rl, such that yi∈1,−1, the primal optimization problem is: min ω,b,ξ,ρ 1 2ωTω−νρ +1 l l X i=1 ξi subject to yi(ωTφ(xi) + b)≥ρ−ξi, ξi≥0, i = 1, . . . , l, ρ ≥0. (3.11) The dual problem is min α 1 2αTQα subject to 0 ≤αi≤1/l, i = 1, . . . , l, eTα≥ν, yTα= 0, (3.12) where Qi,j =yiyjK(xi, xj). (Chang and b. Lin, 2001) showed that problem 3.12 is feasible if and only if ν≤2min(#yi= +1,#yi=−1) l≤1,(3.13) 3.3. DATA MINING TECHNIQUES 49 so the range for νis smaller than (0,1]. ν-Support Vector Regression For ν−SV R also a parameter ν∈(0,1] is used, like in ν−SV C, to control the number of support vectors. Also, the parameter −SV R becomes a parameter here. The problem to solve is: min ω,b,ξ,ξ∗, 1 2ωTω+C(ν +1 l l X i=1 (ξi+ξ∗ i)) subject to (ωTφ(xi) + b)−zi≤+ξi, zi−(ωTφ(xi) + b)≤+ξ∗ i, ξi, ξ∗ i≥0, i = 1, . . . , l,  ≥0. (3.14) The dual problem is: min α,α∗ 1 2(α−α∗)TQ(α−α∗) + zT(α−α∗) subject to eT(α−α∗)=0, eT(α+α∗)≤Cν, 0≤αi, α∗ i≤C/l, i = 1, . . . , l. (3.15) The aproximate function is l X i=1 (−αi+α∗ i)K(xi, x) + b(3.16) As in ν−SV C, the inequality eT(α+α∗)≤Cν can be replaced by an equailty. Because users usually choose a small value for C,C/l may be too small. Then, in LibSVM, the users specifies C=C/l, and LibSVM solves the following problem: 50 3. MATERIALS AND METHODS min α,α∗ 1 2(α−α∗)TQ(α−α∗) + zT(α−α∗) subject to eT(α−α∗) = 0 , eT(α+α∗) = Clν, 0≤αi, α∗ i≤C, i = 1, . . . , l. (3.17) In (Chang and b. Lin, 2002) was proved that −SV R with parameters (C, ) has the same solution as ν−SV R with parameters (lC, ν). 3.3.4 Decision Trees According to (Rokach and Maimon, 2008), in machine learning environment prediction methods are commonly referred to as supervised learning, which stands opposed to unsupervised learning (which refers to modeling the distribution of instances in a input space). Supervised methods are those that attempt to establish a relationship between input attributes (also called independent variables) and a target attribute (also called dependen variable). The discovered relationship is represented in a structure called Model. The Model can be used to predict values of the dependent variable knowing the values of independent variables. Supervised models can be divided into two subclasses: Classification Models and Regression Models. Regression Models can predict a real (continuous) value while Classification Models can map the input values into predefined classes. A decision tree structure can be described as a labeled directed acyclic graph where all nodes except the root have a single parent and can have zero, one or more children. Nodes with no children are called leaf nodes (also known as “terminal” or “decision” nodes). In Figure 3.2 a example of decision tree is shown. A decision tree is a predictive model that uses a tree like structure to predict values from observations. It is useful in the areas of statistics, data mining and machine learning. The predicted values can be a finite set of values (classification trees) or can be continuous values (regression trees). In a decision tree structure each internal node denotes a test of an attribute, each branch represents an outcome of the test and each leaf node holds a class label and the topmost node in a tree is the root node (Han and Kamber, 2000). In the simplest case, each test considers a single attribute and the instance space is partitioned according to the attribute values. In the case of numeric attributes, the conditions refers to a range. 3.3. DATA MINING TECHNIQUES 51 Figure 3.2: Decision Tree Classification trees can be used to clasify objects into a predefined set of predefined classes based on their attributes. Classification trees are frequently used in applied fields such as finance, marketing, engineering and medicine. Given a tuple X, for which the associated class label is unknown, the attribute values of the tuple are tested against the decision tree, tracing a path from the root to a leaf node which holds the class predicted for that tuple. Decision tree construction does not require any domain knowledge, can handle multidimensional data and the representation of acquired knowledge in tree form is intuitive and easy for humans to assimilate. The learning process is simple and fast, and, in general, the accuracy is good, however, success depends on the data used. Tree size or complexity can have a important effect on its accuracy. Normally ,tree complexity is measured by one af these metrics: total number of nodes, total number of leaves, tree depth and number of attrinutes used. Attributes are tipically one of two types: nominal (values of an unordered set) and numeric (real numbers or so). The domain of an attribute can be denoted as dom(ai) = vi,1, vi,2, ...vi,|dom(ai)|where |dom(ai)|is the attribute cardinality. For the target attribute the domain can be represented as dom(y) = c1, c2, ...c|dom(y)|. The instance space can be defined as the Cartesian product of all the input attributes domains: X=dom(a1)×dom(a2)×... ×dom(an); and the universal instance space also includes the target attribute: U=X×dom(y). 52 3. MATERIALS AND METHODS The training set is a set of unordered tuples, which can be represented as S(B) = (hx1, y1i,hx2, y2i, ... hxm, ymi) where xq∈Xand yq∈dom(y). An induction algorithm (or inducer) is an entity capable of construct a model (classifier) that generalizes the relationship between the input attributes and the target attribute. Using this classifier it is possible to prediuct the target value of a tuple xq. The classifier can be used to classify an unseen tuple either by explicitly assigning it to a certain class or by providing a vector of probabilities representing the conditional probability of the given instance to belong to each class. Induction of an optimal decision tree from a given data set is considered to be a difficult task (NP-hard o NP-complete, Rokach and Maimon (2008)), so it is only feasible in small problems. Heuristic methods are necessary in order to construct a optimal decision tree. Methods for this can be classified into two categories, top-down and bottom-up. Examples of top-down (the most frequent) are ID3, C4.5 and CART. Split criterion is used to choose which attribute is used to create child nodes from a given node in a decision tree. Next, split cliterion is explained. Given a set of experiencies Ewe define Mas a measure of disorder in E. If we classify by an attribute Athen achild nodes will be produced, with partial measures Mj (j= 1, . . . , a). With these we calculate the new measure in the parent node, using the ponderate measure: M(A) = a X j=1 prob(AjMj(3.18) where prob(Aj) = number of experiences for jvalue for attribute Aat this node number of experiences at this node (3.19) When classifying with attribute A, the disorder decrement would be: ∆(A) = M−M(A) (3.20) The attribute with the bigger decrement will be the choosen one by the Split Criterion. Next, the complete algorithm for constructing a decision tree is presented. 3.4. METRICS 53 TreeGrowing(S, A, y, SplitCriterion, StoppingCriterion) Where: STraining Set AInput Feature Set yTarget Feature SplitCriterion - the method for evaluating a certain split StoppingCriterion - the criteria to stop the growing process Create a single tree Twith a single root node IF StoppingCriterion(S) THEN Mark Tas a leaf with the most common value of yin Sas label ELSE ∀ai∈Afind athat obtain the best SplitCriterion(ai, S) Label twith a FOR EACH outcome viof a: SET Subtreei=TreeGrowing(σa=viS, A, y) Connect the root node of tTto Subtreeiwith an edge that is labelled as vi END FOR END IF RETURN TreePrining(S, T, y) TreePruning(S, T, y) Where: STraining Set yTarget Feature TThe tree to be pruned DO Select a node tin Tsuch that pruning it maximally improve some evaluation criteria IF t 6=THEN T=prunet(T, t) UNTIL t= RETURN T Algorithm 2: Decision Tree Growing Algorithm 3.4 Metrics Model performance is usually evaluated against a single dataset and the strategy is to select the model with the best performance for this particular dataset. An 54 3. MATERIALS AND METHODS individual error is the difference between a predicted value and the corresponding observed or true value. Next, statistical methods used to measure forecasting errors are introduced. Each individual statistic used to measure the error has advantages, also disadvantages over the others. Several authors have proposed a set of metrics to evaluate general forecasting methods evaluation and error metrics, see for example (Chen and Yang, 2004). In (Voyant et al., 2015) 20 statistical parameters to estimate the short term predictability of the global horizontal irradiation time series are reviewed. The mean absolute log-return (which the author claims to have never been used before in global radiation forecasting) proves to be very efficient. In (B.Viorel, 2008) MBE (Mean Bias Error), RMSE (Root Mean Square Error) and MABE (Mean Absolute Bias Error) are cited as basic error metrics. MBE measures the systematic errors (or bias) while RMSE is mostly a measure of random errors and MABE is less frequently used than the two other statistics. In this thesis several statistics are used so that the joint vision allows provide better assess of the accuracy of prediction models. First step for calculating forecasting errors is to define the residue series, that is, the difference between observed and forecasted values: et=ˆ Zt−Zt(3.21) where etare the residue, ˆ Ztare the forecasted values of the interest variable and Ztare the actual values of the forecasted variable. 3.4.1 MAPE One of the most used method is the Mean Absolute Percentage Errror (MAPE), introduced in equation 3.22. To avoid cancelation from distinct signed values, absolute value is used. Multiplying by 100, this method’s result is expressed as percentage: 3.4. METRICS 55 MAPE =1 n n X t=1 |et Zt |(x100) (3.22) The Advantage of this method is that it allows to directly compare error between different systems because is normalized. The disavantage is that it does not discriminate error from low values, for example, low solar radiation energy values that occur early in the morning or at the sunset when energy is low. 3.4.2 MSE The Mean Square Error (MSE), also used to evaluate the output of data mining models, is estimated according to equation 3.23: MSE =Pm t=1(Xt−ˆ Xt)2 m.(3.23) 3.4.3 RMSE The next method introduced here is Root Mean Square Error (RMSE), which is described in equiation 3.24. It provides a good measure of the model’s precission because it is a cuadratic magnitude, that is, gives more importance to bigger error values. RMSE =v u u t 1 n n X t=1 e2 t(3.24) The main disadvantage of this method is that is expressed in same units as the forecasted variable so is dependent of it’s dynamic margin. This method is suitable when forecasting solar radiation only because it’s dynamic margin is similar in all places. Also, this method is a uncertainty measure widely used in energy and PV performance studies, so it allows comparision with other studies. 62 3. MATERIALS AND METHODS •16 day / daily forecast: for any location, forecast data for 16 days ahead, daily basis, available in JSON, XML or HTML format. •Historical data: historical weather data for more than 20.000 cities. •UV index: current UV index and historical data for any geo location. •Weather map layers: maps with precipitation, clouds, pressure, temperature and wind info. •Weather stations: recent data from more than 40.000 stations around the world, search weather stations close to a geographic location. •Bulk downloading: bulk files with current weather and forecasts from more than 20.000 cities. Figure 3.7 shows the main page of the OpenWeatherMap web site. Figure 3.8 shows the OpemWeatherMap platform conception, structured in three layers, with de data sources in the lower layer and the providing data services in the upper layer. Data are collected from thousands of sensors distributed all around the planet, including data from satellites and radars. Archives content data from MODIS, Landsat 7 and 8. Ground sensor network covers more than 40.000 weather stations. High-end technologies are used to automatically process billions of data points every second. Big data and cloud technologies enables OpenWeatherMap to work with immense data and process them right away providing the user with images, maps and data. For the purposes of this work, current weather data and 5 day / 3 hour forecast data are downloaded and stored in a data base. Current weather data Current weather data can be obtained using an API directly on the OpenWeatherMap website (www.openweathermap.org). The url for the API call look like this: http :// api . openweathermap . org / data /2.5/ weather ? id =2514256& ,→appid=b1b15e88fa797225412429c1c50c122a 3.5. AVAILABLE METEOROLOGICAL DATA 63 Figure 3.7: OpenWeatherMap web site The url includes two key parameters, and API KEY for user identification purposes (example: b1b15e88fa797225412429c1c50c122a) and a station id which designates the station we want the data from, in this case 2514256. The response is in JSON format and looks like this: {” coord ”:{” l on ”: −4 .4 2 ,” l a t ” : 3 6 . 7 2 }, ” weather ” : [ {” id ” :8 00 , ” main ” :” Cl ear ” ,” d e s c r i p t i o n ” : ” c l e a r ,→sky ” , ” i c o n ” : ” 01 n ” }] , ” base ” : ” cmc s t a t i o n s ” ,” main ” :{” temp ” : 2 8 4 . 4 5 3 , ” p r e s s u r e ” : 1 0 1 8 . 2 4 , ” ,→humidit y ” :1 0 0 , ” temp min ” : 2 8 4 . 4 5 3 , ” temp max ” : 2 8 4 . 4 5 3 , ” s e a l e v e l ” : 1 0 3 4 . 6 6 , ” g r n d l e v e l ,→”:1018.24}, ” wind ” : {” sp ee d ” : 2 . 8 6 , ” deg ” : 1 1 6 . 5 0 5 }, ” c l o u d s ” : {” a l l ” : 0 }, ” dt ” : 14 5 78 2 95 3 1 , ” s y s ,→”: {” message ” : 0 . 0 0 2 8 , ” co un try ” : ”ES” , ” s u n r i s e ” :1 45 78 50 66 1 , ” s u ns e t ” :1 457 89 34 13 },” i d ,→”: 2 51 4 25 6 , ” name ” : ” Malaga ” ,” cod ” : 20 0 } This responde includes the following parameters: temperature, pressure, humidity, wind speed and direction and the name of the station, in this case M´alaga. 64 3. MATERIALS AND METHODS Figure 3.8: OpenWeatherMap platform diagram 5 day / 3 hour forecast Similar to current weather, 5 day forecasts can be downloaded using an API on the web site. The url look like this: http :// api . openweathermap . org /data /2.5/ f o r e c a s t ? id =2514256& ,→appid=b1b15e88fa797225412429c1c50c122a Again, the url includes the two key parameters, the API KEY of the user and the station code. The response is again in JSON format: {” c i t y ” :{” i d ” : 2 51 4 25 6 , ” name ” : ” Malaga ” ,” co or d ” : {” l on ”: −4 .4 20 34 ,” l a t ” : 3 6 . 7 2 0 1 6 1 }, ” co un tr y ” : ” ES ,→” , ” p o pu l at i on ” : 0 , ” s ys ” : {” po pul at i on ” : 0} } , ” cod ”: ” 20 0 ” , ” message ” : 0 . 0 0 2 6 , ” cnt ” : 4 0 , ” l i s t ,→” : [ {” dt ”: 14 57 83 80 00 , ” main ”: {” temp ” : 2 8 4 . 4 5 , ” temp min ” : 2 8 4 . 4 5 , ” temp max ” : 2 8 4 . 4 5 2 , ” ,→p r e ss u r e ” : 1 0 1 7 . 3 , ” s e a l e v e l ” : 1 0 3 3 . 7 2 , ” g r n d l e v e l ” : 1 0 1 7 . 3 , ” hu midity ” : 10 0 , ” temp k f ” : 0}, ” ,→weather ”:[{” i d ” :8 00 , ” main ” :” Clear ” ,” d e s c r i p t i o n ” :” c l e a r sky ” ,” i co n ”: ”0 1 n ”}] , ” c l o u d s ” : {” ,→a l l ” : 0 }, ” wind ” : {” s peed ” : 3 . 3 1 , ” deg ” : 8 5 . 0 0 4 5 }, ” s y s ” : {” pod ” : ” n ”},” d t t x t ”:”2016 −03−13 ,→03:00:00”},{” dt ”:145784 88 00 , ” main ”:{” temp ” : 2 8 3 . 9 4 , ” temp min ” : 2 8 3 . 9 4 , ” temp max ,→” : 2 8 3 . 9 4 4 , ” p r e s s u r e ” : 1 0 1 7 . 3 4 , ” s e a l e v e l ” : 1 0 3 3 . 8 , ” g r n d l e v e l ” : 1 0 1 7 . 3 4 , ” hu midit y ” : 10 0 , ” ,→te mp kf ” : 0 }, ” w ea ther ” : [ {” id ” :8 00 , ” main ” :” Cl ear ” ,” d e s c r i p t i o n ” : ” c l e a r sky ” , ” i co n ”:” 01 n ,→”}] , ” c l o u d s ” : {” a l l ” : 0 }, ” wind ” :{” s peed ” : 3 . 8 6 , ” deg ” : 6 5 . 5 0 3 7 }, ” s y s ” : {” pod ” : ” n ”},” d t t x t ,→”:”2016 −03 −13 06:00:00”},{” dt ”: 14 57 85 96 00 ,” main ” :{” temp ” : 2 8 5 . 6 7 , ” temp min ” : 2 8 5 . 6 7 , ” ,→temp max ” : 2 8 5 . 6 7 1 , ” p r e s s u r e ” : 1 0 1 8 . 2 8 , ” s e a l e v e l ” : 1 0 3 4 . 6 2 , ” g r n d l e v e l ” : 1 0 1 8 . 2 8 , ” hum idity ,→” : 1 0 0 , ” temp k f ” : 0 }, ” weath er ” : [ {” i d ” :8 00 , ” main ” :” C lear ” ,” d e s c r i p t i o n ” :” c l e a r sky ” ,” ico n ,→”:”01d”}] , ” c l o u d s ” : {” a l l ” : 0 }, ” wind ” :{” s peed ” : 4 . 0 6 , ” deg ” : 6 3 . 5 0 5 3 }, ” s y s ” : {” pod ” : ” d ”},” ,→d t t x t ”:”2016 −03−13 09:00:00”},{” dt ” :1 45 78 70400 , ” main ” :{” temp ” : 2 8 7 . 1 5 , ” temp min ,→” : 2 8 7 . 1 5 , ” temp max ” : 2 8 7 . 1 5 1 , ” p r e s s u r e ” : 1 0 1 8 . 1 1 , ” s e a l e v e l ” : 1 0 3 4 . 2 3 , ” g r n d l e v e l ,→” : 1 0 1 8 . 1 1 , ” humidity ” : 9 4 , ” t emp k f ” : 0 }, ” weat he r ” : [ {” id ” :8 00 , ” main ” :” C le ar ” ,” d e s c r i p t i o n ,→”: ” c l e a r sky ” ,” ico n ” :” 01 d ” }] , ” c l o u d s ” : {” a l l ” : 0 }, ” wind ” : {” s pe ed ” : 2 . 6 2 , ” deg ” : 1 2 0 . 5 1 5 }, ” ,→sy s ”: {” pod ” : ” d ”},” d t t x t ”:”2016 −03−13 12:00:00”},{” dt ” :1 457881200 , ” main ” :{” temp ,→” : 2 8 6 . 9 7 , ” temp min ” : 2 8 6 . 9 6 7 , ” temp max ” : 2 8 6 . 9 7 , ” p r e s s u r e ” : 1 0 1 6 . 3 8 , ” s e a l e v e l ” : 1 0 3 2 . 3 3 , ” ,→g r n d l e v e l ” : 1 0 1 6 . 3 8 , ” humi dity ” : 9 5 , ” t emp k f ” : 0 }, ” wea th er ” : [ {” i d ” : 8 0 0 , ” main ” : ” C l ea r ” , ” ,→d e s c r i p t i o n ” :” c l e a r sky ” ,” i co n ” :” 01 d ” }] , ” c l o u d s ” : {” a l l ” : 0 }, ” wind ” : {” speed ” : 2 . 4 1 , ” deg ,→”:149.001}, ” s y s ” : {” pod ” : ” d ”},” d t t x t ”:”2016 −03−13 15:00:00”},{” dt ”:145789 20 00 , ” main ”:{” ,→temp ” : 2 8 6 . 1 2 , ” temp min ” : 2 8 6 . 1 2 , ” temp max ” : 2 8 6 . 1 2 4 , ” p r e s s u r e ” : 1 0 1 5 . 4 1 , ” s e a l e v e l 3.5. AVAILABLE METEOROLOGICAL DATA 65 ,→” : 1 0 3 1 . 4 4 , ” g r n d l e v e l ” : 1 0 1 5 . 4 1 , ” h umidit y ” :1 0 0 , ” temp k f ” : 0 }, ” weathe r ” : [ {” i d ” : 8 0 0 , ” main ,→”: ” Cl ear ” ,” d e s c r i p t i o n ” :” c l e a r sky ” ,” i co n ”:” 01 d ” }] , ” c l o u d s ” : {” a l l ” : 0 }, ” wind ” : {” speed ,→” : 1 . 0 6 , ” deg ” :1 57 },” s ys ” :{” pod ” : ” d ”},” d t t x t ”:”2016 −03−13 18:00:00”},{” dt ” :1 45 79 02 80 0 , ” ,→main ” :{” temp ” : 2 8 4 . 5 7 , ” temp min ” : 2 8 4 . 5 7 , ” temp max ” : 2 8 4 . 5 7 2 , ” p r e s s u r e ” : 1 0 1 5 . 6 4 , ” s e a l e v e l ,→” : 1 0 3 1 . 9 , ” g r n d l e v e l ” : 1 0 1 5 . 6 4 , ” humi di ty ” : 1 0 0 , ” tem p k f ” : 0 }, ” weather ” : [ {” i d ” : 8 0 0 , ” main ” : ” ,→Cle ar ” ,” d e s c r i p t i o n ” : ” c l e a r sky ” ,” i co n ” :” 01 n ” }] , ” c l o u d s ” : {” a l l ” : 0 }, ” wind ” :{” speed ,→” : 0 . 9 2 , ” deg ” : 1 4 6 . 0 0 4 }, ” s y s ” : {” pod ” : ” n ”},” d t t x t ”:”2016 −03−13 21:00:00”},{” dt ,→”:1 45 79 13 60 0 ,” main ” :{” temp ” : 2 8 4 . 0 2 , ” temp min ” : 2 8 4 . 0 2 , ” temp max ” : 2 8 4 .0 2 1 , ” p r e s s u r e ,→” : 1 0 1 5 . 4 6 , ” s e a l e v e l ” : 1 0 3 1 . 7 4 , ” g r n d l e v e l ” : 1 0 1 5 . 4 6 , ” hum idity ” : 1 00 , ” t emp k f ” : 0 }, \vdots ”weather ”:[{” i d ” :8 01 , ” main ” :” Clouds ” , ” d e s c r i p t i o n ” :” few c lo ud s ” ,” i co n ” :” 02 n ” }] , ” c l o u d s ” : {” a l l ,→” : 2 4 }, ” wind ” : {” s pe ed ” : 4 . 7 1 , ” deg ” : 9 1 . 0 0 1 7 }, ” s y s ” : {” pod ” : ” n ”},” d t t x t ”:”2016 −03−18 ,→00:00:00”}]} OpenWeatherMap database A simple dababase is used to store OpenWeatherMap’s current and forecasted values. An application has been developed to automatically read data from web site and and store it into the database. 3.5.3 Spanish weather service The Agencia Espa˜nola de Meteorolog´ıa (AEMET) is the official weather service of Government of Spain, (AEMET, 2016). The service provides data about observations as well as forecasting and also historical data can be obtained. Alerts from extreme meteorological conditions (wind, rain, snow) can be retrieved from the web site, represented over a map of Spain terrirory. In the Aemet website there are forecasting data available to download in XML format. These data include forecast for 5 days and include these parameters: temperature, humidity, rain probability and sky cover index. The forecasting values are delivered for 6 hours periods. Data are available to download in XML format via an API, example of forecasting data for Malaga city can be found here: http://www.aemet.es/xml/municipios/localidad_29067.xml. Each meteorological station has a numeric code for identification purposes, this code is used in the URL for retrieveing data. The code for the station used here is 29067. Here is a portion of the XML response for the previous url, this XML response is related to M´alaga city and to the next day weather: <d i a f e c h a=”2016−01−08”> 66 3. MATERIALS AND METHODS Figure 3.9: AEMET Website <p r o b p r e c i p i t a c i o n pe ri od o=”00−24”>85</ p r o b p r e c i p i t a c i o n> <p r o b p r e c i p i t a c i o n pe ri od o=”00−12”>0</ p r o b p r e c i p i t a c i o n> <p r o b p r e c i p i t a c i o n pe ri od o=”12−24”>85</ p r o b p r e c i p i t a c i o n> <p r o b p r e c i p i t a c i o n pe ri od o=”00−06”>0</ p r o b p r e c i p i t a c i o n> <p r o b p r e c i p i t a c i o n pe ri od o=”06−12”>0</ p r o b p r e c i p i t a c i o n> <p r o b p r e c i p i t a c i o n pe ri od o=”12−18”>0</ p r o b p r e c i p i t a c i o n> <p r o b p r e c i p i t a c i o n pe ri od o=”18−24”>80</ p r o b p r e c i p i t a c i o n> <c o t a n i e v e p r o v p e ri od o=”00−24”/> <c o t a n i e v e p r o v p e ri od o=”00−12”/> <c o t a n i e v e p r o v p e ri od o=”12−24”/> <c o t a n i e v e p r o v p e ri od o=”00−06”/> <c o t a n i e v e p r o v p e ri od o=”06−12”/> <c o t a n i e v e p r o v p e ri od o=”12−18”/> <c o t a n i e v e p r o v p e ri od o=”18−24”/> <e s t a d o c i e l o pe ri od o=”00−24” d e s c r i p c i o n=” I n t e r v a l o s nubosos con l l u v i a e sc a sa ”>43</ ,→estado cielo> <e s t a d o c i e l o pe ri od o=”00−12” d e s c r i p c i o n=” Despejado ”>11</ e s t a d o c i e l o> <e s t a d o c i e l o pe ri od o=”12−24” d e s c r i p c i o n=” I n t e r v a l o s nubosos con l l u v i a e sc a sa ”>43</ ,→estado cielo> <e s t a d o c i e l o pe ri od o=”00−06” d e s c r i p c i o n=” Despejado ”>11n</ e s t a d o c i e l o> <e s t a d o c i e l o pe ri od o=”06−12” d e s c r i p c i o n=” Despejado ”>11</ e s t a d o c i e l o> <e s t a d o c i e l o pe ri od o=”12−18” d e s c r i p c i o n=”Poco nuboso ”>12</ e s t a d o c i e l o> <e s t a d o c i e l o pe ri od o=”18−24” d e s c r i p c i o n=”Nuboso con l l u v i a e sca s a ”>44n</ e s t a d o c i e l o> <v i e n t o p e r i od o=”00−24”> <direccion>SE</ d i r e c c i o n> <velocidad>15</ v e l o c i d a d> </ v i e n t o> <v i e n t o p e r i od o=”00−12”> <direccion>O</ d i r e c c i o n> <velocidad>15</ v e l o c i d a d> </ v i e n t o> <v i e n t o p e r i od o=”12−24”> <direccion>SE</ d i r e c c i o n> <velocidad>15</ v e l o c i d a d> </ v i e n t o> <v i e n t o p e r i od o=”00−06”> <direccion>C</ d i r e c c i o n> <velocidad>0</ v e l o c i d a d> </ v i e n t o> <v i e n t o p e r i od o=”06−12”> <direccion>SE</ d i r e c c i o n> <velocidad>15</ v e l o c i d a d> </ v i e n t o> <v i e n t o p e r i od o=”12−18”> <direccion>SE</ d i r e c c i o n> <velocidad>10</ v e l o c i d a d> </ v i e n t o> <v i e n t o p e r i od o=”18−24”> <direccion>S</ d i r e c c i o n> <velocidad>10</ v e l o c i d a d> </ v i e n t o> <racha max p eri odo=”00−24”/> <racha max p eri odo=”00−12”/> <racha max p eri odo=”12−24”/> 3.5. AVAILABLE METEOROLOGICAL DATA 67 <racha max p eri odo=”00−06”/> <racha max p eri odo=”06−12”/> <racha max p eri odo=”12−18”/> <racha max p eri odo=”18−24”/> <temperatura> <maxima>20</maxima> <minima>11</minima> <dato hora=” 06”>11</ dato> <dato hora=” 12”>19</ dato> <dato hora=” 18”>17</ dato> <dato hora=” 24”>16</ dato> </temperatura> <sens termica> <maxima>20</maxima> <minima>11</minima> <dato hora=” 06”>11</ dato> <dato hora=” 12”>19</ dato> <dato hora=” 18”>17</ dato> <dato hora=” 24”>16</ dato> </sens termica> <hu me da d re la tiva> <maxima>90</maxima> <minima>60</minima> <dato hora=” 06”>80</ dato> <dato hora=” 12”>60</ dato> <dato hora=” 18”>75</ dato> <dato hora=” 24”>85</ dato> </ humed ad relat iv a> <uv max>2</uv max> </ dia> For this study only next day forecasted data are used because the desired forecasting horizon is 24 hours. For the next day these are the meteorological data parameters available from AEMET: •Rain probability: is an estimation of the probability that it will rain in that period of time (%). •Snow line: the altitude above which snow and ice can be found (meters). •Sky cover: Sky cover indicates the amount of clouds in the sky. •Wind speed and direction: direction (N,S,E,W) and speed (km/h). •Air temperature: is the ambient air temperature expressed in °C. •Humidity: is the usual relative air humidity (%). These data are provided in a 6 hours period basis (0 to 6, 6 to 12, 12 to 18 and 18 to 24 h.). A data base to store these data has been developed because past data are not avaiable from AEMET and these data are needed to train any algorithm or system to forecast some variables. The data base has been developed over PostgreSQL 9.0. 68 3. MATERIALS AND METHODS These data is read from AEMET and stored in the data base every day, although not all fields are stored, but only the interesting ones: air temperature, humidity, rain probability and sky cover. Available sky cover levels are shown in table 3.1. Code Description 11 Clear 12 Little Cloud 13 Cloudy intervals 14 Cloud 15 Cloudy 16 Covered 17 High clouds 23 Cloudy intervals with rain 24 Cloud with rain 25 Cloudy with rain 26 Covered with rain 35 Cloudy with snow 43 Cloudy intervals and low rainfall 44 Cloud with a little rain 45 Cloudy with a little rain 46 Covered with little rain 51 Cloudy intervals with storm 52 Cloud with storm 53 Cloudy with storm 54 Covered with storm 61 Cloud intervals with storm and low rainfall 62 Cloud with storm and low rainfall 63 Cloudy woth storm and low rainfall 64 Covered with storm and low rainfall 72 Cloud with little snow 73 Cloudy with little snow 74 Covered with little snow Table 3.1: Sky cover codes and descriptions from AEMET 3.6. CONCLUSIONS 69 3.6 Conclusions In this chapter statistical and data mining models that will be used in characterizacion and modeling the hourly solar global radiation and the different metrics to evaluate the performance of proposed models are presented. Moreover, the different sources of data used in this work are described. The data used are both historical record from weather public services as the AEMET (Agencia Espa˜nola de Meteorolog´ıa) and predictions. In the case of forecasted data, the API developed for obtaining these data are presented. All these software will be included in the prediction of solar radation tool developed in the frame of this thesis. 71 Bibliography AEMET, 2016. Agencia espa˜nola de meteorolog´ıa (web site). http://www.aemet. es, accessed: 2015-09-30. Berthold, M. R., Hand, D. J. (Eds.), 2003. Intelligent Data Analysis: An Introduction, 2nd Edition. Springer Verlag. B.Viorel, 2008. Modeling Solar Radiation at the Earths Surface. Recent Advances. Springer. Celebi, M. E., Kingravi, H. A., Vela, P. A., 2013. A comparative study of efficient initialization methods for the k-means clustering algorithm. Expert Systems with Applications 40 (1), 200 – 210. URL http://www.sciencedirect.com/science/article/pii/ S0957417412008767 Chang, C. C., b. Lin, C., Sept 2001. Training v-support vector classifiers: Theory and algorithms. Neural Computation 13 (9), 2119–2147. Chang, C. C., b. Lin, C., Aug 2002. Training v-support vector regression: Theory and algorithms. Neural Computation 14 (8), 1959–1977. Chang, C.-C., Lin, C.-J., 2011. LIBSVM: A library for support vector machines. ACM Transactions on Intelligent Systems and Technology 2, 27:1–27:27, software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm. Chen, Z., Yang, Y., 2004. Assessing forecast accuracy measures. Gorunescu, F., 2011. Data Mining: Concepts, Models and Techniques. Intelligent Systems Reference Library. Springer Berlin Heidelberg. Han, J., Kamber, M., 2000. Data mining: Concepts and techniques. Jain, A., Murty, M., Flynn, P., 1999. Data clustering: A review. ACM Computing Surveys 31(3), 264–323. 78 4. FUNDAMENTALS OF SOLAR GLOBAL RADIATION Figure 4.2: Azimuth angle Ψ, solar elevation eand zenith angle θ 4.3 Extraterrestrial solar radiation The extraterrestrial radiation is the radiation received outside the atmosphere or the radiation at the surface if there were not atmosphere. It can be calculated using the solar constant and the exact distance between Sun and Earth. The solar constant is defined as the amount of total solar energy (including all wavelengths) received by unit time and unit area at the mean Sun-Earth distance. The solar constant was considered effectively a constant but nowadays is recognised that this term is varying continuosly, from minutes to years or even decades, but particularly within the 11-years solar activity cycle. Therefore, this magnitude is referred to as total solar irradiance (TSI) and the term solar constant can be used to describe the long-term average of TSI. For Isc, (Gueymard, 2004) gives a value of 1366.1 Wm−2, years before (Iqbal, 1983) gave a value of 1367 Wm−2 4.4. HOURLY SERIES OF GLOBAL SOLAR RADIATION 79 The extraterrestrial radiation can be calculated using the folloging formulas. The hourly extraterrestrial radiation, G0,h, received on an horizontal surface is obtained using equation 4.6: G0,h =IscE0cos θz=IscE0(sin δsin φ+ cos δcos φcos ωs)(Whm−2),(4.6) where Isc is the solar constant, E0is the eccentricity factor, δis the declination angle, φis the latitude, ωsis the hour angle centered at (ωs−π/24, ωs+π/24). The expressions to estimate E0,δand ωscan be found in (Iqbal, 1983). The daily extraterrestrial radiation, G0,d, received on an horizontal surface is obtained using equation 4.7: G0,d =24 πIscE0(ωsr sin δsin φ+ cos δcos φsin ωsr)(Whm−2),(4.7) where ωsr is sunrise angle. 4.4 Hourly series of global solar radiation Hourly series are used as base to study and analyze solar radiation characteristics and behaviour and as input to the machine learning algorithms in order to forecast solar global radiation. Data recorded by meteorological stations can be recorded at different periods of time, for exaple, every minute, every 10 minutes, every hour, etc. In this case, available data are recorded usually every minute so previous processing is necessary in order to obtain hourly series from it. 4.4.1 Analysis of series It’s desirable that data series (instantaneous or hourly) show no trends (seasonal or daily) so the prediction model can concentrate on the stochastic behaviour 80 4. FUNDAMENTALS OF SOLAR GLOBAL RADIATION Hours 12345678 W/m2 100 200 300 400 500 600 700 800 900 1000 1100 Figure 4.3: Hourly global radiation for several days of year 2013 of the series. Series that avoid showing these trends would be better as input for forecasting methods because these forecasting methods would not have to predict different radiation levels in summer and winter but only manage to predict solar radiation values according to weather conditions. The hourly solar extraterrestrial and solar global radiation at noon throughout the year 2013 are shown in figure 4.4, where the trends in summer and winter can be easily noted. Figure 4.3 shows the global solar horizontal radiation for several day at the same station. It’s clear that along the different year seasons radiation levels can be very different. It’s desirable that for the series to be stationary. Thus, study series (hourly and daily) of global solar radiation can be done following these steps: •Find and characterize observed trends and eliminate them from original series. •Analyze statistical properties of the resulting series: mean values, variance, 4.4. HOURLY SERIES OF GLOBAL SOLAR RADIATION 81 Days 0 50 100 150 200 250 300 350 W/m2 0 200 400 600 800 1000 1200 1400 Extraterrestrial Global hor. Figure 4.4: Hourly global radiation at 12h and extraterrestrial hourly global radiation (left) and hourly clearness index at 12 h (right) for year 2013 distribution functions, etc. •Depending on observed properties, propose the best models to simulate original series and predict new values. Next, some techniques to eliminate daily and annual trends are reviewed. 4.4.2 Methods to remove the seasonal and daily trend in global radiation series In the literature there are several ways to eliminate seasonal and daily trends for solar global radiaton, some of them are: •Average or moving average using Fourier functions •Physical models •Maximum exposure hourly global radiation values 82 4. FUNDAMENTALS OF SOLAR GLOBAL RADIATION Average or moving average using Fourier functions In (Balouktsis and Tsalides, 1986) authors obtain Fourier coefficients for each hour of the day along the year and use this function: GF h,d =A0,h +A1,h cos 2π(d−1) L+B1,h sin 2π(d−1) L(4.8) where GF h,d is the value of the Fourier series for day dand hour h;Bi,Aiare the Fourier coefficients and Lis the considered hour period. Using these GF h,d values, an index Fis calculated as follows: Fh,d =Gh,d GF h,d (4.9) However, this series is not normal and it is necessary to use gaussin techniques. Physical models A “clear sky” model is a model capable of estimate the solar irradiance received at the Earth’s surface. It may produce estimates of direct (Ebnc), diffuse (Edhc) and global (Eghc) solar irradiance (Engerer and Mills, 2015). For example, (Ineichen and Perez, 2002) describes a model for direct and global radiation as: Ebnc =b∗Eextn∗exp(−0.09 ∗AM ∗(TL−1)) Eghc =a1∗Eextn∗cos(Θz)∗exp(−a2∗AM)∗(fh1+fh2∗(TL−1)) (4.10) These formulations include empirical adjustements for altitude coefficients and incorporate also turbidity information via Linke turbidity coefficient. (Ineichen, 2008) proposed a simplified version of the Solis model, which was computationally expensive and required sparsely measured input. The simplified 4.4. HOURLY SERIES OF GLOBAL SOLAR RADIATION 83 model is easier to use and is capable of producing estimates for direct, diffuse and global radiation: Ebnc =E0 ext ∗exp(−τb/cos(Θz)b) Eghc =E0 ext ∗exp(−τg/cos(Θz)g)∗cos(Θz) Edhc =E0 ext ∗exp(−τd/cos(Θz)d) (4.11) where E0 ext, τb, τg, τdare all dependent on the aerosol optical depth, water vapor and atmospheric pressure are also required as inputs. The Esra Model was developed for the Europan Solar Radiation Atlas (Rigollier et al., 2000). The model requires the Linke turbidity as input and uses and air mass based parametrization for Rayleigh optical thickness (δR). The model can estimate direct and diffuse solar radiation: Ebnc =Eextn∗exp(−0.8662 ∗TL∗AM ∗δR) Edhc =Eextn∗TRd(TL)∗Fd(Θz, TL)(4.12) The REST2 model (Gueymard, 2008) is separated into two bands representing the broadband components of two separate series of spectra and incorporates transmission estimates for Rayleigh scattering (TRi), uniform gas (TGi), ozone (TOi), nitrogen dioxide (TNi), and water vapor (TWi) absortion and aerosol extinction (TAi). The beam estimate is: Ebnci=Eextni ∗TRi∗TGi∗TOi∗TNi∗TWi∗TAi(4.13) The diffuse clear sky estimate (Edi) is divided into two layers, in the upper layer component (Edpi) Rayleigh scattering, ozone and uniform gas absortion are estimated. For the bottom layer (Eddi), aerosol, water vapor, and nitrogen dioxide absortion and, separately, aerosol scattering (Tasi), ground and sky albedo (ρGi, ρSi) are used: Edpi=Eexthi ∗TGi∗TOi∗TNi∗TWi ∗[BRi∗(1 − ∗TRi)∗T0.25 Ai+BA∗Fi∗TRi∗(1 − ∗T0.25 asi )] Eddi=ρGi∗ρSi∗(Ebni∗cos(Θz) + Edpi)/(1 −ρGi∗ρSi) (4.14) 84 4. FUNDAMENTALS OF SOLAR GLOBAL RADIATION Maximum exposure hourly global radiation values This method is based on the maximum expected values of hourly exposition and can be also used for daily exposition series. (Boileau, 1983) describes two models for forecasting solar radiation, with first, Model B, it studies the increments of daily global radiation from one day to the next and shows that the fluctuations can be modeled using a white noise ai: vi=ai−Θai−1(4.15) It was demonstrated that the same value for parameter Θ (Θ = 0.75) can be used for the three sites considered. With the second model, Model C, better results are achieved: ˜ Ii=Ii− hIii(4.16) so forecasts are calculated with one past day only using a Markov model with the centered variable ˜ Ii. 4.4.3 Clearness Index Clearness index is commonly used to analyze and characterize solar global radiation as this index allows the seasonal and daily trends observed in solar radiation to be removed. Using this index instead of clear sky index is proposed as it is possible to estimate the first one only using well known expressions without knowing the specific weather conditions of the locations, sometimes necessary for some of clear sky models, (Zhong and Kleissl, 2015). The main problem of using Clearness index is that obtained series are not truly stationary timeseries. The clearness index is defined as the ratio of the horizontal global radiation to the radiation received outside the atmosphere or extraterrestrial irradiance according to Expression 4.17: Kt=Gt Gt,0 ,(4.17) where Gtis the solar global radiation recorded for time tand Gt,0is the extraterrestrial solar global radiation for this period. 4.5. SOLAR FORECAST METRICS 85 According to (Woyte et al., 2007), the instantaneous Clearness index Ktor transmission coefficient over a arbitrarly oriented surface is: Kt=G S0Ecos I(4.18) where S0is the solar constant, Eis the eccentricity correction factor and Iis the angle of incidence. This formulation of Ktis dependent on the zenith angle. 4.5 Solar forecast metrics To quantify the accuracy of solar forecasts various metrics can be used. Determining which one is most appropiate depends on the user: while system operators need metrics that accurately reflects the cost of forecast errors, researchers require indicators of relative performance of different forecast models and indicators for a single model under different conditions. Also, an appropiate test dataset and analysis procedure are very important. Test dataset must exclude all data used to train models so that evaluation is performed on new or unseen data. Data has to be checked to remove observation errors and hardware and operation issues so that forecast is done over good quality data without forecasting deterministic data. Performance metrics can be categorized according to three types of forecasting error: bias, variance and correlation. In (Hoff et al., 2012) several absolute and relative statistical metrics for errors in forecasting are analyzed and it shows that a large number of metrics are needed to provide a clear picture of the forecasting accuracy of any method. In 3 several general error methods are introduced: MSE, RMSE,MAE, etc. Bias characterizes the balance between over and underprediction. Mean Bias Error is the most commonly used bias measure and is defined as MBE =1 N N X t=1 (I(t)−ˆ I(t)) (4.19) where I(t) is the measured irradiance at time t,ˆ I(t) is the forecasted irradiance at time tand Nis the number of data points in the data set. If the case of a 86 4. FUNDAMENTALS OF SOLAR GLOBAL RADIATION perfect forecast (ˆ I(t) = I(t)) this metric returns 0 but also for situations where positive and negative errors simply cancel out by summing to 0. The coefficient of determination R2is able to measure how well forecast values predict trends in measured values, it is a comparison between the variance of the errors and the variance of the data to be modeled: R2= 1 −σ2(ˆ I−I) σ2(I)(4.20) where σ2is the variance of the dataset. For perfect forecasting R2= 1. The value of R2is directly related to the RMSE by noticing that: R2≈1−RMSE2 var(I)(4.21) To evaluate variance of forecast errors two metrics are commonly used: the root mean square error (RMSE) and the mean absolute error (MAE). RMSE is related to the standard deviation of the errors. RMSE is calculated as: RMSE =v u u t 1 N N X t=1 (ˆ I(t)−I(t))2(4.22) Tipically night values are removed in the above calculations of R2and RMSE. It is usual to normalize values relative to energy produced or energy capacity; utilities tend to use the second one (more favorable) while scientists tend to prefer the first one. None of these metrics take into account the variability in the irradiance time-series data. (Perez et al., 2010) found that RMSE error is lower in places with sunnier (less variable) weather conditions. 4.5.1 Persistence model In general, a persistence model is a model that assumes that the next value will be equal to the current observed value. For example, a persistence model for the 4.5. SOLAR FORECAST METRICS 87 exchange rate between Euro and Dollar would predict that tomorrow the value of the exchage rate will be the same as today. Another example, a persistence model for the wind speed would say that the wind speed for the next hour will be the same speed as now. For hourly solar irradiance, a persistence model can be built taking into account that solar irradiance varies within the day hours. It makes no sense assuming that irradiance will be the same along the different day hours. A basic model can be built assuming that for a given day, the irradiance at a certain hour will be the same to the irradiace at the previous day at the same hour. For daily solar irradiance, a simpler model can be built assuming that the daily irradiance for the next day will be the same as for current day. 4.5.2 Forecast Skill The forecast skill over 24 hour persistence forecasts, s, proposed in (Coimbra and Kleissl, 2013), is estimated according to equation 4.23: s= 1 −RMSEmodel RMSEpersistence .(4.23) The forecast skill is independent of the specific meteorological or climatological characteristics of the site under consideration, (Coimbra and Kleissl, 2013). The forecast skill as defined above is such that when s= 1 the solar forecast is perfect and when s= 0 the forecast uncertainty is as large as the variability. By definition, the persistence model should have a forecast skill s= 0, because the error is equal to the variability, so, the ratio expressed in 4.23 is a measure of improvement over the persistence forecast. If the value of sis negative, then, the model performance is worse than the persistence forecast. So, forecast models with better performance than persistence model will have a value for sbetween 0 and 1, with higher values for better performance. 94 5. FORECASTING WITH CPDF 0 10 20 30 40 50 60 70 80 90 100 Probability distribution function 0 1 2 3 4 5 6 7 8 9 10 2010-11-01 0 10 20 30 40 50 60 70 80 90 100 CPDF 0 0.2 0.4 0.6 0.8 1 2010-11-01 Figure 5.4: probability distribution values (left) and CPDF (right) for day 201011-01 0 10 20 30 40 50 60 70 80 90 100 Probability distribution function 0 1 2 3 4 5 6 7 8 9 10 2011-02-03 0 10 20 30 40 50 60 70 80 90 100 CPDF 0 0.2 0.4 0.6 0.8 1 2011-02-03 Figure 5.5: Accumulated values (left) and CPDF (right) for day 2011-02-03 0 10 20 30 40 50 60 70 80 90 100 Probability distribution function 0 1 2 3 4 5 6 7 8 9 10 2011-10-02 0 10 20 30 40 50 60 70 80 90 100 CPDF 0 0.2 0.4 0.6 0.8 1 2011-10-02 Figure 5.6: Accumulated values (left) and CPDF (right) for day 2011-10-02 5.4. CLUSTERING CPDF 95 5.4 Clustering CPDF There are different Cumulative Probability Distribution Functions depending on the observed values of clearness index estimated each 20 minutes. The proposal is to obtain the minimum number of CPDF that represent all observed curves using a clustering technique that cluster them and to prove that these curves are statistically equals to the rest of estimated CPDF. To achieve this goal the use of a clustering technique is proposed. Specifically, the K-means method has been used. CPDF vectors are used as the inputs for K-means clustering. The number of clusters need to be fixed in advance in order to use K-means. We checked four different number of clusters, from 4 to 7, taking into account the expected different CPDF’s. For each execution of K-means we obtained the centroid of each cluster. Using these centroids we checked the equality of CPDF’s of the cluster in each test using the Kolmogorov-Smirnov two sample test. Table 5.1 shows the results obtained for the different numbers of checked clusters when using α= 0.05 (significance level). Number of CPDF’s for which % CPDF’s % of clusters clusters Dn,m > cαDn,m > cα 4 29 4.4 100 5 26 3.9 80 6 21 3.1 33 7 19 2.9 29 Table 5.1: Results of Kolmogorov-Smirnov two sample test for different number of clusters Taking these results into account, 6 clusters are selected instead of 7 clusters as only 2 of 6 clusters in both cases have some CPDF’s (approximately 3%) that are significantly different from theirs centroids and therefore 6 clusters are enough to capture the different CPDF’s observed. The days included in each cluster are from different months in all cases which means that the clustering method allows us to capture the different CPDF’s observed along the year and not the season of the year. Figures 5.7, 5.8 and 5.9 shows CPDF’s distribution for each cluster when using K= 6. 96 5. FORECASTING WITH CPDF 0 20 40 60 80 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kt(x100) CPDF 0 20 40 60 80 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kt(x100) CPDF Figure 5.7: CPDF’s for clusters 1 (left) and 2 (right) 0 20 40 60 80 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kt(x100) CPDF 0 20 40 60 80 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kt(x100) CPDF Figure 5.8: CPDF’s for clusters 3 (left) and 4 (right) 5.5 Relation between Kdvalues and clusters After clustering all the observations, the relationship between the daily clearness index corresponding to each CPDF and its cluster is analyzed. Figure 5.10 shows these values. As can be observed, the daily clearness index is related to the cluster to which the CPDF for the given day belongs. These results can be used to decide the cluster to which one day belongs. It has been checked that for CPDF’s that can belong to two different clusters (see Fig.5.10) the Dn,m between these each CPDF and the centroid of each possible cluster is always lower than the critical value cα. Assigning the cluster depending on the daily clearness index value using the expression 5.1 is proposed. 5.5. RELATION BETWEEN KDVALUES AND CLUSTERS 97 0 20 40 60 80 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kt(x100) CPDF 0 20 40 60 80 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kt(x100) CPDF Figure 5.9: CPDF’s for clusters 5 (left) and 6 (right) Number of cluster =                3 if Kd≤0.22 2 if 0.22 < Kd≤0.42 4 if 0.42 < Kd≤0.55 1 if 0.55 < Kd≤0.62 5 if 0.62 < Kd≤0.7 6 if Kd>0.7 (5.1) The observed relationship between daily clearness index and cluster suggests that K-means could produce clusters with days with a similar hourly solar radiation profile. Moreover, the relationship between the solar radiation distribution during a day and the daily clearness index value has been pointed out in (Bendt et al., 1981). However, solely one hourly profile cannot be used for all days due to the observed differences in different clusters. Using these two facts, using all the days of each cluster to estimate the solar radiation hourly profile for that cluster is proposed. Therefore, the hourly clearness index mean value and its standard deviation have been calculated for each cluster. Figures 5.11, 5.12 and 5.13 show these values. As can be observed, the standard deviation for most clusters and hours is not large and the results show that these values change with solar time, particularly for the hours at the start and end of the day. These results agree with those previously obtained in (Aguiar and Collares-Pereira, 1992). Moreover, these values decrease significantly for the hours with more radiation (central hours of day). This is the reason for proposing the use of these mean values as the hourly profile model for each cluster. 98 5. FORECASTING WITH CPDF 1 2 3 4 5 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Cluster Kt Figure 5.10: Daily clearness index, Kd, vs cluster to which the day belongs 5.6 Forecasting With the daily clearness index and the selected cluster obtained using equation 5.1, the hourly Khprofile of the cluster can be selected and used to forecast the hourly radiation values simply by multiplying each hourly Khvalue by its corresponding extraterrestrial radiation value, that are well known values because they are deterministic, only depending on the Sun position relative to Earth. This process is described in algorithm 3. Input :KdDaily clearness index; clearness index hourly profiles. Using Kdselect the cluster. Estimate solar global radiation hourly values using the clearness index hourly profile for the selected cluster and the extraterrestrial solar global radiation hourly values (Eq.4.6). Output: Hourly solar radiation values for day d Algorithm 3: Procedure for obtaining solar global radiation hourly values. Following the procedure described in algorithm 3, the estimated hourly radi- 5.6. FORECASTING 99 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kh,med Hour Kh,med StdDev 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kh,med Hour Kh,med StdDev Figure 5.11: Mean hourly clearness index and hourly standard deviation for clusters 1 (left) and 2 (right) 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kh,med Hour Kh,med StdDev 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kh,med Hour Kh,med StdDev Figure 5.12: Mean hourly clearness index and hourly standard deviation for clusters 3 (left) and 4 (right) ation values for all the recorded data are estimated for the testing dataset. The metrics for evaluating the hourly estimated modeled profiles are: •The standard deviation of daily profile hourly values for each cluster. •The difference between the total hourly solar radiation estimated, ˆ Gh, and the total hourly solar radiation received for the whole period of data normalized to the total hourly solar radiation received, Gh, according to the expression: Errorrad =Pm i=1 |Gh,i −ˆ Gh,i| Pm i=1 Gh,i 100(%) (5.2) 100 5. FORECASTING WITH CPDF 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kh,med Hour Kh,med StdDev 5 10 15 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Kh,med Hour Kh,med StdDev Figure 5.13: Mean hourly clearness index and hourly standard deviation for clusters 5 (left) and 6 (right) Eq.5.2 is used to estimate the energy error for each cluster to check the accuracy of these predictions. A naive persistent model that assumes that the hourly profile for a day is the same as the profile for the previous day is also used in order to evaluate whether the proposed model improves this naive model. Table 5.2 shows the obtained results. EE proposed EE naive EE random Cluster model (%) model (%) % energy clustering∗(%) % energy 1 10.5 16.5 26.3 24.9 15.5 2 36.8 71.0 4.9 29.7 15.7 3 49.1 193.8 1.3 25.3 15.4 4 25.0 35.9 10.5 24.2 17.0 5 5.0 10.8 37.1 17.9 17.1 6 4.0 12.3 19.9 22.6 19.3 All clusters 10.5 20.6 24.0 Table 5.2: Energy error (EE)(%) for each cluster when forecasting hourly solar global radiation with proposed model and with a persistent naive model and percentage of energy received in all days included in each cluster respect to total energy received.(∗clusters are randomly built) As can be observed, the proposed method is able to estimate the daily profiles of hourly global radiation with an error lower or equal to 5% for the 57% of energy received; this total increases to 84% with an error less than 11%. The highest error occurs for only 1.3% of energy received. The total error for all the clusters is 10.5% that is less than both naive model errors, 20.6%, and errors reported in previous works that range between 20 and 40%, as set out in the Introduction 5.7. CONCLUSIONS 101 section. These results indicate that it is possible to use the obtained daily profiles of hourly solar radiation distribution to forecast the hourly values of this variable. 5.7 Conclusions In this chapter a model to forecast hourly global radiation is presented. This model is based on clearness index, CPDF functions and K-means algorithm. Clearness index values are calculated and then CPDF funcions are constructed for every day in the data set. Then CPDF functions are clustered using K-means to create groups of similar days (days with similar solar radiation profile). As explaied in section 5.5, a relation between daily Kdand clusters can be established and this relation can be used to forecast hourly global radiation with a low error (5% for the 57% of energy received and 11% for the 84%). 103 Bibliography Aguiar, R., Collares-Pereira, M., 1992. T.a.g: A time dependent autoregressive gaussian model for generating synthetic hourly radiation. Solar Energy 49(3), 167–174. Bendt, P., Collares-Pereira, M., Rabl, A., 1981. The frequency distribution of daily insolation values. Solar Energy 27, 1–5. 110 6. FORECASTING MODELING •SVM-C + ANN •SVM-C + SVM-R All these models are estimated using a data subset, usually known as training data. The models obtained using training data are checked using test data to estimate the error of each combination. Training data set has 80% of samples while the test data set has the remaining 20% samples. The samples of each set are chosen randomly. Used data set ranges from October 2010 to December 2013. MATLAB 2014b implementations are used for the ANN and DT methods, while LibSVM library from Chih-Chung Chang and Chih-Jen Lin are used for both SVM-C and SVM-R methods (Chang and Lin, 2011). The ANN predictor is configured with input, hidden and output layers, and three neurons in the hidden layer. The hyperbolic tangent sigmoid transfer function is used in the input and hidden layer, while linear transfer function is used in the output layer. For the DT method, the MATLAB fitctree function is used with default options. SVM-C and SVM-R used here are implemented in an external library, which allows regression and multiple class classification with SVM to be performed. For the SVM-C method the paremeter nu has been set to 0.01 while for the SVM-R method the parameter γfor the kernel function has been set to 0.001. These values have been obtained after checking several values to obtain good results, as they are data dependent. 6.2.3 Hourly solar radiation forecasting Forecasting the hourly solar radiation values for the next day is performed in two phases that allow the next-day clearness index and the next-day type (cluster) to be determined; the hourly profiles of solar radiation are estimated using these two predicted parameters. The process for modeling hourly values of global solar radiation for a day is similar (note the difference between modeling and forecasting). The only difference is in the input data used; the values of the different meteorological parameters of the previous day and predictions of these parameters are used in forecasting while the values of these parameters for the current day are used in modeling process. Once the models for forecasting daily clearness index and cluster for next day have been built, the process for forecasting hourly values of global solar radiation 6.2. PROPOSED MODELS 111 for next-day is shown in Figure 6.3. First, the cluster and the daily clearness index, kd, are predicted using the selected models and the proposed independent variables. Then, the centroid of the predicted cluster (hourly values of index k∗ h,d) is selected and the predicted value of kdis added to each of the eight hourly values (the centroid) of k∗ h,d to obtain the values of hourly clearness index, kh,d. Finally, the hourly global horizontal solar radiation is calculated using Eq. 4.17 and Eq. 4.6. Input: Experiment 1: kd−1, Td−1, T9−12,d−1, T12−15,d−1, Hd−1, H9−12,d−1, H12−15,d−1, Pd−1, P9−12,d−1, P12−15,d−1 Experiment 2: kd−1, ˆ Td, ˆ T9−12,d, ˆ T12−15,d, ˆ Hd, ˆ H9−12,d, ˆ H12−15,d, ˆ Pd, ˆ P9−12,d, ˆ P12−15,d Using selected models: Estimate ˆ kd. Estimate cluster ˆcfor day d Use the centroid of ˆcas ˆ k∗ h,d for h= 8 to 15 Estimate ˆ kh,d. ˆ kh,d =ˆ k∗ h,d +ˆ kd Estimate ˆ Gh,d. and G0,h(Eq.4.6) ˆ Gh,d =ˆ kh,d ×G0,h Output ˆ Gh,d for h= 8 to 15 (Hourly values of global solar radiation for next day) Figure 6.3: Proposed procedure for forecasting next-day values of hourly solar radiation. 112 6. FORECASTING MODELING 6.3 Model and forecast evaluations 6.3.1 Experiments Two experiments have been performed using different input data sets depending on the used independent variables. In the first one, the values of the meteorological variables for the previous day are the independent variables of the analyzed models (forecasting). In the second one, the values of these variables for the same day to be forecasted are the independent variables except for the value of daily clearness index, where the value used corresponds to the previous day (modeling). The second experiment serves as test for using the forecasted meteorological data for the same day provided by a forecasting service. Table 6.1 summarizes the data used in each experiment (columns 1 and 2) including the input variables used for checking the models built in Experiment 2. The predicted data were collected from the Spanish Weather Service (AEMET) from November to December of 2015. In all cases, the dependent variables are the cluster and the daily clearness index. Independent variables Forecasting input Experiment 1 Experiment 2 variables (Exper.2) kd−1kd−1kd−1 Td−1, T9−12,d−1, T12−15,d−1Td, T9−12,d, T12−15,d ˆ Td,ˆ T9−12,d,ˆ T12−15,d Hd−1, H9−12,d−1, H12−15,d−1Hd, H9−12,d, H12−15,d ˆ Hd,ˆ H9−12,d,ˆ H12−15,d Pd−1, P9−12,d−1, P12−15,d−1Pd, P9−12,d, P12−15,d ˆ Pd,ˆ P9−12,d,ˆ P12−15,d Table 6.1: Independent variables used to build the models (cols 1 and 2) and forecasting input variables. Note: d−1 refers to the previous day and dcurrent day, 12 −15 refers to the mean value estimated using values from 12 to 15 hours, and 9−12 is for the mean value estimated using values from 9 to 12 hours, ˆ Xrefers to the forecasted value of X,Tis temperature, His humidity and Patmospheric pressure. 6.3. MODEL AND FORECAST EVALUATIONS 113 6.3.2 Error metrics Error metrics used here were described in 3.4. MSE,rMSE,MAE,rMAE and Forecast Skill are used to compare the different proposed models. The persistence model assumes that the conditions at the time of the forecast will not change, that is, the forecast for tomorrow are the values of today. The value of sis a measure of improvement over the persistence (higher values indicate better forecast skill, see 4.5.2). The forecast skill is independent of the specific meteorological or climatologic characteristics of the site under consideration (Coimbra and Kleissl, 2013). 6.3.3 Evaluation of clustering The first step is to cluster the observations using K-means to obtain four types of daily profiles. Each observation consists of 8 hourly values of variable k∗ h,d (Eq.6.1) estimated from the global solar radiation for each day. Table 6.2 shows the amount of elements (in percentages) in each cluster obtained with K-means and the percentage of total radiation received in the days included in the cluster. Cluster % Elements % Total radiation received in in cluster days included in cluster 1 66 72 2 14 11 3 12 10 4 8 7 Table 6.2: Number of observations in each cluster and percentage of radiation received in days included in cluster Clusters centroids are shown in Figure 6.4. Cluster 1 groups clear days, cluster 2 represents days with overcast variable among the day, cluster 3 represents days that are mainly overcast in the morning and cluster 4 represents days that are mainly overcast in the afternoon. Figure 6.5 shows the observations in each cluster. As it can be observed, most of the days are similar to the centroid of the cluster they belong to. These results confirm the hypothesis that it is possible to group hourly global radiation profiles throughout the day using 4 different classes. 114 6. FORECASTING MODELING Hours 12345678 K* h,d -0.3 -0.2 -0.1 0 0.1 0.2 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Figure 6.4: Values of index k∗ h,d (centroids) for each cluster. 6.3.4 Evaluation of models The second step in the data analysis is to determine the relationships between the daily clearness index and the meteorological parameters using two data mining techniques. This analysis is performed using two different procedures as explained in Section 6.2.2 and using two different input data sets (see Table 6.1). In the first procedure, the models to obtain clusters and daily clearness index are fitted separately and all training set observations are used to fit these models. In the second procedure, once the models to obtain the cluster are fitted and the clusters built, different models are fitted for the observations in each cluster to obtain the daily clearness index. Each procedure has been used for two different data sets. The first set of independent variables consists of variables of the previous day (Experiment 1). The second one corresponds to variables of the day to forecast or to model (Experiment 2), except for the daily clearness index, where the value of the previous day is used. The cluster estimating results for the test data set, when DT is used, are that 59.6% of data are correctly classified while the 65.6% of data are correctly classified 6.3. MODEL AND FORECAST EVALUATIONS 115 Hours 8.30 9.30 10.30 11.30 12.30 13.30 14.30 15.30 Kh,d* -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Hours 8.30 9.30 10.30 11.30 12.30 13.30 14.30 15.30 Kh,d* -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Hours 8.30 9.30 10.30 11.30 12.30 13.30 14.30 15.30 Kh,d* -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Hours 8.30 9.30 10.30 11.30 12.30 13.30 14.30 15.30 Kh,d* -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Figure 6.5: Values of index k∗ h,d for the days included in each cluster. Each line corresponds to the values of this index for a day. when SVM-C is used. A joint assessment of the models to estimate the cluster and the models to estimate daily clearness index has been performed. The proposed models have been evaluated using the MAE,rMAE and RMSE for hourly global solar radiation. First, models are evaluated using training set. Errors estimated for each experiment for training set are shown in Table 6.3. The combination of SVM-C+SVM-R in which a different SVM-R model for each cluster to estimate the daily clearness index is used, allows to obtain lower errors for both data inputs (Experiments 1 and 2). The value of rMAE for Experiment 1 is 13.3% while for Experiment 2 is 10.5%. The RMSE values are 22.2% and 16.8% for Experiment 1 and 2 respectively. All these values were obtained for training sets. 116 6. FORECASTING MODELING Experiment 1 Experiment 2 Method MAE(rMAE)RMSE MAE(rMAE)RMSE cluster kdW/m2(%) W/m2(%) W/m2(%) W/m2(%) DT ANN(1) 102 (19.0) 143 (26.7) 71 (13.3) 103 (19.3) DT SVM-R(1) 94 (17.2) 139 (25.4) 68 (12.5) 102 (18.9) SVM-C ANN(1) 97 (18.2) 137 (25.6) 68 (12.6) 98 (18.3) SVM-C SVM-R(1) 90 (16.4) 133 (24.3) 65 (11.9) 97 (17.9) DT ANN(2) 92 (17.5) 135 (25.6) 67 (12.6) 100 (18.8) DT SVM-R(2) 79 (14.3) 128 (23.3) 61 (11.3) 96 (17.8) SVM-C ANN(2) 87 (16.5) 129 (24.3) 63 (11.8) 94 (17.7) SVM-C SVM-R(2) 73 (13.3) 122 (22.2) 57 (10.5) 91 (16.8) Table 6.3: Training set: errors for each analyzed combination of methods for Experiment 1 (values of meteorological variables for the previous day) and Experiment 2 (values of meteorological variables for the same day to be modeled are the independent variables except the value of daily clearness index). (1) The same model is fitted for all observations (Procedure 1). (2) One model is fitted for the observations of each cluster (Procedure 2). 6.3.5 Evaluation of forecasts Test sets were used to evaluate the forecasts of different models. Models obtained in Experiment 1 were evaluated with the test set obtained from the measurement not used to build the models while models obtained in Experiment 2 need to be evaluated using forecasted data of the meteorological variables instead of the measured data. The data used in this case are the prediction obtained from the AEMET (Spanish Weather Service) from November to December 2015. Estimated errors for each experiment for test sets are shown in Table 6.4. The lowest rMAE is 16.7% in the case of Experiment 1 while it is 15.2% in the case of Experiment 2. These results are obtained when SVM-C and SVM-R are used for estimating the cluster and the daily clearness index respectively and when a different model is estimated for each cluster. Tables 6.5 and 6.6 show detailed MAE,rMAE and RMSE values for each experiment, procedure, cluster and data set (training and test). Not in all cases the SVM-C+SVM-R is the best model for each cluster, but is the best model taking into account all clusters. 6.3. MODEL AND FORECAST EVALUATIONS 117 Experiment 1 Experiment 2 MethodMAE(rMAE)RMSE MAE(rMAE)RMSE cluster kdW/m2(%) W/m2(%) W/m2(%) W/m2(%) DT ANN(1) 107 (19.9) 150 (27.9) 126 (24.3) 170 (32.7) DT SVM-R(1) 101 (18.3) 148 (26.9) 111 (21.3) 153 (29.4) SVM-C ANN(1) 104 (19.3) 145 (27.0) 121 (23.3) 159 (30.6) SVM-C SVM-R(1) 97 (17.6) 144 (26.0) 103 (19.8) 141 (27.1) DT ANN(2) 112 (21.1) 161 (30.1) 120 (23.1) 180 (34.6) DT SVM-R(2) 107 (19.3) 163 (29.3) 106 (20.3) 160 (30.8) SVM-C ANN(2) 102 (18.5) 150 (27.2) 105 (20.2) 152 (29.2) SVM-C SVM-R(2) 97 (16.7) 147 (25.3) 79 (15.2) 119 (22.9) Table 6.4: Test sets: errors for each analyzed combination of methods for Experiment 1 (values of meteorological variables for the previous day) and Experiment 2 (forecasts of meteorological variables for the same day to be modeled are the independent variables except the value of daily clearness index). (1) The same model is fitted for all observations (Procedure 1). (2) One model is fitted for the observations of each cluster (Procedure 2). MethodMAE(rMAE)RMSE cluster kdW/m2(%) W/m2(%) Cluster 12 3 4 1 2 3 4 Experiment 1 DT ANN(1) 94 (0.17) 83 (0.19) 125 (0.25) 158 (0.31) 134 (0.22) 110 (0.22) 168 (0.34) 204 (0.41) DT SVM-R(1) 84 (0.14) 76 (0.17) 125 (0.25) 163 (0.31) 127 (0.21) 103 (0.21) 169 (0.34) 210 (0.42) SVM-C ANN(1) 93 (0.16) 77 (0.18) 116 (0.23) 138 (0.27) 133 (0.22) 103 (0.21) 158 (0.32) 180 (0.36) SVM-C SVM-R(1) 82 (0.14) 70 (0.16) 116 (0.23) 144 (0.27) 126 (0.21) 94 (0.19) 159 (0.32) 187 (0.38) DT ANN(2) 84 (0.14) 77 (0.18) 122 (0.29) 139 (0.33) 130 (0.21) 105 (0.21) 155 (0.32) 177 (0.35) DT SVM-R(2) 70 (0.11) 67 (0.15) 110 (0.26) 124 (0.28) 128 (0.21) 94 (0.19) 141 (0.29) 158 (0.32) SVM-C ANN(2) 82 (0.14) 74 (0.17) 110 (0.26) 119 (0.29) 127 (0.21) 99 (0.20) 144 (0.29) 154 (0.31) SVM-C SVM-R(2) 68 (0.11) 61 (0.14) 97 (0.23) 101 (0.23) 126 (0.21) 86 (0.17) 128 (0.26) 130 (0.26) Experiment 2 DT ANN(1) 61 (0.11) 67 (0.15) 91 (0.20) 128 (0.26) 91 (0.15) 93 (0.19) 123 (0.25) 164 (0.33) DT SVM-R(1) 56 (0.10) 67 (0.15) 91 (0.20) 130 (0.26) 88 (0.14) 93 (0.19) 125 (0.25) 167 (0.34) SVM-C ANN(1) 60 (0.10) 63 (0.14) 83 (0.18) 115 (0.24) 88 (0.15) 87 (0.18) 114 (0.23) 148 (0.30) SVM-C SVM-R(1) 55 (0.09) 63 (0.14) 84 (0.18) 118 (0.24) 86 (0.14) 87 (0.18) 117 (0.24) 151 (0.30) DT ANN(2) 55 (0.09) 73 (0.17) 94 (0.22) 119 (0.27) 87 (0.14) 98 (0.20) 123 (0.25) 150 (0.30) DT SVM-R(2) 50 (0.08) 64 (0.15) 88 (0.20) 112 (0.26) 87 (0.14) 89 (0.18) 115 (0.23) 142 (0.29) SVM-C ANN(2) 54 (0.09) 65 (0.16) 83 (0.19) 104 (0.24) 85 (0.14) 88 (0.18) 111 (0.23) 135 (0.27) SVM-C SVM-R(2) 48 (0.08) 59 (0.14) 78 (0.180) 94 (0.22) 85 (0.14) 83 (0.17) 104 (0.21) 124 (0.25) Table 6.5: Training set: errors by cluster for each analyzed combination of methods for Experiment 1 (values of meteorological variables for the previous day) and Experiment 2 (values of meteorological variables for the same day to be modeled are the independent variables except the value of daily clearness index). (1) The same model is fitted for all observations (Procedure 1). (2) One model is fitted for the observations of each cluster (Procedure 2). 118 6. FORECASTING MODELING MethodMAE(rMAE)RMSE cluster kdW/m2(%) W/m2(%) Cluster 12341234 Experiment 1 DT ANN(1) 98 (0.17) 91 (0.21) 136 (0.28) 167 (0.32) 140 (0.23) 120 (0.24) 179 (0.37) 211 (0.42) DT SVM-R(1) 89 (0.15) 86 (0.19) 139 (0.28) 174 (0.32) 135 (0.22) 117 (0.24) 185 (0.38) 221 (0.44) SVM-C ANN(1) 94 (0.17) 90 (0.20) 133 (0.27) 165 (0.32) 135 (0.22) 117 (0.24) 174 (0.36) 208 (0.42) SVM-C SVM-R(1) 85 (0.15) 85 (0.19) 136 (0.27) 171 (0.32) 129 (0.21) 113 (0.23) 180 (0.37) 219 (0.44) DT ANN(2) 94 (0.16) 122 (0.26) 159 (0.38) 193 (0.50) 142 (0.23) 163 (0.33) 196 (0.40) 236 (0.47) DT SVM-R(2) 88 (0.14) 123 (0.26) 152 (0.35) 187 (0.47) 148 (0.24) 168 (0.34) 189 (0.39) 227 (0.46) SVM-C ANN(2) 100 (0.18) 107 (0.23) 170 (0.33) 80 (0.25) 148 (0.24) 132 (0.27) 204 (0.42) 103 (0.21) SVM-C SVM-R(2) 96 (0.16) 105 (0.23) 165 (0.32) 74 (0.25) 156 (0.26) 129 (0.26) 202 (0.41) 90 (0.18) Experiment 2 DT ANN(1) 68 (0.12) 77 (0.17) 111 (0.25) 149 (0.30) 101 (0.17) 104 (0.21) 148 (0.30) 191 (0.38) DT SVM-R(1) 63 (0.11) 76 (0.17) 107 (0.24) 150 (0.29) 97 (0.16) 105 (0.21) 145 (0.30) 195 (0.39) SVM-C ANN(1) 63 (0.11) 76 (0.17) 113 (0.25) 151 (0.30) 93 (0.15) 103 (0.21) 148 (0.30) 192 (0.39) SVM-C SVM-R(1) 59 (0.10) 75 (0.17) 110 (0.24) 151 (.030) 89 (0.15) 103 (0.21) 146 (0.30) 196 (0.39) DT ANN(2) 66 (0.11) 97 (0.22) 120 (0.30) 162 (0.41) 104 (0.17) 132 (0.27) 154 (0.32) 202 (0.41) DT SVM-R(2) 62 (0.11) 92 (0.20) 113 (0.28) 149 (0.37) 104 (0.17) 129 (0.26) 146 (0.30) 187 (0.37) SVM-C ANN(2) 74 (0.14) 79 (0.18) 115 (0.24) 95 (0.24) 114 (0.19) 106 (0.21) 145 (0.30) 114 (0.23) SVM-C SVM-R(2) 71 (0.13) 73 (0.16) 107 (0.22) 92 (0.24) 113 (0.19) 95 (0.19) 137 (0.28) 112 (0.22) Table 6.6: Test sets: errors by cluster for each analyzed combination of methods for Experiment 1 (values of meteorological variables for the previous day) and Experiment 2 (forecasts of meteorological variables for the same day to be modeled are the independent variables except the value of daily clearness index). (1) The same model is fitted for all observations (Procedure 1). (2) One model is fitted for the observations of each cluster (Procedure 2). 6.3. MODEL AND FORECAST EVALUATIONS 119 Locations rMAE(%) RMSE(%) s(%) Southern Spain (Lorenz et al., 2009) 12.2-20.4 20.8-37.1 2-36 persistence 16.6 32.1 Southern Spain (Perez et al., 2013) 13-21 22-29 17-37 persistence 19 35 Rome (Pierro et al., 2015) 16.3-22.3 26.4-34.1 26-29 persistence 25 47 SVM-C + SVM-R (Experiment 1) 16.7 25.3 23.0 persistence 25.8 32.7 SVM-C + SVM-R (Experiment 2) 15.2 22.9 43.9 persistence 24.0 40.8 Table 6.7: Comparison of the performance of the SVM-C+SVM-R models. Finally, to compare the results obtained for the model that gives the best results for both experiments with previously proposed models, the forecast skills over 24 hour persistence forecasts (s) were also estimated. Table 6.7 compares the benchmark accuracy range previously found in (Lorenz et al., 2009) and (Perez et al., 2013) for the southern Spain and (Pierro et al., 2015) for Rome (Italy) with the performance obtained by proposed model. Regarding the results obtained in Experiment 1 the proposed model gives values of rMAE and RMSE inside the benchmark of the values obtained for Spanish locations and similar to those obtained for Rome, except in this case for RMSE that is slightly lower. The value of forecast skill over 24 hour persistence obtained is 23% and is also inside the benchmark of Spanish locations but again lower that values for Rome. It should be noted that the RMSE of the persistence model for this location is greater than the value obtained for Spanish locations. This suggest that the data used for Spanish locations correspond to a period with more sunny days, as pointed out by authors. The errors obtained in Experiment 2 are also inside the benchmark of the values obtained for Spanish locations, rMAE is 15.2% and RMSE is 22.9%. The value of parameter sis 43.9%, greater than the previosly reported. This could be due to the fact that only data from November and December were used where the number of sunny days is lower and the performance of persistence model is worse. 126 7. MONITORING PV PLANTS Figure 7.1: Main system architecture enables the existence of applications with different functionality using the same OPC servers: synchronize data, display data from OPC servers, evaluating data, etc. Scheduled tasks are used to read and evaluate data from plants running on a machine and as a result the users (plant owners, staff responsible for maintenance, etc.) receive e-mail messages and SMS messages that keep them informed of the status of their facilities, depending on the level of access that each one is assigned. 7.2 Assessing photovoltaic facilities Assessment models can be used to supervise the performance of the solar power plants. As proposed in (Mart´ınez Marchena I, 2014), evaluation parameters can be implemented easily as performance indexes: daily final yield and daily enegy balances (daily energy produced and daily energy forecasted). The daily final yield, Yf,day, is defined as the useful output energy of the PV per kWpinstalled: Yf,day =Eday PST C (7.1) 7.2. ASSESSING PHOTOVOLTAIC FACILITIES 127 where PST C is the nominal power of the installed photovoltaic array at standard test conditions (1 kM/m2of solar irradiance and 25°cell temperature). The useful daily output energy or daily energy supplied by a solar power plant is defined as: Eday =Zday PAC(t)d(t)≈ n X j=1 Pj AC∆t(7.2) where nis the number of measurements along the day and Pj AC are the recorded values of the power generated at the inverter output. The proposal is to assess photovoltaic facilities in two phases. In the first phase, the evaluation is done by comparing the values obtained for the daily final yield and the useful daily output energy when they are estimated using two different methods. On the one hand, these values are estimated directly using the recorded values in the inverter. In the other hand, these parameters are estimated using a model for the performance of photovoltaic plant and the forecasted values of hourly global radiation. In this phase only daily values are compared. Taking into account the results of this comparison, a second phase could be activated in order to compare hourly values. 7.2.1 Models for estimating the performance of a PV system As it has been reported previously, see for instance Ayompe et al. (2011), the inverter power output has a linear relationship with solar irradiance if it is not considered the effect of the temperature. This fact can be observed in Fig. 7.2. We propose to include this effect in the model that allows to estimate the hourly energy generated at the output of the inverter (E∗ h) as an extension of the estimation of the power generated (P∗ AC), according to the expression: E∗ h=Eh,ST C G∗ h,β 1000(1 + γ(Th,mod −25))GL (7.3) where, G∗ h,β is the forecasted hourly global radiation on the surface of the modules, βis the inclination of the modules, γis the temperature coefficient of 128 7. MONITORING PV PLANTS Figure 7.2: Instantaneous values of solar radiation and photovoltaic power Pm,Th,mod is the mean hourly module temperature and GL is the global losses coefficient of the system. Eq.7.3 is obtained from the expression proposed by Osterwald (1986) considering hourly values. The proposed expression includes not only the losses produced by the temperature but also other losses (soiling, spectral losses and so on). The predictions of hourly global radiation on the surface of the modules are estimated in the following way: 1. Forecasting the hourly global radiation values on horizontal surface with the model proposed in Chapter 6.2.2 and meteorological variables recorded by weather agencies. 2. Obtaining the values of direct, diffuse and reflected hourly global radiation on horizontal surface. 3. Obtaining the values of global radiation on the surface of the modules. The results obtained with this model can be used for estimating the daily energy supplied to the grid, E∗ day, which is calculated using the expression 7.4: 7.2. ASSESSING PHOTOVOLTAIC FACILITIES 129 E∗ day =Zday P∗ AC(t)d(t)≈ n X j=1 E∗(j) h(7.4) where nis the number of hours for a day. The proper operation of plants can be evaluated by comparing the values of E∗ day with the values of daily energy produced, Eday, obtained with Eq.7.2. Similarly, for detecting problems in plants operation the value of daily final yield, Eq.7.1, has been compared to the estimated daily final yield, Y∗ f,day, calculated with the expression 7.5: Y∗ f,day =E∗ day PST C (7.5) 7.2.2 Statistical models for assessing solar plants With the estimated values described in the previous section and the corresponding measured values, the system is capable of checking the performance of the plants by using the methodology described in (Mart´ınez Marchena I, 2014). The mean values of previously described parameters for photovoltaic solar plants of the same technology (modules and inverters) has been estimated as the standard deviation of them. For evaluating the operation of each plant, a statistical analysis of the differences between the values of estimated parameters and measured parameters is implemented. Thus, for each new recorded or calculated value of some parameter Xand its corresponding estimation, X∗, at hour ithe difference among them is analyzed according to the criteria (significance level 5%): d(i) X=X(i)−X∗(i); if d(i) X/∈[−1.96ˆσ, +1.96ˆσ]then mark (i),(7.6) where ˆσis the sample standard deviation of X. All the marked values (i) will indicate problems of operation. For the generated energy, E(i) d, and its corresponding estimation using Eq.7.3, 130 7. MONITORING PV PLANTS E∗(i) d, the criteria is as follows: d(i) PAC =E(i) d−E∗(i) d; if d(i) PAC /∈[−1.96ˆσ, +1.96ˆσ]then mark (i),(7.7) For the daily final yield, Yf,day, the expression used to decide if there is any operation problem taking into account the expression Eq.7.6 particularized for this parameter is the following: d(i) Yf,day =Y(i) f,day −Y∗(i) f,day; if d(i) Yf,day /∈[−1.96ˆσ, +1.96ˆσ]then mark (i),(7.8) Fig. 7.3 shows the algorithm description for the daily assessment of a PV plant. Recorded measures are used as input and estimated measures are calculated from these inputs. Then a test is performed to check if the difference between measured and estimated daily yield is statistically significant. If these values are statistically equals (for a specified significance level) the algorithm ends without detecting problems, while if they are different, each hourly value obtained from recorded power values is compared with the corresponding estimated and if the difference is statistically significant, the corresponding hour is marked. The algorithm ends with operational problems detected and the marked hours are returned. 7.3 Example scenario From the architecture described in the previous section, monitoring and evaluating systems for plants are developed, integrating in the software of each plant components corresponding to those subsystems which are present. As example of developed systems, this section presents an application which has been developed for the Agencia Municipal de la Energ´ıa del Ayuntamiento de M´alaga. The application allows remote monitoring and evaluation of 40 different solar power facilities, mostly located in public schools and government buildings. All facilities are located in M´alaga. Inverter technologies present in the different plants include Sunways, Ingeteam, Mastervolt, SMA, PowerOne, SolarMax and Fronius. On each power plant floor 7.3. EXAMPLE SCENARIO 131 start Input: Recorded measures Estimate virtual measurements Y(i) f,day −Y∗(i) f,day ∈ [−1.96ˆσ, +1.96ˆσ] i= 0 (hours) E∗(i) h−E(i) h∈ [−1.96ˆσ, +1.96ˆσ] mark hour i i < total measures Operation problems detected No problems detected end no no yes yes yes Figure 7.3: Flowchart for assessing solar energy plants. 132 7. MONITORING PV PLANTS the existing monitoring system consists of a datalogger that depends on the characteristics of the inverter manufacturer (Meteocontrol, Solarlog, Mastervolt, etc.). With the objective of standardize the acquisition of data, OPC servers have been developed to allow standard access to available data on each technology to integrate them into the system. The values obtained for serveral days without any detected problems are shown in Fig. 7.4. Figure 7.4: Daily parameters In figure 7.5 and 7.6 operational problems detected with developed software in two different plants are shown. In the first one, connection problems detected in one of the inverters are shown, while in the second figure, joint analysis of operation of two inverters from the same power plant allows to detect a malfunction of one of them. 7.4 Conclusions This chapter presents an assessing model for photovoltaic facilities. This model has been integrated in a monitoring framework for solar photovoltaic plants. With 7.4. CONCLUSIONS 133 Figure 7.5: Detecting problems in the inverter the solution proposed applications that integrate a single tool monitoring photovoltaic systems connected to network with inverters or data acquisition systems from different technologies, meter reading and production analysis procedures, evaluation, fault detection and alarm generation of the plant. This proposal enables evaluation of photovoltaic plants remotely using a single program, and without relying on software developed by inverters manufacturers. In addition, this assessment allows rapid action in a plant when malfunctions are detected in the same, thanks to the sending alerts system and remote access to plant data which facilitates maintenance tasks and increases the profitability of photovoltaic systems. 134 7. MONITORING PV PLANTS Figure 7.6: Detecting problems in one plant with two inverters 135 Bibliography Ayompe, L., Duffy, A., McCormack, S., Conlon, M., 2011. Measured performance of a 1.72 kw rooftop grid connected photovoltaic system in ireland. Energy Conversion and Management 52 (2), 816 – 825. URL http://www.sciencedirect.com/science/article/pii/ S0196890410003730 Ciesielska, J., Concas, G., Despotou, E., Fontaine, B., Garbe, K., Fraile-Montoro, D., et al., 2011. Global market outlook for photovoltaics until 2015. European Photovoltaic Industry Association (EPIA). Figueiredo, J. M., da Costa, J. M. G. S., May 2008. An efficient system to monitor and control the energy production and consumption. In: 2008 5th International Conference on the European Electricity Market. pp. 1–6. Guozhen, H., Tao, C., Changsong, C., Shanxu, D., 2009. Solutions for scada system communication reliability in photovoltaic power plants. In: Power Electronics and Motion Control Conference, 2009. IPEMC’09. IEEE 6th International. IEEE, pp. 2482–2485. Holley, D. W., 2004. Understanding and using opc for maintenance and reliability applications. Computing & Control Engineering 15 (1), 28 – 31. URL http://0-search.ebscohost.com.jabega.uma.es/login.aspx? direct=true&db=a9h&AN=12650871&lang=es&site=ehost-live&scope=site Kalaitzakis, K., Koutroulis, E., Vlachos, V., 2003. Development of a data acquisition system for remote monitoring of renewable energy systems. Measurement 34 (2), 75 – 83. URL http://www.sciencedirect.com/science/article/pii/ S0263224103000253 Mart´ınez Marchena I, Sidrach-de-Cardona M, M.-L. L., 2014. Framework for monitoring and assessing small and medium solar energy plants. Sol. Energy Eng.