Advanced forecast and scheduling of power systems with highly variable sources

Advanced forecast and scheduling of power systems with highly variable sources PhD Program in Sustainable Energy Systems Pedro Miguel Neves da Fonte, M.Sc. Dissertation submitted to the Faculty of Engineering of University of Porto in partial fulfilment of the requirements for the degree of Doctor of Philosophy Supervisor: Professor Fernando Pires Maciel Barbosa, Ph. D. Co-supervisor: Professor Cláudio Domingos Martins Monteiro, Ph. D Department of Electrical and Computer Engineering, Faculty of Engineering, University of Porto January, 2015 © Pedro Miguel Neves da Fonte, 2015 À Paula, à Carolina e ao Rodrigo iii Acknowledgments It was a long, challenging and sometimes difficult path that allowed me to conclude this thesis. This thesis was only possible thanks to many people and institutions, who helped me in several ways to reach this goal. I owe my thanks to all of them. First of all, I would like to express my gratitude to my scientific adviser, Professor Fernando Maciel Barbosa; to his commitment and dedication as well as his valuable advice. It was his constant support and advice that kept me moving forward and steadfast, even in times of great frustration. A very special thanks to my co-adviser Professor Cláudio Monteiro, for all his support from the beginning and throughout this work and, in particular for having placed at my disposal, his insight and scientific ability. Without his advice, I would not have reached my goals. I would like to highlight both his human and scientific qualities. To the São Miguel’s power company, EDA – Electricidade dos Açores, a special acknowledge for providing the data used in this work. To the Smartwatt team, a great acknowledgement for their team spirit, and all support. I have to give special thanks to Bruno Santos for his infinite patience in helping me during the implementation of the case studies. Without his efforts, this work would have had a different ending. With his scientific and human abilities, I predict that he will have a future full of success. To my PhD colleagues, Pedro Almeida, Manuel Rocha, Julia Vasiljevska, Antero Moreira da Silva, Joel Soares and Pedro Melo, who have demonstrated a spirit of friendship and companionship. This spirit will never be forgotten. My thanks to professor João Peças Lopes, for his kindness and comprehension revealed when he received me for the first time in the PDSEE program. To all my colleagues of the ISEL’s Electrical Machines Group. First my thanks goes to Professor José Carlos Quadrado, who influenced me to follow a teaching and researching career. Very special thanks to Rita Pereira and Ricardo Luís, for all their incredible support, their tolerance and encouragement during hard times. To them I wish the biggest success in their PhD research. To Professor Sérgio Abrantes Machado whose technical and human competencies are a source of inspiration and an example to follow. A special thanks for the remaining colleagues of the ISEL’s Electrical Machines Group, for all the words of encouragement to pursue the work. The most special acknowledgment goes to my family. To my parents and brother, who throughout my life inspired my accomplishments, on a personal and professional level. I would iv like to make my greatest thanks to the encouragement and trust that they have in me. I hope I deserve that trust. Lastly I thank Paula, my always present wife, my life mate for more than twenty years, who has always supported me in both good and bad times, and always in her own detriment. During all these years she has been my pillar of support. Finally to my children, Carolina and Rodrigo, hoping that, someday, they will be able to understand all the moments during which I was not present and I hope that this work can one day serve as an inspiration to their long lives. To all those that I did not mention but are in my mind I give great thanks. A todos, um grande Obrigado. v Abstract Management of electric energy production is one of the most important issues in electric grids operation. With this management it is intended to feed the load ensuring the generation/consumption equilibrium in a most economical way, simultaneously respecting technical constraints. In addition to this technical/economical management it is mandatory to ensure the system reliability in order to safeguard the continuity of service in case of some fault, or catastrophe of some of the power system components, whether they are on generation level, transportation or loads. All this management is done through the unit commitment and economic dispatch, to be able to decide which generation units have to be connected to the power grid as well as the allocated load to each one. This process, based in mature and well dominated technology, became a new challenge with the massive introduction of electric energy generation based in renewable power sources. These, are predisposed to depend on variable and hard to control resources, adding a large amount of uncertainty to the decision process. This problem is enhanced in systems with low rated power, as is the case of islands without connection to the large continental grids, without storage capacity or quick starting generation units. Generally this problem is addressed by two different approaches which, in the end, complement each other. One is at a forecasting level; researching improvements in performances of forecasting as well as the apprehension of the renewable production uncertainties. The second pursues the development of scheduling models which potentiate information obtained by the forecast, as well as the increment of velocity. Most recent researching works focus on the application of stochastic programming models in order to decide the unit commitment of the units as well as in production and reserve’s allocation and security levels. This work intends to develop a complete methodology including the forecasting of renewable generation and load, with respective characterization of the uncertainty, complemented with scheduling based in a risk assessment model. The case study refers to the production power system of São Miguel Island (Azores, Portugal). During the development of this work, an overview about the techniques and mathematical formulation of scheduling models (deterministic and stochastic) as well as some state of the art solving techniques was done. It is also presented an overview concerning the renewable production and load forecasting, together with models to represent the uncertainty. It was done a detailed study about some of these techniques, namely concerning the choice of forecasting models and explanatory variables as well as a short sensibility analysis about the influence of explanatory variables. Thus, it was developed a model of aggregate forecasts in order to indicate the thermal production necessities. Then it is presented an original contribution of a scheduling vi model based on risk assessment based on formulation of the production system of São Miguel. Also, characterization of thermal production is done, introducing the concept of equivalent optimal generation unit. In parallel, it is presented a metaheuristic based on a cloud of particles to solve the economic dispatch problem, with its performances being compared with some results presented in biography. At the end it is shown a real case study where detailed explanation is presented, under a real context, the developed methodology. Finally, it is done a comparison with the scheduling proposed by the system operator of São Miguel Island and conclusions are drawn. In short, this work proposes a complete methodology since the forecasting up to the generation scheduling for an isolated system with large penetration of renewable generation. Keywords: Metaheuristics; Power generation scheduling; Power forecasting; Risk assessment; Power system operation in islands vii Resumo A gestão da produção de energia elétrica é uma das questões mais importantes na operação de redes de energia. Com esta gestão pretende-se alimentar a necessidade de carga mantendo o equilíbrio produção/consumo da forma mais económica e ao mesmo tempo respeitando as restrições técnicas. Para além desta gestão técnico/económica é necessário garantir a fiabilidade do sistema de modo a garantir a continuidade de serviço no caso de mesmo no caso de alguma falha, desde que não catastrófica, de alguns componentes do sistema de potência, sejam eles a nível da geração, transporte ou cargas. Toda esta gestão é feita através do comissionamento de grupos e despacho económico de modo a decidir quais a unidade que devem estar ligadas à rede bem como a alocação de carga a cada uma delas. Todo este processo baseado em tecnologia madura e bem cimentada tornou-se num desafio com a introdução em massa de geração de energia elétrica com base em energias renováveis. Estas tendencialmente dependem de recursos variáveis e difíceis de controlar. Nesse sentido foi adicionado uma grande quantidade de incerteza ao processo de decisão. Este problema é potenciado em sistemas com baixa potência instalada como é o caso de ilhas sem ligação às grandes redes continentais e sem capacidade de armazenamento de energia ou unidades produtoras com arranques rápidos. Este problema geralmente é abordado por duas perspectivas diferentes que no fim se complementam. Uma a nível da previsão, com a procura de uma melhoria nas performances da previsão e captura da incerteza da produção renovável. A outra no desenvolvimento de modelos de comissionamento para potenciar a informação obtida pelas previsões bem como o incremento da sua rapidez. Os trabalhos mais recentes apostam na aplicação de modelos de programação estocástica com vários níveis de modo a decidir o comissionamento de grupos, bem como a alocação de produção e reservas, e níveis de segurança. Neste trabalho desenvolveu-se uma metodologia completa que passa pela previsão da produção renovável e carga com a respetiva caracterização da incerteza e pelo comissionamento da geração baseado numa análise de riscos. O caso estudado é referente ao sistema produtor da ilha de São Miguel nos Açores (Azores/Portugal). Durante o desenvolvimento deste trabalho foi feito uma visão geral acerca das técnicas e formulações matemáticas dos modelos de agendamento (determinístico e estocástico) bem como de algumas técnicas de resolução presentes no estado-da-arte. É também apresentada uma visão geral acerca das técnicas de previsão de produção renovável e carga, bem como dos modelos para representar a sua incerteza. Foi feito um estudo detalhado acerca de algumas dessas técnicas, nomeadamente acerca da escolha dos modelos de previsão e das variáveis explicativas bem como uma pequena análise de sensibilidades acerca da influência das variáveis explicativas. xiv Figure 3.8 – Resulting “theoretical” power curve .................................................................... 80 Figure 3.9 – Wind power limit, wind power production and theoretical power production .... 80 Figure 3.10 – Probability distribution in function of smoothing parameter h .......................... 81 Figure 3.11 – Conditional probability distribution of wind power forecast ............................. 82 Figure 3.12 – Hourly measured and forecasted wind power generation .................................. 83 Figure 3.13 – Hourly average hydro power production and hourly forecasted precipitation .. 84 Figure 3.14 – Daily average hydro power precipitation and HPP............................................ 85 Figure 3.15 – Average hydro power production and HPP ....................................................... 86 Figure 3.16 – Measures vs hydro power forecasted and HPP .................................................. 87 Figure 3.17 – Hourly measured and forecasted hydro power .................................................. 88 Figure 3.18 – Hourly average geothermal power production all over 2012 ............................ 89 Figure 3.19 – Measured and forecasted geothermal with uncertainty ..................................... 90 Figure 3.20 – Behaviour of load for each hour all over 2012 (all days) .................................. 91 Figure 3.21 – Behaviour of load for each hour all over the weekends of 2012 ....................... 91 Figure 3.22 – Behaviour of load for different first weeks of 2012 .......................................... 91 Figure 3.23 – Relation between load and temperature (Crete Island) ...................................... 92 Figure 3.24 – Relation between load and temperature (S. Miguel Island) ............................... 92 Figure 3.25 – Probability distribution of load .......................................................................... 93 Figure 3.26 – Measured and forecasted load with uncertainty................................................. 94 Figure 3.27 – Measured and forecasted net load with uncertainty .......................................... 95 Figure 3.28 – Diagram of the net load forecasting................................................................... 96 Figure 3.29 – Examples of Beta pdf fitting to the probability distributions ............................. 97 Figure 4.1 – Proposed methodology ...................................................................................... 102 Figure 4.2 – Specific fuel consumption for each type of thermal units ................................. 103 Figure 4.3 – Fuel consumption for each type of thermal units .............................................. 104 Figure 4.4 – Power production limits of the GENSET’s ....................................................... 105 Figure 4.5 – Equivalent optimal generation unit .................................................................... 108 Figure 4.6 – SCO’s flowchart ................................................................................................ 109 Figure 4.7 – Example of quadratic regression from cloud of particles fitness values (one dimension) .............................................................................................................................. 111 Figure 4.8 – Evaluation of ϕ(i,j), ............................................................................................ 112 Figure 4.9 – Contour plot of sphere function with 2 dimensions .......................................... 113 Figure 4.10 – Particle cloud behaviour for the variable x1 ..................................................... 114 Figure 4.11 – Particle cloud behaviour for the variable x2 ..................................................... 114 Figure 4.12 – Inverted sigmoid behaviour for different values of ϕ ∆ ................................... 115 Figure 4.13 – Inverted sigmoid behaviour with K=1,5 .......................................................... 116 xv Figure 4.14 – SCO’s performance solving Rastrigin’s function ............................................ 118 Figure 4.15 – Behaviour of cloud variance ............................................................................ 118 Figure 4.16 – Comparative minimum, maximum and average solutions (first case study) ... 121 Figure 4.17 – Convergence behaviour of SCO for a 6-units problem ................................... 122 Figure 4.18 – Comparative minimum, maximum and average solution (second case study) 123 Figure 4.19 – Convergence behaviour of SCO for a 15-units problem ................................. 124 Figure 4.20 – a) non-convex cost function with valve point effect; b) its derivative ............ 125 Figure 4.21 – Convergence behaviour of SCO for a 40-units problem ................................. 126 Figure 4.22 – Comparative minimum, maximum and average solutions (third case study) .. 126 Figure 4.23 – Convergence behaviour of SCO for a 10-units problem ................................. 127 Figure 4.24 – Example of single period unit commitment based on risk assessment ............ 129 Figure 4.25 – Uncertainty associated to a specific committed GENSET (Load-RES) .......... 131 Figure 4.26 – Cumulative distribution function associated to a specific GENSET (Load-RES)131 Figure 4.27 – Pdf of net load with (L-(H+GEO)) and without wind curtailment (L-RES) .... 133 Figure 4.28 – Cdf of net load with (L-(H+GEO)) and without wind curtailment (L-RES) .... 133 Figure 4.29 – Probability of thermal units work below the minimum (pdf) .......................... 134 Figure 4.30 –Probability of thermal units work below the minimum (cdf) ........................... 134 Figure 4.31 – Inverse cdf of net load with (L-(H+GEO)) and without wind curtailment L-RES ................................................................................................................................................ 135 Figure 4.32 – State space diagram of a repairable component ............................................... 137 Figure 4.33 – Unit commitment via forward dynamic programming .................................... 140 Figure 4.34 – Flowchart of proposed methodology ............................................................... 141 Figure 5.1 – Estimated wind power without curtailment and measured ................................ 145 Figure 5.2 – Measured and predicted values of net load for the case study .......................... 146 Figure 5.3 – Reliability diagram of the probabilistic net load forecast .................................. 147 Figure 5.4 – Geothermal power production (January 1st up to March 3rd)............................ 147 Figure 5.5 – Reliability diagram (with the removed values) .................................................. 148 Figure 5.6 – Net load forecasting reliability diagram (deviations from “ideal” reliability) ... 148 Figure 5.7 – Sharpness diagram of the probabilistic net load forecast .................................. 149 Figure 5.8 – Resolution of the probabilistic net load forecast ............................................... 150 Figure 5.9 – Net load forecasting for 24 hours ahead ............................................................ 152 Figure 5.10 – Risk associated to GENSET 1GS_1GB operation ............................................ 153 Figure 5.11 – Percentage of risk costs associated to GENSET 1GS_1GB .............................. 154 Figure 5.12 – Cost associated to GENSET 1GS_1GB ............................................................ 154 Figure 5.13 – Risk assessment for all GENSET’s at 8:00...................................................... 155 Figure 5.14 – Percentage of total risk costs for all GENSET’s at 8:00 .................................. 155 xvi Figure 5.15 – Total risk cost to all GENSET at 8:00 ............................................................. 156 Figure 5.16 – Single-period unit commitment ....................................................................... 157 Figure 5.17 – Upward and downward forecasted reserves .................................................... 158 Figure 5.18 – Single-period unit commitment and perspective measurement of net load ..... 158 Figure 5.19 – Net load point forecast and perspective of measurement without curtailment 159 Figure 5.20 – Upward and downward measured reserves ..................................................... 160 Figure 5.21 – Forecasted and measured upward reserves ...................................................... 160 Figure 5.22 – Forecasted and measured downward reserves ................................................. 160 Figure 5.23 – Measures of net load and committed GENSET limits .................................... 161 Figure 5.24 – Measures and perspective of wind power production (WPP) with limitation . 162 Figure 5.25 – System operator upward and downward reserves ............................................ 162 Figure 5.26 – Comparison between upward reserves ............................................................ 163 Figure 5.27 – Comparison between downward reserves ........................................................ 163 Figure 5.28 – Values outside the GENSET limits, in percentage of total net load ................ 164 Figure 5.29 – Risk costs for both approaches ........................................................................ 165 Figure 5.30 – Comparison between total costs for both approaches ...................................... 165 Figure 5.31 – Multi-period unit commitment ......................................................................... 167 Figure 5.32 – GENSET power limits and measured net load ................................................ 168 Figure 5.33 – Upward and downward in multi-period approach ........................................... 169 Figure 5.34 – Risk assessment and system operator downward reserves for a multi-period approach ................................................................................................................................. 170 Figure 5.35 – Multi-period unit commitment ......................................................................... 171 Figure 5.36 – Wind energy curtailed and produced below the minimum by system operator172 Figure 5.37 – Wind energy curtailed and produced below the minimum by the risk assessment approach ................................................................................................................................. 172 Figure 5.38 – Multi-period operation costs ............................................................................ 172 Figure 5.39 – Risk costs in function of net load .................................................................... 174 Figure AI.1 – Specific consumption of each GENSET.......................................................... 197 xvii List of tables Table 2.1 – Overview of techniques of scheduling with uncertainty ....................................... 29 Table 2.2 – Overview of stochastic programming ................................................................... 41 Table 3.1 – Rated power of each power source ....................................................................... 70 Table 3.2 – Bandwidth hj for wind power forecast Kernel’s ................................................... 82 Table 3.3 – Small hydro power capacity in São Miguel Island ............................................... 83 Table 3.4 – Parameters used to hydro power forecast .............................................................. 87 Table 3.5 – Bandwidth for hydro power forecast Kernel’s ...................................................... 87 Table 3.6 – Geothermal power capacity in São Miguel Island ................................................ 89 Table 3.7 – Bandwidth for geothermal power forecast Kernel’s ............................................. 89 Table 3.8 – Bandwidth for load forecast Kernel’s ................................................................... 93 Table 3.9 – Performances of point forecasts ............................................................................ 94 Table 4.1 – Specific consumption for thermal units .............................................................. 103 Table 4.2 – Possible combinations of GENSET’s ................................................................. 105 Table 4.3 – Parameters of the SCO to solve Rastrigin’s function ......................................... 117 Table 4.4 – Generating unit’s data (first case study) .............................................................. 121 Table 4.5 – Generating unit’s prohibited zones (first case study) .......................................... 121 Table 4.6 – Results obtained (6-units 1263 MW) (Best individual) ...................................... 122 Table 4.7 – Generating unit’s prohibited zones (second case study) ..................................... 123 Table 4.8 – Generating unit’s data (second case study) ......................................................... 123 Table 4.9 – Results obtained (15-units 2630 MW) (Best individual) .................................... 124 Table 4.10 – Results obtained (40-units 10500 MW) (Best individual) ................................ 126 Table 4.11 – Results obtained (10-units 2700 MW) (Best individual) .................................. 128 Table 5.1 – Parameters for each case study............................................................................ 151 Table 5.2 – Probability and risk costs for GENSET 1GS_1GB .............................................. 155 Table 5.3 – Measured and forecasted energy below and above the GENSET limits ............ 159 Table 5.4 – Number of hours operating outside the GENSET’s limits .................................. 164 Table 5.5 – Values outside the GENSET limits, in percentage of total net load ................... 164 Table 5.6 – Comparison between total costs for both approaches ......................................... 166 Table 5.7 – Number of start-ups for both approaches ............................................................ 167 Table 5.8 – Average reserves in multi-period UC.................................................................. 169 Table 5.9 – Values outside the GENSET limits for multi-period UC (in percentage of total net load ) ...................................................................................................................................... 171 Table 5.10 – Number of hours operating outside the GENSET’s limits ................................ 171 Table 5.11 – Comparison between total costs for both approaches ....................................... 173 Table 5.12 – Comparison regarding load shed ....................................................................... 173 xviii Table 5.13 – Comparison regarding wind curtailment ........................................................... 173 Table 5.14 – Comparison regarding thermal units working below the minimum ................. 174 Table AI.1 – Technical characteristics of thermal units ......................................................... 197 Table AII.1 – Generating unit’s data (third case study) ......................................................... 199 Table AII.2 – Results obtained (40-units 10500 MW)(Best individual) ................................ 200 Table AII.3 – Generating unit’s data (fourth case study) ....................................................... 200 Table AII.4 – Generating unit’s data (fourth case study) ....................................................... 201 Table AII.5 – Cost functions parameters of GENSET’s ........................................................ 202 Table AII.6 – Cost functions parameters of GENSET’s ........................................................ 203 Table AII.7 – Cost functions parameters of GENSET’s ........................................................ 204 Table AIII.1 – Unit commitment for single period (system operator) ................................... 205 Table AIII.2 – Unit commitment for single period (Risk assessment) .................................. 205 Table AIII.3 – Unit commitment for multi-stage period (Risk assessment) .......................... 206 xix List of abbreviations and symbols AGC – Automatic Generation Control ALADIN – Aire Limitee Adaptation Dynamique development InterNational ANFIS – Artificial Neuro Fuzzy Inference System ANN – Artificial Neural Networks ANN_CC – Cascade Correlation ANN_CG – Conjugate-Gradient ANN_LM – Levenberg-Marquardt ARMA – Auto Regressive Moving Average ARIMA – Autoregressive Integral Moving Average ARMAX – Auto-Regressive with Exogenous Input BMGP – Below Minimum Generation Probability BP – Back propagation CCPSO – Chaotic sequences and Crossover operation PSO CEP – Classical Evolutionary Programming CFD – Computational Fluid Dynamic CGRG – Central Geotérmica Ribeira Grande CGPV – Central geotérmica do Pico Vermelho CHTN – Central Hídrica dos Túneis CHTB – Central Hídrica dos Tambores CHFN – Central Hídrica da Fábrica Nova CHCA – Central Hídrica do Canário CHFR – Central Hídrica Foz da Ribeira CHRP – Central Hídrica Ribeira da Praia CHSC – Central Hídrica do Salto do Cabrito COPSO – Crossover Operation PSO CSPSO – Chaotic Sequences PSO DED – Dynamic Economical Dispatch DP – Dynamic programming ED – Economic dispatch EDA – Eletricidade dos Açores EP – Evolutionary Programming ESO – Evolutionary Strategy Optimization EV – Expected Value FS-SCUC – Full Scenario-Security Constrained Unit Commitment FEP – Fast Evolutionary Programming xx GA – Genetic Algorithms GENSET – Generators Set HPP – Hydrological Power Potential IFEP – Improved Fast Evolutionary Programming IPSO – Improved Particle Swarm Optimization ISO – Independent System Operator ITS – Improved Taboo Search KDE – Kernel Density Estimator LHS – Latin Hypercube Sampling LS – Load Shed LOLE – Loss Of Load Expectation LOLP – Loss Of Load Probability LOWP – Loss Of Wind Probability MAE – Mean Absolute Error MAPE – Mean Average Percentage Error ME – Mean Error MIP – Mixed Integer Programming MILP – Mixed Integer Linear Programming MLP – Multi layer perceptron MOS – Model Output Statistics MPSO – Modified PSO MSE – Mean Squared Error MRI – Meteo Risk Index NO – Normal Operation NOC – Normal Operation Cost NMAE – Normalized Mean Absolute Error NMAPE – Normalized Mean Absolute Percentage Error NPRI – Normalised Prediction Risk Index NPSO – New Particle Swarm Optimization NPSO – LRS - New Particle Swarm Optimization with Local Random Search NRMSE – Normalized Root Mean Squared Error NW – Nadaraya-Watson NWP – Numerical Weather Prediction ORR – Outage Replacement Rate PE – Percentage Error PSO – Particle Swarm Optimization xxi PSO-CEP – PSO embebbed in CEP PSO-LRS – Particle Swarm Optimization with Local Random Search QR – Quantile Regression RBF – Radial Based Functions RCED – Reliability Constraint Economic Dispatch RCUC – Reliability Constraint Unit Commitment RES – Renewable Energy Sources RMSE – Root Mean Squared Error SCADA - Supervisory Control And Data Acquisition SCED – Security Constrained Economic Dispatch SFEP – Swarm direction Fast Evolutionary Programming SF – Skill Forecasting SCO – Sensing Cloud Optimization SCUC – Security Constrains Unit Commitment SHPP – Small Hydro Power Plant SONARX – Self-Organizing Nonlinear Auto-Regressive model with eXogenous input SSCUC – Stochastic Security Constrained Unit Commitment s.t. – subject to UC – Unit commitment WPF – Wind power forecast WPP – Wind power production WS – Wait-and-See W2P – Wind to Power WRF – Weather Research Forecasting xxii xxiii List of variables A Matrix of coefficients in linear programming A Incremental response of the power generation to the precipitation a i Coefficient of cost function of thermal unit i a ik Coefficient of cost function of thermal unit i with fuel k A i Non-load cost, thermal unit i AS i,t Boiler cool down coefficient of unit i at time t B Hydro generation decay in dry days b i Coefficient of cost function of thermal unit i b ik Coefficient of cost function of thermal unit i with fuel k b Vector of coefficients in linear programming formulation BM i,t Base value of maintenance cost of unit i at time t BS i,t Boiler start-up cost of unit i at time t C Vector of coefficients in linear programming formulation c i Coefficient of cost function of thermal unit i c ik Coefficient of cost function of thermal unit i with fuel k ( ) d i Ck Shut-down cost, thermal unit i, period k ( ) p i Ck Production cost, thermal unit i, period k ( ) ,ps i Ck Production cost, thermal unit i, period k, scenario s ( ) u i Ck Start-up cost, thermal unit i, period k cw j (k) Wind generation shed, wind unit j, period k C BO Blackout cost CC i Cold starting cost, thermal unit i C GENSET GENSET cost C WC,h Cost of wind curtailment, hour h C LS,h Cost of load shed, hour h C min|WC,h Cost of production below the generators minimum, hour h D Recourse matrix of stochastic programming D(k) Load, period k f t D Load at time t D i,t Number of hours that the unit i is off-line DR i Decrement ramp rate of unit i e i Coefficient of cost function of thermal unit i e ik Coefficient of cost function of thermal unit i with fuel k Chapter 1 - Introduction 2 1. Introduction Power systems have been changing significantly during the last decades and they will keep changing in the near future. Changes occur due to several reasons, namely: environmental obligations, security of supply, new generation technologies, technological development mainly in communications and control, and the need of new market opportunities and deregulation. Although in real life the forecasting procedures imply some uncertainty around the load and renewable production, such as wind/solar forecasts (caused by forecast errors and variability), up to a very recent past, a single expected value for each forecasting horizon, called deterministic, spot or point forecast are provided and were used in the generation dispatch and commitment procedures [1]. As load and generation can deviate from their forecasts, it becomes increasingly unclear (particularly, with the increasing penetration of renewable resources) if the system will be able to meet the conventional generation requirements within the look-ahead horizon, such as reserves, ramps, minimum up times, minimum down times and power balance. Additional balancing efforts are needed as it gets closer to the real time and additional costs will be incurred by those needs. Although there are several works developing the incorporation of these uncertainties into power system operations, large proportion of these efforts are limited to wind generation uncertainties and ignore the fact that there are additional sources of uncertainty, such as other renewable and intermittent generation and unexpected generation or transmission lines outages. Increasingintroduction of electric energy production with renewable sources, mainly those with high variability, has created several challenges to the market operator and/or energy networks operators, predominately in the scheduling chapter. Both have to do their scheduling, with different detail levels, depending on technical or economic objectives in order to minimize the operation cost, taking into account the problem restrictions: reserve managing, reliability guarantee, management of consumption and management of independent power producers. These optimization problems take the form of unit Commitment (UC) and Economic Dispatch (ED) problems. When speaking about UC and ED in a generic way it must be considered that there are different levels such as Security Constrains Unit Commitment (SCUC), Reliability Constraint Unit Commitment (RCUC), Security Constrained Economic Dispatch (SCED) and Reliability Constraint Economic Dispatch (RCED). This problem is enhanced in low power networks, especially in islands without any connection to continental networks. Due to its large implementation, generally, wind generation is considered the main source of variability on renewable generation but there are other sources that can introduce much faster variations, as solar generation [2]. On the other Chapter 1 - Introduction 3 hand, hydro generation, with easier operation control, may depend on economic strategies, disconnected from available resource. If there is a small storage capacity, the production can be temporally disconnected from the rainfall. Still in island context, a large variation on renewable production can lead to stability problems in the network, which can create generation and/or load shed and, at limit, blackouts. This way, and for security, the scheduling is generally done by a conservative way, with low risk but sometimes far away from an optimal operation. Consequently, reserve levels are greater and the possibility of wasting renewable production leads to a more expensive operation. An efficient use of accurate short-term probabilistic forecast of renewable energy sources (RES) and load can allow optimal committed/dispatched thermal generators, avoiding generation, load shed or even blackouts. Motivated by the importance of these topics, there are a number of research works related to this area, [1],[3]–[16] , among many other cited throughout the text: • Characterisation of spatial and temporal uncertainties related to renewable sources (wind, hydro and solar):I is necessary to understand the variability and uncertainty associated with those resources, considering different levels of temporal and spatial resolution and for different levels of aggregations. It is also required the evaluation of optimal trade-off between cost and details of information used in these new scheduling approaches; • Forecast models for wind, hydro and solar generation: There is a considering number of research works for wind forecast but a better comprehension of characteristics (such as: spatial resolution, temporal granularity, time horizon, optimal forecast refreshment, error and information requirements) is needed. These studies are necessary to characterise the best configurations of the forecasts to be used in the scheduling. • Performance of the actual forecasts services in use: Improvement of studies of real time series of forecasts, characterizing the accuracy of the different forecast types (wind, solar, small-hydro and consumption) is needed as well as evaluation of the value of these different forecasts for the system. • Modelling and processing of uncertainties related with forecasts: There is a portion of researches in this area but difficulties in the integration of this information on the decision process associated with the scheduling are an issue • Scheduling specifications in market environment approaches: Scheduling resulting from the different market sessions is itself a specification of the initial and base Chapter 1 - Introduction 4 solution for the scheduling optimization. It is necessary the integration of this specification on the scheduling optimization structure. • Scheduling optimisation solutions must be improved in order to increase velocity, precision, temporal resolution and time horizon. Work on new structures of the scheduling process which use several parallel scheduling runs, using different resolutions and horizons in order to use incremental information procedures is a necessity. These approaches will improve the performance for realistic scheduling operations. • Faster and more precise scheduling optimization algorithms with the aim of solving larger problems with high number of variables. 1.1 Selected research questions There are several questions that should be answered before the formalization of the final case study. The research must be split in smaller problems, which will contribute to get the global problem solution. 1.1.1 Which is the most adequate probabilistic forecasting model to each forecast horizon? Due to the increment of the renewable power sources in the electric systems, there is a remarkable concern about the forecasting of renewable production. Firstly, the efforts were focused in spot forecasts which evolved to models incorporating uncertainties. Nowadays, researchers try to improve models, pursuing forecasts which will better capture uncertainties related to the resource. There are researches in areas such as solar [17], small hydro, wind or even waves [18] but stronger efforts are focused on wind generation. Nevertheless, the uncertainty models that are used in the wind forecast area can be used, subjected to their own specifications, to other renewable resources forecasts. Hypothesis Considering the above, in this work, the studies related with wind forecast will serve as base to remaining forecasts. The Numerical Weather Prediction (NWP) is the main source for forecast uncertainty and, secondarily, amplification and damping effect of the non-linear relation between wind speed and wind power [1]. Following the same author, there are several approaches to model these uncertainties by probabilistic techniques. Following a parametric approach, as demonstrated in [19] , through the 4º moment (kurtosis), the shape of parametric probability density functions (pdf) of forecast errors changes with the forecast horizons. It is Chapter 1 - Introduction 5 also shown that a Beta distribution can be used to model wind power generation errors [19]. However, errors are not always following a Beta distribution [20],[21]. In [22] a study comparing the Beta distribution with the Extreme Value distribution is presented and the authors concluded that Beta distribution outperforms the Extreme Value distribution to low and medium ranges of wind power, being the Extreme Value more adequate to high values of wind power. Aiming wind speed forecasts, several distributions such as Weibull, Rayleigh, Beta, Gama and Normal were tested to model the uncertainty [23]. On the other hand several authors such as [9],[13],[24] and [25], among others, prefer a non-parametric approach, not having to take any decision in advance regarding which distribution should be chosen. This way, the understanding of the probabilistic models and the selection of one sufficiently robust, to represent the uncertainties in a satisfactory way, is a necessity. 1.1.2 What is the aggregation role of RES for the scheduling process? Several studies have demonstrated the decrease of forecast errors and power variability with the renewable power sources aggregation. However, due to its large implementation, the spatial wind power aggregation has been the most studied [26]–[28]. Some studies concerning not only wind aggregation but also other renewable resources were done but generally in small scale or in hybrid solutions [29],[30]. Hypothesis In this work, a mix of renewable power sources aggregation will be studied, with the objective of understanding its impact on variability and forecast errors. It is expected that this reduction should depend on the aggregated power sources and the revelation of each one. In the scheduling process, on a renewable context, aggregation of renewable energy sources (RES) is typically subtracted from load (Load-RES) in order to produce the net load which then can be used to compute the thermal generator requirements [31],[32]. Considering that the forecasts are independent from each other, aggregations will be done by the convolution. 1.1.3 Which optimization tools should be used? Unit commitment and ED are the base of the power systems scheduling processes. Although being solved at the same time, the UC intends to define which and how long the generation units should be online and the ED proposes to define the production of each on-line unit to meet the load at the minimum operation costs. Improvements on unit generation scheduling can lead to significant cost savings, simultaneously ensuring operational restrictions are not violated. The dynamic economical dispatch (DED) is an extension of conventional ED problem, taking into account the ramp rate limits of generating units [33]. Chapter 1 - Introduction 6 Traditional approach to the ED problem considers, for simplicity, that the cost function for each unit is approximately represented by a single quadratic function. The essential assumption is that the incremental cost curves were monotonically, increasing piecewiselinear function, which could be solved by conventional programming methods and optimization techniques, as the Gradient method, Lagrangian function, Lambda-iteration method, base point and participation factor methods, dynamic programming, Newton’s method Linear and Quadratic Programming and Interior Point method, among others [34]–[36]. Generally, these mathematical methods require the derivative information of cost functions. However, the generation cost functions of recent thermal units are not continuous, not convex, neither differentiable due to valve-point loading effect, multi-fuel burn systems and operational prohibited zones. Thus, the ED problem becomes a non-convex optimization problem with constrains, which cannot be solved directly by some of the traditional mathematical methods. A deep survey regarding these issues can be found in [37]. Regarding UC, generally its formulation contains binary variables, defining which units must be online. The traditional algorithms demand the usage of integer programming as mixed integer linear/non-linear programming which can be time consuming. In the case of scheduling with uncertainty, when the uncertainty is modeled by scenarios, it is mandatory to solve each one independently multiplying the computation time. Hypothesis With respect to units cost functions with non-convex, the problem solution becomes much more complex. To overcome this problem, over the past years, metaheuristic tools have gained more and more importance in optimization problems as unit commitment and economic dispatch. Heuristic usually refers to a procedure that seeks an optimum solution but does not guarantee that it will find it or even if it exists. Metaheuristics are general frameworks for heuristics in solving hard problems. Meta-heuristics do not stop in the first local optimum as a simple heuristic does and can be classified into two groups: those performing a single walk using special procedures, trying not to be trapped in a local optimum and those performing multiple walks. As a result of this, several heuristic methods were proposed to solve this kind of problems, such as Genetic Algorithms (GA) [34], Simulated Annealing, Taboo Search, Evolutionary Programming (EP), Evolutionary Strategies, Particle Swarm Optimization (PSO), Bacteria Foraging Optimization [38],[39], Ant Colony Optimization [33],[38],[39], and Artificial Neural Networks (ANN) approach with Hopfield Networks and hybrid artificial intelligence methods [40]. From the base algorithms several improved approaches and hybrid were proposed, as Improved Taboo Search [34], Fast Evolutionary Programming and Improved Fast Evolutionary Programming [41], Improved Chapter 1 - Introduction 7 Particle Swarm Optimization [40] and hybrids as PSO with Evolutionary programming [42], PSO with crossover operations [40], Fast Evolutionary programming with Swarm Direction [33], New Particle Swarm Optimization with Local Random Search(NPSO-LRS) [40] and Real Coded Genetic Algorithm – Ant Colony Optimization [43], among many others [44]. Understanding which optimization model should be used is intended at the end of this work. Concerning the scheduling, it can be accelerated by decoupling the UC and ED in order to avoid the necessity of run the ED to each set of feasible solutions of each hour of each scenario. 1.2 Challenges This work will be developed under real context and with real information. One of the challenges will be dealing with the quality of information. The data set will incorporate measured values of load, wind, hydro and geothermal production, these values can be “polluted” by wrong measures, absence of measured values, unexpected production profiles provoked by outages, malfunctions, amongst others. Also, there is the natural dynamic of power systems with outages due to maintenance, which change the profile of production systems and can skew the datasets used in forecasts. In some of these situations it will be necessary to do some pre-processing to some data set values in order to prevent deviations. On the scheduling context, results will be compared with those obtained by the system operator, implying the discovery of a comparing platform. The dataset is composed by hourly average values, consequently it is not possible to analyse some fast dynamics inside the hour. 1.3 Chosen methodologies Generally, the scheduling process has complex formulations. Due to the dimension and required precision, it becomes extensive and with a big computational effort. Considering this, all the optimization process can be accelerated only using the sufficient calculus resolution to the existing uncertainty level in each moment. Accordingly, it is verified that the information about the uncertainty is useful to accelerate the calculation procedure as well as the forecast information, which is also essential for the scheduling process. At the end it is expected that the process can be accelerated with the new optimization algorithm. The chosen methodology for solving the problem and verifying the hypothesis will be: • Test the forecasting techniques to be applied in the scheduling; • Evaluate the uncertainty to be applied in the scheduling; • Develop adequate optimization methods to be applied in the scheduling; • Evaluate performances of the optimization methods; Chapter 1 - Introduction 8 • Evaluate the benefits of new scheduling method. 1.4 Thesis objective The main objective of this thesis is to develop an integrated set of mathematic techniques aiming the optimisation of electric energy systems operation with significant penetration of high variability resources. In particular, it is intended to specify the characteristics of forecasting systems which best suit the purpose, namely defining types of input data, mathematical models to include uncertainty and the representation of this uncertainty. A new generation scheduling process based on risk assessment for fast solving the scheduling with convex or non-convex cost function is also an objective. With this thesis, the contribution for the mitigation of the economic and environmental impact from the uncertainty related to the operation variables, which affect the scheduling, is expected to be significant. The final results of this research work will be an approach of advanced scheduling integrating: • Forecast of renewable sources, based on state-of-the-art models; • Uncertainties associated with forecasts; • Innovative optimization approaches; • Scheduling optimization tools. The final result should be a complete method built to give a full answer to the energy scheduling systems with high integration of variable power resources. An application will be developed to study the management of the information from a power production perspective. The innovation of this thesis is spread in several parts of the problem but the main added value will be the aggregation of all the components in the global scheduling approach. 1.5 Thesis outline The research work developed within the scope of this thesis is structured in 6 chapters: Chapter 1 In this chapter it is outlined the motivation and conceptual lines as well as the main research questions and challenges. Some hypotheses for solving the research questions and chosen methodologies are described. Finally the thesis objectives are outlined. Chapter 2 This chapter begins with an overview concerning the power generation scheduling, first under a classical deterministic approach and posteriorly under uncertainty. Chapter 1 - Introduction 9 The mathematical formulations are shown and the importance of the reserves is highlighted. It is also presented a brief overview concerning the latest researches in this area. The second part consists of several aspects concerning the stochastic optimization with the manifold approaches, namely recourse problems, distribution problems, chance-constraints problem and worst-case constraints problems are represented. In the third part of the chapter it is done an overview about renewable power production forecast, namely wind, hydro and solar. Lastly, some techniques for power forecasts uncertainty estimation are presented. Chapter 3 This chapter presents an introduction to the case study and an explanation of the necessity of accurate forecasts. In the second part all the steps for achieving renewable and load forecast with associated uncertainty are described. It is explained how the prediction models are chosen, as well as the choice of explanatory variables and some evaluation criteria. Finally the results from wind, hydro and geothermal power forecasts coupled with the load and their aggregation are presented. Chapter 4 This chapter presents the full methodology for the generation scheduling under uncertainty, based on risk assessment. In the first part it is described the case study thermal generation characterization and it is introduced the concept of equivalent optimal generation unit. Next, it is announced an original metaheuristic based on cloud of particles in order to solve non-convex problems. Finally it is presented the scheduling formulation based on risk assessment. Chapter 5 - In this chapter the case studies centered on the proposed methodology are presented. A net load forecasting is done followed by a single-period and a multi-period unit commitment based on risk assessment. The São Miguel Island’s system operator proposed scheduling is presented and the results are compared with those from the proposed methodology. Chapter 6 - In this chapter the overall conclusions are drawn and the original contribution is presented. Conclusions concerning the research questions and formulated hypothesis are addressed and some future research perspectives are also discussed. Annex - In annex are presented some intermediate results which, due to their length, were not included in remain chapters. Chapter 1 - Introduction 10 Equation Chapter 2 Section 1 CHAPTER 2 Scheduling – General Overview Contents This chapter is a general overview concerning the schedule process regarding the unit commitment problem and economic dispatch under deterministic and stochastic approach. It is also an outline of stochastic programming methods. Finally, it is done an overview concerning the renewable power forecast and uncertainty modelling Chapter 2 - Scheduling – General Overview 18 The tertiary reserve has as main objective to guarantee the regulation control, changing the generation by adjusting the scheduling, whose process should be done between 15 minutes and 1 hour after the contingency. This additional reserve is calculated by adding the variation in load with the variation in variable generation. Although being not so important as the primary and secondary reserves, as the contingency reserve must be maintained permanently, the regulating reserve is an additional reserve with additional cost [45]. Traditionally, it is not common to use renewable power sources as reserves, mainly as primary and secondary reserves, since there is resource waste, but it can give a support to the tertiary reserve. There are several approaches to define the amount of reserves. In a deterministic point of view it can be defined as a given percentage of forecasted load, a percentage of renewable production forecasts, equal to the amount of the most loaded unit [47],[52] (generally for the primary reserve) or a mix of some of the previous rules [52]. In the case of probabilistic approach, it uses a function of the probability of not having enough generation to meet the load (due to load and production forecasting errors), or also using approaches based on standard deviation of load forecasting errors, among others. 2.2.2 Economic dispatch The economical dispatch problem is one other important issue in the power system scheduling. Fundamentally, it is intended to evaluate the value that each on-line unit should generate with the lowest cost respecting the technical and load constraints. The ED uses as a basis the UC solution, excluding from the optimization the generation units that are offline. Opposing the UC, which solutions can result from constant costs (in the case of simple ranking of priority), ED characteristics’ of production costs can be nonlinear. Consequently, the optimum is an allocation of generation between the units. On the majority of publications that strictly analyze the subject of ED, the online generation units are already known and commonly the case studies are done using a single period (sometimes to present new solution algorithms). The multi-period analysis is less addressed since it presents a more difficult solution. In reality, ED is done in a dynamic way (Dynamic Economical Dispatch) [33],[53] as it takes into account the variation of demand over time, as shown in equation (2.13), where FCi,t is the cost function of each unit i during interval T, and NG represents the number of on-line units [39]. ( ) ,, 1 min G N T it it ti FC P = ∑∑ (2.13) Chapter 2 - Scheduling – General Overview 19 The basic formulation includes the production limits of each unit as shown in (2.5) or (2.9) and balance equations as (2.8) or (2.12). As for the UC, the ED can include several kinds of constraints, such as SCED which deal with reserves and network constraints and RCED, including constraints about generation and network reliability. Merging the problem of UC with the problem of ED using Mixed Integer Programming optimization algorithms is a widely adopted approach. One of the advantages is the fact that the non-linear details in the ED can justify a change in the solution of the UC. On the other hand, the ED can integrate UC costs considering the fix costs of the generation [45]. Over the past decades, many methods have been developed to solve the ED problem. There are the traditional methods such as Gradient, Lagrangean, Lambda-iteration, Dynamic Programming, Newton’s, Linear Programming and Interior Point, among other methods [34]. Some of these methods are used under the assumption that the thermal unit’s costs functions are convex, continuous, differentiable and monotonically increasing, along the domain. Real thermal units can present cost functions with different characteristics from those mentioned above. For instance, steam turbines can present valve-point effects and some units’ burn different type of fuels resulting in a different function for each fuel. In this case the cost function can take the form equation (2.14), for k types of fuel. ( ) ( ) ( ) ( ) 2 min min 1 , 1 , 1 1 1 1, 1, 1 1 2 min min 2 , 2 , 2 2 2 2, 2, 1 2 , 2 min ,, , sin , for fuel 1, sin , for fuel 2, sin i it i it i i i i t i t i i i i it i it i i i i t i t i i i it ik i t ik i t ik ik ik ik t i aP bP c e f P P P P P aP bP c e f P P P P P FC aP bP c e f P P + ++× × − ≤≤ + + + × × − <≤ = + ++× × −  ( ) ( ) max ,1 , for fuel , k t ik i i kP P P −       <≤   (2.14) Due to possible vibration in the shafts or problems with the continuous start and stop of the coal mills, there are prohibited operation zones were the thermal units cannot work at steady state. Taking this into account, the constraints can be non-continuous as shown in equation (2.15), min , ,1 , ,1 , , max ,, 2,3,.., pz LB i it i UB LB it i j it i j pz LB iN it i P PP P P PP j N P PP − ≤≤  ∈ ≤≤ =  ≤≤  (2.15) where ,LB ij P are the lower bound of the jth prohibited zone of unit i and ,1 UB ij P− is the upper bound of the (j-1)th prohibited zone of the same unit. Thus, the ED problem becomes a non-convex optimization problem with constraints, which cannot be solved directly by some of the traditional mathematical methods. Dynamic programming can solve this kind of problem, but Chapter 2 - Scheduling – General Overview 20 can suffer with the dimension and the time needed to solve it [35],[36] and [54]. On the other hand, commercial tools, which are able to solve economical dispatch for thermal units, always require convex cost functions. This condition can be attributed to the limitations of the optimizing tool or the need of rapidity and non-convex algorithms tend to be slow. To overcome this problem, sometimes the technique is to split the space solution in convex subspaces and then use conventional algorithms. This technique may create a vast number of solutions, some possible, others not, and the best solution must be found inside the set of feasible results. Alongside, several heuristic methods were proposed to solve this kind of problem such as Genetic Algorithms, Simulated Annealing, Taboo Search, Evolutionary Programming, Evolutionary Strategies, Particle Swarm Optimization, Bacteria Foraging Optimization, Ant Colony Optimization, Artificial Neural Networks approach with Hopfield Networks, and hybrid artificial intelligence methods. From the base algorithms several improved approaches and hybrid were proposed, as Improved Taboo Search [35], Fast Evolutionary Programming and Improved Fast Evolutionary Programming [42], Improved Particle Swarm Optimization [55] and hybrids as PSO with Evolutionary programming [54], PSO with crossover operations [38] and Fast Evolutionary programming with Swarm Direction [39], among many others. 2.2.3 Unit commitment and economic dispatch under uncertainty In a classical approach (deterministic), without the integration of renewable source power plants with their intermittent and variability profile, the source of uncertainties is only related to the load forecast or some unexpected unit or line outage. As a consequence, security of a power system refers to its ability to survive to contingencies, while avoiding any undesirable disruption of service. As a security measure, the so called N-1 security criterion is commonly used, where the system is considered to be N-1 secure if any single component outage does not lead to an overloaded component or to other operational violations. When renewable power sources are included, the amount of uncertainty increases, due to the variability and errors introduced by the forecast process. As a consequence scheduling becomes a much more challenging problem. However, even when considering renewable power sources, if a spot forecast (point forecast) is assumed, the problem can be considered as deterministic too. The difference between uncertainty and variability is a matter that must be taken into consideration; the variability is related to the type of resource, which can, by its nature, change without any capacity of management like in wind or solar sources. While uncertainty is linked with the forecasting, which, considering the nature of the resource, can be hard to forecast with a appropriate certainty. The scheduling process has to deal, in distinct ways, Chapter 2 - Scheduling – General Overview 21 with both cases. Regarding variability the problem impact can be reduced with fast scheduling processes, giving solutions in advance. Forecasting errors, characterized as the difference between the values used to the unit commitment and the real ones available in the real time dispatch, can cause serious difficulties for the system operator, who must balance the positive or negative deviations of an intermittent production. To deal with uncertainties, various stochastic analyses have been developed. In a stochastic approach, the main question, cited from [4], “What is the level of operating reserves that should be imposed at the UC stages, to take into account the additional uncertainty from the renewable power sources ?” must be addressed. Nowadays, even when considering the increase of solar and small hydro capacity, wind resource is clearly the most important issue for scheduling systems with high penetration of renewable energy. Correspondingly, the vast majority of publications are focused on wind power generation. Works addressing this matter focus on three main concerns: • How to improve forecasting techniques; • How to model the uncertainty and its impact on the reserves and consequently on the scheduling; • How to make the scheduling process faster with less computational burdens. There are many models concerning UC, which fundamentally differ from how the constraints are formulated to capture the dynamic performances of generation units or in the cost models. The overview below gives a general idea about the more recent research lines regarding these themes. In [56], a comparison between different probabilistic forecasting and scenario reduction methods to solve the UC and assess operating reserves is done. The scenarios are generated from probabilistic density functions created by Quantiles Regression (QR) and Kernel Density Estimators (KDE) (based on Nadaraya-Watson (NW) estimator), using the wind power forecast as explanatory variable. The results from stochastic and deterministic UC are compared. It is concluded that a higher number of scenarios improve the performance of the stochastic UC strategy in spite of increasing several times the computational efforts. Their case study showed that the random scenario reduction is in line with other techniques like those resulting from Kantorovich distance. It is also concluded that the dynamic reserve derived from NW outperforms the one resulting from QR. The scenarios created from NW outperform the scenarios created with QR, in the case of UC. Similarly, in [16], and in [4], the probabilistic method was based in scenarios where the focus was the impact of wind power uncertainty on power systems operation as UC, ED, and Chapter 2 - Scheduling – General Overview 22 reserves management. It is done a comparison between a stochastic approach and a deterministic approach with different levels of reserves and conclusions are that the deterministic formulation from a certain value of reserves requirement onwards, reaches results comparable with the stochastic approach. It is also concluded that wind power forecast errors have a great impact on the scheduling of generation units in a day-ahead market with implications on the real-time economic dispatch. They do not address other uncertainties, such as load uncertainty and transmission line or generator outages (It is considered that the reserves based on classic criteria are sufficient to deal with those uncertainties). In [15] it is presented a deterministic formulation for the UC and its extension to stochastic programming formulation. Wind generation is considered as the source of uncertainty, where the wind speed uncertainty is estimated with the use of ensemble approach, while the load has no uncertainty. The network model and forced outages are not considered. The authors sustain that when using stochastic formulation to represent the wind uncertainty, the requisite of a previous value of reserves can be strongly reduced. In fact, a comparison of robustness is done between an explicit value of reserves and the implicit amount of reserves that result from the stochastic formulation. In [57] an SCUC with wind power generation is analyzed to test how the power system reacts by redispatching thermal units in real time when the actual power is different than the forecasted. It is pointed out that the ramping of the thermal units is crucial for accommodating the uncertainty and variability of wind power generation. In this case the wind uncertainty was formulated with generated scenarios, including a Latin Hypercube Sampling (LHS) technique in the simple Monte Carlo simulation. Concerning the uncertainty, the assumption is wind power errors follow a normal distribution and the mathematical formulation does not include load shed neither wind curtailment. Transmission network constraints are also included in the problem. In figure 2.2 the flowchart of proposed SCUC is shown. The SCUC is solved as a master problem with the wind power generation spot forecasting and Figure 2.2 –Security-constrained unit commitment with wind power Chapter 2 - Scheduling – General Overview 23 a second problem is solved simulating scenarios to represent the wind power uncertainty, both solved by Mixed Integer Programming (MIP). The large scale mixed integer UC is solved as a master problem followed by a network security check sub-problem. In case of any violation of the sub-problem, the Benders cut is formed and added to the master problem. As in [57], paper [58] presents a combined model of optimal reserve ED considering uncertainty of wind power generation. The SCUC and reserve dispatch problem based on forecasted wind power is settled in the master problem. In the sub-problems the obtained generation dispatch has to meet the requirements due to the uncertainty of wind power. The volatility of wind power was simulated by scenarios generated by LHS technique in the Monte Carlo simulation (based on normal distribution). The difference is the introduction of wind curtailment and load shed into the original formulation of [57]. It was verified that incrementing wind curtailment costs, the amount of wind energy scheduled increases and, naturally, the cost of reserves increases too. Reference [59] presents a full-scenario security-constrained unit commitment (FS-SCUC) when considering wind generation and load variability. Following the line of [57], the Benders decomposition of the whole problem into a main problem and two subproblems is proposed, aiming the reduction of the scale of UC and improve the computational speed using Mixed Integer Linear Programming (MILP). The main problem solves the UC without network security constraints. Subsequently, the first sub-problem checks whether the commitment and dispatch solution of the master problem can satisfy the network security constraints or not. The second sub-problem checks if the worst-case security constraints can be satisfied with the obtained schedule on the master problem. In the formulation, active and reactive powers are incorporated. The wind power generation and node load are treated as volatile loads (the wind power is considered as a negative load) and the uncertainty is modeled by forecasting intervals which will originate scenarios. Reserves are implicitly determined by a specific commitment schedule corresponding to the full-scenario security constraints. As in [57] and [58] it is concluded that cost decreases if wind power curtailment is allowed and the curtailment compensation is small. In figure 2.3 a flowchart of overall FS-SCUC is represented. Figure 2.3 – Overall procedure of SCUC Chapter 2 - Scheduling – General Overview 24 In [60] it is proposed a Constrained Ordinal Optimization (COO) method for solving a scenario-based Stochastic Security Constrained Unit Commitment (S-SCUC). Consideration of Ordinal Optimization aims the facilitation of the scenario-based solution method, since its goal is to seek good enough solutions with high probability instead of searching the best solutions with absolute certainty. In this specific case, it is used a Weibull distribution for the wind speed, which is converted into power by a wind turbine power curve, and load forecast errors are modeled by normally distributed functions with constant standard deviation. Random outages of generation units and transmission lines are both considered. The wind curtailment is allowed and load shedding is also considered but only as the last resort to maintain the system reliability. As in previous methods, the corresponding Security Constrained Economic Dispatch (SCED) is repeated for all scenarios until a feasible UC is obtained. In figure 2.4 the flowchart of a feasibility model is shown. Generally, the final solution may not be the optimal for the stochastic SCUC problem but the purpose of the proposed COO method is to find a good enough solution with a high probability instead of an optimal solution. As conclusion, the proposed COO demonstrated to be computationally more efficient than an MILPbased approach, showing a remarkable decrease of computational time. Commonly, in many approaches, it is required that the load must be fed, even in the worst cases, leading to high levels of reserve. In [31] a Chance-Constrained Unit Commitment is proposed, when considering that the demand should be met in any plausible scenario, leaving out extreme scenarios in order to reduce the cost of dispatched energy. The chance constraint problem is formulated subjected to stochastic demand for the scheduling of energy and spinning reserves with a α-quantile and n-K criteria. The proposed chance constraints are the risk measure of not meeting the demand with a certain confidence level α. As in [59], wind generation is considered as a negative load creating the variable net load, resulting from the subtraction of wind generation from load, both defined by normal distributions with fixed standard deviations. The scenarios of net load are generated by a LHS technique based on the normal distribution of net load. In a single period case, it is evaluated the total cost versus wind energy penetration. Figure 2.4 – Flowchart of the feasibility model Chapter 2 - Scheduling – General Overview 25 On the formulation, some probability of generation outages is also considered but the network model is not considered. The formulation is done in two stages: firstly it solves the UC with up/down reserve scheduling, and secondly the real demand has to be met under a certain probability. There are reserves to deploy the real demand, but load shedding or wind curtailment are allowed and the chance constraints are defined as reliability constraints of system, limiting the load shedding and wind curtailment. In [61] a Chance-Constrained Unit Commitment is also proposed. The problem is formulated as a chance-constrained two-stage stochastic programming and a combined sample average approximation (SAA) algorithm was developed to solve the model efficiently The idea behind SAA is to approximate the real distribution of random variables by a Monte Carlo sampling empirical distribution. The authors state that the model ensures, with high probability, that a large portion of the wind power output at each operation time is used. Wind power is the only source of uncertainty since the load is considered deterministic. The wind power uncertainty is represented by scenarios obtained from a normal distribution, with standard deviation equal to a percentage of the expected values. Compared with [31] the network is also considered in the model but the outages are omitted. In [62] a comparison between a scenario-based and interval optimization approaches to stochastic SCUC is done. The scenario based stochastic SCUC problem is decomposed into a main problem and three sub problems. The main problem is an UC without network security constraints and serves as base case. Then a first sub problem addresses the hourly network evaluation of the main problem solutions for the base case. A second sub problem does the hourly feasibility check for each scenario, checking possible violations of the main UC solution in each scenario. Finally, the third sub problem checks the optimality of master UC solution in each scenario. In figure 2.5 a scheme of the proposed scenario-based approach is shown. Figure 2.5 – Scenario-based approach The reserve requirements are implicitly represented by deviations in the dispatch solutions of the base case and scenarios and are optimally determined via preventive and corrective actions. Chapter 2 - Scheduling – General Overview 26 Different from the scenario-based approach, the interval optimization doesn’t need any preknowledge of wind generation probability distribution. In this case only the wind power forecast has uncertainty. The interval optimization requires less computational efforts to generate upper and down bounds to the objective value. However, the choice of the uncertain interval coverage rate can dramatically change the optimal solutions. It is also not appropriated to simulate discrete random variables as units or transmission lines outages. The wind speed scenarios generation are based on a Weibull distribution and converted to power by a turbine power curve. In the case of intervals, the uncertainty is considered as a percentage of installed capacity. The authors conclude that the scenario-based approach provides a more suitable solution. Yet, it can be a large scale problem with high computation burdens. Although the scenario case needs more computational efforts, it is less sensitive to the starting values. On the other hand the uncertainty interval reveals to be faster than scenario-based, being equivalent to two scenarios, but is very sensitive to the assumed uncertainty interval. In [32] it is stated that traditional UC with deterministic spinning reserve requirements are inadequate, given the variability and uncertainty of wind power. Thus, a new probabilistic model of SCUC is proposed to minimize the energy cost, spinning reserve and load shedding (Expected Energy Not Served (EENS)), using the determination of additional spinning reserve due to the integration of wind generation. The formulation of EENS takes into account the probability distribution of forecast errors of wind and load, and the outage replacement rate (ORR) of several generation units. There is not any assumption about spinning reserve constraint because its value is based on an internal cost/benefit analysis. The net load (load minus wind power) forecast errors are assumed to be normally distributed with zero mean and a fixed standard deviation. The proposed approach determines the optimal amount of spinning reserve, which minimizes the total cost of system operation (balancing energy costs, start-up cost, reserve and expected cost of load shedding), reaching the trade-off between economy and reliability of the system. It is concluded that, even increasing the amount of EENS, comparing with those obtained by the traditional constant reserve UC, the operating cost is lower as well as the total cost. The reserves are strongly dependent on value of loss of load and ORR and less dependent from load and wind forecasts errors. In paper [63], unlike the previously mentioned, it is analyzed a different approach.. Pumpedstorage and wind unit generation are coordinated and optimized with a stochastic SCUC model through several coordination strategies. The proposed optimization problem is formulated as a Mixed Integer Problem (MIP). The Benders decomposition technique is used to decompose the original large-scale problem into a more tractable master MIP problem and several Linear Programming (LP) subproblems. The subproblems check the power flow Chapter 2 - Scheduling – General Overview 27 which results from the master problem solution in the base case and all scenarios. If there is any violation, the corresponding feasibility Benders cuts are generated and fed back to the master problem, for the solution of the next iteration. The stochastic SCUC for coordinated scheduling of wind-pumped-storage units is shown in figure 2.6. Forecast errors of wind and load and random outages of generators and transmission lines are taken into consideration. The load forecast errors are represented by a truncated normal distribution with a fixed standard deviation and the average is equal to the hourly power forecast. The wind power forecast is modeled by Auto Regressive Moving Average (ARMA) and the uncertainty is equal to a fixed percentage of the wind power forecast. The scenarios are generated with LHS technique in the Monte Carlo simulation and the forced outage rate of transmission line and generation units are defined by a low percentage value. It is shown that a correct coordination of wind-pumped storage, due to the reduction of variation of thermal generation commitment, may lead to lower total operation costs, wind curtailment and corrective actions in scenarios costs. Figure 2.6 – Stochastic SCUC for coordinated scheduling of wind pumped-storage units A different approach from stochastic programing is proposed in [64], namely the robust scheduling based on Robust Optimization (RO). In this approach, the solution is considered robust if it remains “close” to optimal for all scenarios and robust model if it remains “almost” feasible for the same scenarios. Being so, there are not unfeasible solutions since RO finds the solution which violates the constraints by the last amount. With the extreme Chapter 2 - Scheduling – General Overview 34 The presented model is quite complete, taking into account the power ramp limits, minimum up/down times and the contributions of each unit to the spinning reserve. Clearly, there are other modified formulations that follow the same line, considering the spinning reserve as a percentage of the total production, assign cost to the reserves, allow load shed, among others. The ramp limits have a specific role allowing the study of the response of thermal units in the presence of changes in generation and load or even forced generation outages. Nevertheless, the model can be improved introducing the network model and consequent power flow, evolving to complete security and reliability formulations. 2.2.5 Stochastic scheduling formulation The deterministic formulation presented in 2.2.4, can be generalized to stochastic formulation running several scenarios of wind power or load, instead of only one. The differences from the deterministic formulation are basically in the objective function and in the notation. There is the introduction of the notation s, in the formulation (2.40) up to (2.52), which correspond to the scenario number [4],[15],[46],[47] and [56]. ( ) ( ) ( ) , 1 11 min G N SK s ps u d i ii s ki prob C k C k C k = = =  ++   ∑ ∑∑ (2.40) s.t. ( ) ( ) ( ) 11 WG NN ss ji ji pw k pt k D k = = += ∑∑ (2.41) ( ) ( ) ( ) ( ) 1 . G Nss i i i pt k pt k D k r k =  −≥   ∑ (2.42) ( ) ( ) ( ) ,s s fs jj j pwk cwk PW k+= (2.43) ( ) ( ) ( ) ( ) ,,, 1 . L ps s i i j li li l C k Au k MC k k δ = = + ∑ (2.44) ( ) ( ) ( ) , 1 . L s ss j i i li l pt k P u k k δ = = + ∑ (2.45) ( ) , , sli li k δ ≤∆ (2.46) ( ) , 0 s li k δ ≥ (2.47) ( ) ( ) ( ) . ss iii i PTukpkptk≤≤ (2.48) ( ) ( ) 0. si ii pt k PT u k≤≤ (2.49) Chapter 2 - Scheduling – General Overview 35 ( ) ( ) ( ) ( ) ( ) ( ) 1 . 1 . 1 .1 ss i i i ii i i i i ptk ptk RUuk SU uk uk PT uk≤ −+ −+ − − + −       (2.50) ( ) ( ) ( ) ( ) .1. 1 si i i ii i ptk PTuk SDuk uk≤ ++ − +   (2.51) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 . .1 1 sj i i i ii i i i i pt k pt k pt k RD u k SD u k u k PT u k−− ≤ + + −− + −       (2.52) The minimum up/down time equations (2.34) up to (2.39) do not suffer any changes in the formulation. Further, this scheduling formulation can be extended, considering corrective actions on the thermal generation units taking into account the uncertainty of wind power plants, loads, generation and transmission outages. The problem can be transformed into two (or more) stages stochastic programming as proposed in [15],[57]–[63] and [65]. In the first stage it is intended to decide the value of the thermal power output and the commitment profiles over the entire scheduling horizon. The power outputs are considered as nonanticipatory (here and now) because it is assumed that the forecasted values are known. The formulation of second-stage is done regarding multi realizations of the forecasted variables in order to meet scenario dependent thermal power output (wait-and-see) [15]. In order to keep the problem computationally tractable, the commitment variables are not scenario dependents. If the committed unit’s remains on/offline in the second stage, it means that, in practice, an economic dispatch is being run to each scenario. 2.3 Stochastic optimization In an environment without uncertainty, many deterministic model-based approaches have been developed. Depending on the objectives, decision variables and constraints, the deterministic optimization problems can be formally classified as Linear Programing (LP), Integer programming (IP), Mixed Integer LP (MILP), Non-Linear Programming (NLP), and Mixed Integer NLP (MINLP), among others. These model-based approaches are often impossible to implement in “real world” due to the high modelling complexity or, if there is some kind of uncertainty. This uncertainty can be introduced due to the stochastic nature of objective functions, variables or constraints, especially in the case of dynamic and complex systems. These may have parameters where, usually, the uncertainty is generating large effects on the objective functions and constraints. This is the case of power production scheduling, where there is the necessity to take decisions for the future, depending on the load and power production forecast uncertainty unexpected outages. Hence, it is always present the great challenge of taking decisions future operations. As a consequence,it is necessary to implement models of optimization under uncertainty. Chapter 2 - Scheduling – General Overview 36 The Stochastic Programming or stochastic optimizing problems are a sort of optimization where the stochastic properties of the uncertainties on the data and the model are taken into account. Stochastic programming is a mathematical technique which explicitly incorporates uncertainty of some parameters, underlying the optimization model. It can be used in financial planning, supply chain management, transportation logistics, telecommunications, network designs and energy systems planning, among others [67]. There are two main modelling issues in the stochastic programming, namely the optimal resources allocation model and randomness model, as presented in figure 2.7 [67]. The optimum decision model, together with the constraints, constitutes the core of the problem, which must be solved and it depends on the specific characteristics of each application problem. The model of the randomness and scenario generation is the major issue in the application of stochastic programming, it is the representation of the underlying random process. Figure 2.7 –Breakdown of stochastic programming In figure 2.8 is presented a classification of stochastic programming models which can be applied both to linear and non-linear programming [67]. They are classified considering the way which the uncertainty is defined and how the problem is adapted to the optimization model. Figure 2.8 –Classification of stochastic programming problems Chapter 2 - Scheduling – General Overview 37 To illustrate this kind of problems the analysis can be initiated with the formulation of a linear problem defined by (2.53), where × ∈mn RA , , n xR∈C and ∈ m Rb [67]. min .. 0 x st x x ≥ ≥ C Ab (2.53) Lets also consider a discrete probability space represented by ( ) ,,PΩℑ and the realization of the uncertain parameters denoted by ( ) ξω with ω ∈Ω . For each event ω, the realization of the parameters A, b and C is defined as ( ) ,, ωω ξ =AbC . The probabilities associated with these realizations can be denoted by ( )( ) P ξω or only ( ) P ω [67]. 2.3.1 Distribution problems Distribution problems are broadly known as those which provide the distribution values of cost function for different realizations of the random parameters and also for the expected values of such parameters. The distribution problems can even be split in Expected Value (EV) and Wait and See (WS). Expected value The expected value model is built by changing the random parameters with their expected values as shown in (2.54). This way the model of EV can be considered as a linear problem, in fact, the uncertainty is handled before its introduced into the underlying optimization model. Despite not representing the full distribution and uncertainty, the usage of EV formulation, can be used in order to gain some sensibility regarding the decision problem. ( ) ( ) ( ) ( ) min .. 0 x st E xE x ξω ξω  ≥  ≥ C Ab (2.54) Wait-and-see In wait-and-see problem it is assumed that the decision-makers are capable of delaying their decisions, waiting until an observation is made on the random element and then solve it as a deterministic problem. Therefore, this approach is based on perfect information about the future [67]. In such cases when additional measurement information on the uncertainties becomes available, the operational strategy can be adopted. This strategy requires the solution Chapter 2 - Scheduling – General Overview 38 of several deterministic optimization problems in order to find the deterministic optimal decision at each scenario or random sample. However, wait-and-see strategy does not consider the uncertainty properties and has some drawbacks, such as the actions, which are always taken a posteriori. Furthermore, a feedback control cannot ensure constraints on openloop variables [68]. Since the complete future realizations are rarely known, this approach is not suitable for handling time-varying processes, as power generation scheduling. These models are often used to analyse the probability distribution of the objective values and belong to linear programming models family, being each one connected with an individual scenario [67]. The formulation for WS problem can be depicted as in (2.55). ( ) ( ) ( ) ( ) min .. 0 x st x x ξω ξω ≥ ≥ C Ab (2.55) 2.3.2 Recourse problems (here-and-now problems) The recourse problems are also known as here-and-now problems. While in the distribution problems approach, the expected values of the uncertainties are used in the problem formulation, the here-and-now problems involve the definition of both the objective function and constraints in terms of some probabilistic representation (expected value, variance, and quantiles). Moreover, the decision variables are separated from uncertain parameters. In the recourse formulation, it is allowed the constraints violation, but are penalized through a penalty term in the objective function. (This approach is only recommended when the objective function and constraints are able to be described by the same measurement). If the cost model is hard to model or when the constraints are associated with safety requirements, it is better to not compensate for violations by additional costs. In these cases it is encouraged to maintain a high level of reliability. This means that constraints have to be satisfied at least with a probability exceeding some pre-selected value. As it is a recourse problem it can be solved in several stages. The single stage, solution of the objective function by stochastic programming model can be formulated as (2.56) [67], where xF∈ . ( ) ( ) min Ex FF ω ω ξω ∈Ω   = C  (2.56) Chapter 2 - Scheduling – General Overview 39 The optimal value of (2.56) represents the minimum expected cost of the stochastic problem, and the optimal solution * xF∈ hedges against all possible events of ω that may occur in the future [68]. The standard formulation of the two-stages stochastic programming model with recourse is done by (2.57)[68]. The intention is that the sum of first stage costs and the expected value of the random second stage, or recourse costs, is minimized [69],[70]. The objective is to choose the first stage variables in a manner that the sum with expected value of the second stage cost is minimized. ( ) ( ) min , .. 0 x Qx st x x ξω + ≥ ≥ C Ab (2.57) where, ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , min .. 0 ω ω ξω ξω ξω ξω ξω ξω ξω ω  = = + ≥ ∈Ω Qx E st x qy Dy h B y (2.58) The vectors A and b from (2.57) are known without any uncertainty. The function Q(x, ξ (ω)) represents a non-linear term, which is referred to as the recourse function. The matrix B( ξ (ω)), D( ξ (ω))(recourse), the vectors h( ξ (ω)) and q( ξ (ω)) may be random. For a given first stage decision x, the corresponding recourse actions y( ξ (ω)) are obtained solving the sub-problem associated with recourse function Q(x, ξ (ω)). The future unfolds in several sequential steps and subsequent recourse actions are taken dealing with the generalization of the two-stage recourse problem, known as multistage stochastic programming problem with recourse. A decision made in stage t should take into account all future realizations of the random parameters and such decisions only affect the remaining decisions in stages t+1..T [67],[68]. The general formulation of a multistage recourse problem is set out in equation (2.59), Chapter 2 - Scheduling – General Overview 40 { 12 3 11 2 22 33 32 1... 2 11 1 1 21 1 22 2 2 31 1 32 2 33 3 3 11 22 33 min min min ...... min .. TT TT xx x T T T TT T T tt Cx E Cx E Cx E Cx st x xx xxx xxx x lxu ξξξ ξξ ξ −    + + ++       ≥ ≥ ≥ ≥ ≤≤ Ab AA b AAA b AA A A b  (2.59) where t =1,…,T represents the stage in the planning horizon, and the vectors ξt = (bt,ct,At1,…,AtT) with [ ] 2,...,tT∀∈ are random vectors on a probability space ( ) ,,PΩℑ . This approach already gets better values than those obtained by the expected value approaches of distribution problems. If any feasible solution obtained by the expected value exists, it is already contemplated in the here-and-now model. The two-stage recourse problem was used in works presented in [15],[58] and [63], while in [65] a three-stage problem was used (all referenced in table 2.1). Another approach for obtaining worthy solutions is the sample average approximation (SAA) method [69]. A sample ξ1,…, ξN of N realizations of random vector ξ(ω) is generated and with this, the real distribution ξ is replaced by an empirical distribution corresponding to a Monte Carlo sampling (2.60). This technique was explicitly used in [15] and [31]. ( ) ( ) 1 1 min , .. 0 N i i x Qx N st x x ξω = + ≥ ≥ ∑ C Ab (2.60) 2.3.3 Chance-constraints problem Other method of stochastic programming is the probabilistic or chance-constrained which focuses on the system reliability. The reliability is the system ability to remain feasible in an uncertain environment. It can be expressed as the minimum requisite of the probability of satisfying the system constraints. This optimization technique deals with random processes where one (or several) constraints or an objective function must be satisfied with high probability, as defined in equation (2.61). Therefore there are two main reasons for the problem to become intractable; the formulation of the constraint can be hard and the feasible Chapter 2 - Scheduling – General Overview 41 search space limited by the chance-constraint is generally not convex even if the constraints are convex in x for each realization ξ of ω [31]. In the case of α equal to 1 the problem is equivalent to a deterministic one. ( ) ( ) ( ) ( ) ( ) ( ) { } min .. 0 x st PA b x ξω ξω ξω α ≥≥ ≥ C (2.61) This approach was used in [31] and [61] to solve an UC with wind power uncertainty. 2.3.4 Worst-case constraints To deal with the a priori unknown operating reality, two general methods are widely used, the worst-case and the base-case. While in the base case it is used the nominal (mean) value of the uncertainty variables, the worst-case is a simplified approach for evaluation of the robustness and feasibility. Generally, it is applied to problems where the distribution is unknown. It is assumed that all the variations can occur simultaneously in the worst combination possible. It is a very conservative analysis, since it considers that worst cases of variables or parameters deviations will occur simultaneously. Nevertheless, and although the reachable “low profit”, the worst-case is widely used in the optimization areas due to its simplicity and reliability in ensuring the constraints. This is the case of [59] and [64] where it is defined a set of extreme scenarios for description of wind power uncertainty. After this study some conclusions could be obtained and are presented in table 2.2. Table 2.2 – Overview of stochastic programming Expected value Chance-constraint Worst case Solution not robust Robust solution Absolutely robust solution Solution at low cost Not too expensive Extremely expensive Easy to solve Often difficult to solve Cannot even exist So far it is verified that there are several sources of uncertainty in the power systems scheduling. Some authors include uncertainty related with generation unit’s outages or transmission lines, uncertainty in load forecast and/or in renewable power forecast in their formulations. However, it is a constant presence in different formulations or solutions techniques, the uncertainty related with renewable power production, namely wind power generation. The renewable power production forecast is based on several techniques which may depend on the type of the renewable source, the available information, forecasting Chapter 2 - Scheduling – General Overview 42 horizon, among others. To understand these issues an overview concerning renewable production and load power forecast is done. 2.4 Power forecast methodologies Despite not existing a unique definition for classifying the temporal horizon of forecasting, it is current to consider three temporal scales: very-short-term (up to 6h or even 9h) [71], shortterm (6h to 72h) and mid-term (3 days to 10 days) [71]–[73]: Very short-term forecast: The application of this time horizon depends on the market rules with the forecasts being useful for trading in intraday markets. For the system operator, the usefulness of these forecasts is related to the ancillary services management of the power system, as well as for UC and ED refreshment. It is also useful in the management of rapid conventional power plants (very usual in isolated systems as islands); Short-term forecast: Forecast for a time horizon between 6 and 72 hours. This time horizon is strongly related with the power system scheduling, namely UC and ED. This time horizon is of great importance for the input in the electricity daily market. In this case the forecast horizon is defined by the requirement of market operator. The short-term forecasts can also be used for maintenance scheduling, particularly when the time horizon is 72 h; Mid-term forecast: From 72 hours up to 7 days allowing programming maintenance operations. In an energy system integration level, the power forecast can have the most varied applications, namely [71],[72]: • Optimization of energy grid management at the level of economic dispatch or predispatch of power plants, dynamic security assessment, reserves allocations, power flow with neighbouring systems, water storage in power plants reservoir, etc. The forecasting horizons depend on the production system size and conventional power plants type. • Optimization of energy trade in the energy market environmental. The producers at market, energy suppliers, energy traders and independent producers define generation schemes for the time horizons set by the market rules, usually with 48 hours in advance, suffering penalties for deviations from these plans; Chapter 2 - Scheduling – General Overview 43 • Power plants or transmission lines maintenance planning, for which longer forecasting horizons can reveal interest, though the acuteness of the meteorological predictions strongly decreases from 5 to 7 days. 2.4.1 Reference models Any new model or forecasting method to predict power can only be considered satisfactory if it obtains better results, i.e. smaller errors than the reference considered methods. The simplest method is the persistence where it is considered that the forecasted value to instant t+k is equal to the one measured at instant t (2.62). ˆ t t kt pp + = (2.62) Despite of its simplicity, this method is hard to beat in very short-term horizons [71]–[74]. The generalization of persistence leads to the moving average method, which provides a future value with the average of n past values. At the limit, the average is all available past data. Another method considered as reference is based in climatological concepts [72]-[76], and use the average of the meteorological statistics accumulated during several years, to a specific location during a defined time interval. This method combines the persistence and the mean p as shown in equation, where the weight ak is a function of the correlation between the last measured value pt and the previous values. ( ) ˆ1 kt k t kt p ap a p += +− (2.63) Typically, this method obtains better performances than the persistence method in cases of forecasts from 12 to 18 hours [73],[74]. The drawback of this method is the need to estimate ak which has to be done under some assumptions. 2.5 Wind power forecast models In the prediction of very short-term and short-term, fundamentally there are two paths; one uses physical models and the other uses statistical models. However, there are systems using the combination of both methods, since, in reality, both are necessary to the forecasting success [72],[74]. 2.5.1 Physical models The physical models try, as much as possible, to only use physical considerations to achieve the best estimates to a specific place and, in certain cases, to use statistical models as Model Chapter 2 - Scheduling – General Overview 50 Forecasted precipitation Measured precipitation Very short-term (ANN) Short-term (ANN) Other input data Combiner Forecast after 1 h Forecast after 2 h Forecast after n h In [90] it was applied an ANN of MLP with backpropagation learning algorithm using as input the date (month and day), the inflows and the level of storage to provide the daily production for the next month. In [82] it is presented a study for daily river inflow forecasting, without considering power production, where a study has been done for comparing two forecast methods, namely Self-Organizing Neural Networks and Auto Regressive Moving Average (ARMA) based in time series. In another different approach, [91] studied the use of ANN of MLP to rainfall prediction for 6 h ahead. It is intended to know the inflows and avoid flash floods, mainly in small watersheds. It even was stated that in very small watersheds, sensible to flash floods, the predictions by numerical method may not be the most advisable. Many other topics based in ANN were studied, as Echo State Network, Self-Organizing Nonlinear Auto-Regressive model with eXogenous input (SONARX), ANN Radial Based Functions (RBF) and Artificial Neuro Fuzzy Inference System (ANFIS)[86]. In [92], the author also used ANN with several learning algorithms, as: Backpropagation (ANN_BP), Conjugate-Gradient (ANN_CG), Cascade Correlation (ANN_CC) and Levenberg-Marquardt (ANN_LM) to short-term forecast diary stream flow. In [93], are presented inflow forecasts based on gradient descendent method, Resilient Backpropagation, Scaled Conjugate Gradient and Levenberg-Marquardt which results were compared with the ARIMA method. Similar studies were presented in [94] where the daily forecast of inflows using ANN type MLP and RBF are presented. In this case, it was used as network inputs, past values of rainfall and runoff to forecast the inflow with MLP. As conclusion, in [92] and [95] it is sustained that, in average, 90% of hydrological studies applications use MLP with BP. 2.6.3 Ensemble of models Depending on the time horizon and available information and resources, there are several forecasting methods which can be used. However, they all have limitations, such as the need for large amount of historical data or the low ability to handle with non-linearity. With this, selecting combined models with their different characteristics and applications in order to contribute for the improvement of forecasting processes is a must. With the exposed, it is verifiable that there are several works proposing inflow forecasts with different time horizons, from hours up to months, but without proposing any kind of power Figure 2.14 – Reservoir inflow forecasting for a mini-hydro power plant based on rainfall forecasting Chapter 2 - Scheduling – General Overview 51 production strategies. These issues were recently outlined in [96] where it is presented an original short-term forecasting model for hourly average electric power production of small-hydro power plants. The proposed method takes into account operation strategies of the small-hydro power plants and starts with the estimation for the daily average power production, followed by hourly average power production. 2.7 Solar photovoltaic forecast models The electric power production in photovoltaic devices depends mainly on two variables, irradiance and cells temperature. The irradiance has direct influence on the cells current, keeping the voltage relatively constant during a large interval of irradiance values. As to the cell temperature, its increase causes a decreasing in the voltage, keeping the current almost constant. There is also another meteorological factor that indirectly has influence on the power production; the wind speed cools the cells’ surface consequently decreasing the cells’ temperature. Externally, there are other factors contributing to the power production, as the geographic position, the panel’s angle, among others. However, it can be considered that the solar power production depends mainly on global irradiance, which is composed by three components: direct irradiance, due to the direct incidence of solar radiation, diffuse irradiance due to the radiation that comes from atmosphere (although not being directly from the sun) and albedo resulting by the direct irradiance reflected from soil around. 2.7.1 Global irradiance forecast models Images from geostationary satellites have been used to determine and forecast the solar irradiance conditions in certain locations. This method is based on clouds shapes’ structure during time intervals, proceeding later to the extrapolation of its motion. The result is a prediction of clouds’ position and hence a prediction of solar irradiance to the studied site [97]. Solar predictions are done to very short-term based in satellite images and short-term with numerical weather predictions. In figure 2.15 a simplified process model is shown. Other approaches, based on ANN are proposed in [97] where ANN_BP, ANN_LM, ANN_RBF, recurrent networks and ANFIS to hourly solar irradiance forecast are used. In [98] are used ANN_BP, ANN conjugate gradient algorithm, ANN quasi-Newton algorithm and ANN_LM to the daily irradiance forecast. Chapter 2 - Scheduling – General Overview 52 Satellite images Forecast of cloud-index images Cloud-index images Forecast of global irradiance Forecast: • Motion vector fields • Smoothing Heliosat method Several other authors as [99]–[103] used time series models as ARMA and ARIMA with satisfactory results for short time intervals. In [104] it is considered that for very short-term, satellite images and ARMA and ARIMA can be used, in spite of the modelling of the daily irradiance, which cannot incorporate the clouds influence. 2.7.2 Irradiance splitting forecast models As mentioned before, all presented models predict the global irradiance in a horizontal plan, which is necessary to split into their components, namely, direct, diffuse and albedo. This step can be surpassed if the predictions are directly made in their components. In [105] it was proposed the forecast of daily and monthly diffuse components of global irradiance. Correlations between diffuse portion (quotient between diffuse and global irradiance) and clarity index (quotient between extra-terrestrial and global irradiance) were developed in a daily scale. In [106] several experiments were done with ANN, with several combinations of inputs to predict the hourly diffuse irradiance. 2.7.3 Power production forecast models To avoid the need of more calculations after the solar irradiance forecast, several authors presented models to predict directly the power production as those shown in figure 2.16. Starting with the numerical weather forecasts, with the predictions locally refined by statistical models associated with local meteorological stations. Next it is interpolated to solar power plants localization. After, it is necessary to simulate the photovoltaic system, considering its orientation and angle to convert the horizontal incident irradiance to an angular plan and split it in the irradiance components. Finally, the complete electrical system is modelled. Figure 2.15 – Example of cloud position prediction Chapter 2 - Scheduling – General Overview 53 In [107] the predictions are done in a perspective of energy production and the errors between individual power plant forecasts and aggregated forecasts are compared. It is concluded that the RMSE to one day ahead forecast is almost three times higher than to a national aggregation. 2.8 Load forecast models Load forecasting is another crucial component for the scheduling and operation of power systems. In operational framework such as power system scheduling, there is the need of knowing the load, which has to be fed by the units to be committed, as well as to estimate load flows for decreasing the risk of overloading, which leads to improvement of net load reliability and decreases the probability of occurrences of lines or units outages and blackouts. In planning framework it helps to make important decisions, namely decisions on purchasing electric power, load switching or infrastructure development. Said this, its accurate forecast is extremely important for energy suppliers, ISO, financial institutions and other participant’s on electric energy generation, transmission, distribution and markets. Comparing with some renewable power production, namely wind and solar, generally the load profile presents much less variability. It follows a well-defined pattern during the days with clear variations along the day between peak and off-peak periods. Over the years, generally, the great differences are in the load amplitude and some temporal shift linked with seasons and Winter/Summer hour. As referred in point 2.4, there is not a unique definition for classifying the temporal horizon of forecasting. In [108] it is stated that load forecasts can be split into three time frameworks, Midrange weather forecast Forecast of irradiance on horizontal plane Hourly site specific forecast of horizontal irradiance PV power forecast Location of PV system Post processing PV system orientation Model for irradiance on horizontal plane PV system characteristics PV simulation model Figure 2.16 – Power production forecasted model Chapter 2 - Scheduling – General Overview 54 namely short term forecast (form 1 hour to one week), medium-term forecast from a week up to a year and long-term forecast (longer than a year). In all time frames the majority of forecasting methods are based on statistical techniques as regression, or artificial intelligence algorithms such as ANN, fuzzy logic and expert systems [108]–[110]. In the particular case of medium and long-term forecasting, the so-called enduse and econometric approach or their combination is broadly used. These models use description of appliance used by costumers, size of houses, technology changes, equipment age, age costumer behavior and population dynamics. Also, economic factors as employment levels, per capita incomes and electricity prices are used. Although the variety of methods include several applicable to medium and long-term load forecasting, in this work only the short-term load forecasting methods will be addressed. In the case of short-term load forecasting there are three main aspects, such as time factors, weather data, and costumers classes [108]. The time factors include the week of the year, the day of the week, and the hour of the day. The week of the year gives information about the period of the year under study, as an example: it is harder to forecast during the holidays than non-holidays due to consumption profiles not following the frequent pattern. The day of the week describes the differences between working days and weekends as well as the different behaviors in working days, namely the Mondays and Fridays, which present a structurally different load from the rest of the days. This behavior has higher impact during summer. Concerning the weather, the temperature is the most important data, followed by the humidity and wind speed. The binomial temperature/humidity has more importance during summer and temperature/wind speed during the winter. The costumer classes is more significant at a local level since it gives significance to the type of costumers, such as residential, commercial, and industrial. Following [108], there is a large variety of statistical intelligence techniques which have been developed for short-term load forecasting such as similar-day approach, regression methods, time series, ANN, expert systems and fuzzy logic. The similar-day approach is based on searching historical data for days with similar characteristics as weather, day of the week and date. With this, the load of similar day is considered as a forecast. Regression methods based on time series are widely used as statistical techniques. Different authors introduce different explanatory variables such as, weather, type of the day, costumer classes, among others, in order to find their relation with load consumption. Methods such as ARMA, ARIMA, ARMAX and ARIMAX are often used [108]. In [110] semi-parametric approach to model nonlinear relationship combining temperature with time and type-of-day using ARIMA is addressed. Chapter 2 - Scheduling – General Overview 55 In [109] it is proposed the use a feedforward backpropagation ANN, which combines load profiles and temperature in order to find the non-linear relation between them and predict the load 24 hours ahead. There is other ANN which can be used, such as Radial Basis Functions (RBF), Hopfield, Boltzmann machine, among others. However, the feedforward backpropagation is still the most used. The expert systems are molten with fuzzy logic, since both are based in heuristic techniques, where the forecasting is based in rules and proceeding is defined by human experts. In these cases there is not any associated mathematical model, being the relation between the explanatory variables, weather and time, among others, being the load defined by linguistic rules. 2.9 Uncertainty estimation Thus far, in all the forecast techniques, whether they are renewable production or load, only a single forecasted value is provided. The main disadvantage is that no information about possible deviations from the predicted value is available. For the decision makers the benefits are quite limited, mainly in the applications based on risk assessment or stochastic optimization, as it could be seen in point 2.3. For this, it is much desirable for the decisionmakers having an idea about the uncertainty associated with each forecasting period. When a forecast is done, there are always associated errors which can result from several factors, such as incorrect or incomplete models, wrong parameters, extreme events, incorrect starting conditions, variations of sources, dynamics over the forecast period, amongst others. The models developed for load forecast generally get good results with lower deviations from the measured values because load profiles follow a characteristic pattern. Comparing with load forecast, the prediction of RES, due to its variability, presents much bigger challenges. In the RES chapter, due to its high installed power capacity all over the World and high variability and uncertainty, wind power forecast gathers the majority of the attention of researchers and the major number of published works. Though addressing other types of RES (namely solar and hydro), wind generation is presented as the main source of generation uncertainty in power systems scheduling, grid operation and market environment. In wind power forecast, there are three major factors which have influence on the uncertainty, namely, the NWP, the conversion of wind to power (due to the nonlinearity of the power curve) and terrain complexity. On the other hand, the NWP and the clouds dynamic are the main source of uncertainty in the case of solar photovoltaic since conversion is well defined. In the case of hydro power forecasting systems, the uncertainty generally propagates from the Chapter 2 - Scheduling – General Overview 56 NWP model through the rainfall-run-off model. The rainfall-run-off models are limited by their representation of flow dynamics, whose main problem is not the representation of the dynamic but knowing the local parameters [111]. In this analysis one must keep in mind that there is a difference between how to capture the uncertainty of forecasts and the way to represent the uncertainty of those predictions (by probabilistic models or scenarios). A very comprehensive analysis can be found in [1] and a very extended and complete in [71]. The following examples are focused in wind power production but can easily be extrapolated to other power sources. In figure 2.17 the based on NWP point forecast approach is depicted. In this case only a single point forecast is done to each look-ahead time, (generally wind speed and direction, atmospheric pressure, precipitation, temperature, among others). As represented in figure 2.17, with this data two paths are available, or it applies the data directly to the probabilistic model or it previously converts the NWP to wind power spot forecast (WPF) using a wind-topower (W2P) model, calculates the errors and only then applies the probabilistic model. As a single spot forecast is not enough to characterize the uncertainty, on this approach it is mandatory having a historical data set of power production, or errors and corresponding explanatory variables, which must be updated by the time each new value is known. Since wind power generation is a nonstationary process, a time adaptive and recursive estimation method can be applied [11],[71],[112] and [113]. These upgrades are very important when in presence of changes in the production profiles due to the installation of new generators, long outages owed to big maintenances or even something as natural as changes in the vegetation on the surroundings of the wind facility. There are several probabilistic models that can be used in order to represent the uncertainty of forecasts, which will be discussed in point 2.10. Figure 2.17 – Approach based on NWP point forecast Other approach for determining the probabilistic forecast is to produce ensembles. This method is different from the previous: instead of a spot forecast an ensemble of forecasts is provided. There are, fundamentally, two methods: one consists on different runs of a NWP model with different initial conditions or different numerical representation of the atmosphere, the other consists on using a different NWP modeling, or different forecasts made at different times [115],[116]. Basically, there are three different methods to attain the probabilistic forecast from NWP ensembles, the “filtering approach”, the “direct approach” Chapter 2 - Scheduling – General Overview 57 and the “dimension reduction”. In figure 2.18 is shown the filtering approach. After the knowledge of the NWP ensembles, it is done the conversion to wind power ensembles. If the ensembles are resultant from the perturbation of a single NWP model, only one conversion model is generally considered, since the ensembles may be considered, in general, indistinguishable. Otherwise, in the case of different NWP modeling or different forecasts models made at different times, it is used one conversion model for each ensemble member. Figure 2.18 – Filtering approach The output values of conversion models must be calibrated by a post-processing method in order to convert the uncelebrated power ensembles into probabilistic forecasts [1]. In figure 2.19 and 2.20 the direct and dimension reduction approaches are respectively shown. As in the case of figure 2.17, in the direct approach the probabilistic model can be fed directly with NWP outputs. However, in this case, it is used an ensemble instead of a single spot forecast. The main shortcoming is the increase of model complexity (when the number of input variables is big) without remarkable increment of results, even with the growth of sample number [1]. Figure 2.19 – Direct approach To reduce the complexity and transform it in a more tractable problem, a dimension reduction approach can be done before feeding the probabilistic model, as depicted in figure 2.20. In this case, the number of ensemble members can be reduced by aggregation or converting the ensembles into two values (mean and variance, for instance). Figure 2.20 – Dimension reduction approach Chapter 2 - Scheduling – General Overview 58 2.10 Uncertainty models As explained above, despite the advances in power forecast, there always are associated errors which depend of the resources that are forecasted, the prediction models, forecasting horizons and extreme conditions. The uncertainty created by these errors has a great impact on power systems scheduling since the forecasted values at the beginning of the scheduling process can be quite different from those in the operation periods. In a traditional and conservative point of view, generally, the uncertainties are compensated using conservative decisions, like overdesigning the equipment or overestimating the operational parameters basing them on worstcase of the uncertainty parameters. This approach, though being secure, may lead to significant results’ deterioration from the optimization problem. To overcome this situation, more complete uncertainty models must be provided, based on probabilistic forecasts, scenarios or risk indices. The probabilistic forecasts consist of estimating the future uncertainty of power resource expressed as a probability measure. The power production uncertainty can be described using random variables, which may be expressed by many forms, such as [1],[71]: • Moments of distributions (mean, variance, skewness, kurtosis); • A set of quantiles and interval forecasts; • Probability mass function (pmf); • Probability density functions (pdf) or cumulative distribution functions (cdf) (parametric and non-parametric). The use of any of the previous uncertainty representations is a case dependent, which should be determined by the nature of the decision-making problem or by the end user’s request. 2.10.1 Moments of distributions The moments of distributions can be used to represent parametrically the uncertainty in decision-making problems [71]. The moments are an acceptable way to the determination of the expected errors and are widely used due to its simplicity and easy calculation from the data sets. The data is shortened in classes of forecasted power and then the standard deviations are calculated [114]. The drawback of this method is the discontinuity between classes and the parameterization (number of bins and width). On the other hand, the standard deviation by itself does not give any information about the probability of a forecast error falling out within a specific interval. Chapter 2 - Scheduling – General Overview 59 2.10.2 Quantiles and intervals forecast The probabilistic forecast based on quantiles is a non-parametric approach which avoids the errors introduced by the wrong choice of a parametric distribution [115]. To introduce the concept of quantile forecast let’s consider ft+k as the probability density function of the random variable Pt+k as well as Ft+k as the cumulative distribution function. The forecasted quantile to look-ahead time t+k done at time t, ˆtkt q α + with nominal proportion [ ] 0,1 α ∈ can be defined as ( ) 1 ˆ ˆ tk tk tt qF α α − ++ = [71],[112]. When a single quantile forecast ˆtkt q α +with nominal coverage α is defined, it only gives the information that the random variable Pt+k has a probability α of being less than ˆ t kt q α + . As such, it is not given any information concerning any confidence interval. However, with more than one quantile it is possible creating confidence intervals (or prediction intervals) ( ) ˆ t kt I α + bounded by the forecasted quantiles, as shown in equation (2.64), [5],[112],[116] and [117]. ( ) ( ) ( ) ˆˆˆ , lu t kt t kt t kt I qq αα α + ++  = (2.64) These intervals define a range of possible values within which it is projected that observed values pt+k falls with a certain probability. This probability is defined by the nominal coverage (1α ) with α defined within [0,1]. These prediction intervals, centered on the median, are defined by its lower and upper forecasted quantiles, with lower and upper nominal proportions αl and αu respectively, and calculated by (2.65) and (2.66) [71]. 1 ul αα α −=− (2.65) 1 2 l α α − = (2.66) Generalizing this process and calculating n interval forecasts with various nominal coverage rates, allows the definition of predictive distributions. Thus, a probabilistic forecast made at time t for leading time t+k is given by the set of corresponding 2n predictive quantiles as shown in equation (2.67), or in a compact form in (2.68) [8]. ( ) ( ) ( ) ( ) ( ) ( ) { } 11 1 /2 /2 1 /2 1 /2 1 /2 ˆˆ ˆ ˆˆ ˆ ˆ , ,...., , ,...., , nn n n n t kt t kt t kt t kt t kt t kt t kt f q q qq q q αα α α α α −− − −− + + + ++ + + ≡ (2.67) ( ) ( ) ( ) { } 2 12 ˆˆˆ ˆ , ,...., n t kt t kt t kt t kt f qq q α αα + ++ + ≡ (2.68) In order to show the intervals forecast evolving along the forecasting intervals, some quantile regression techniques can be addressed. There are fundamentally three approaches, namely Chapter 2 - Scheduling – General Overview 66 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Power [p.u.] scenarios Forecasted Measured Figure 2.25 – Wind power forecast scenarios (23 scenarios from 1000) With the creation of these scenarios it is possible to represent the generation ramps’ limits between two time steps. This analysis is clearly explained in [112] and [113]. 2.10.5 Risk Indices or Skill Forecasting In the case of power forecasts based on meteorological phenomena, even with adequate and precise forecasting methods, when the weather stability is low, a new source of errors is introduced due to the powerlessness of the forecasting methods to deal with extreme weather conditions. This is the case of small hydro power production when facing wind power under unstable atmospheric conditions or unexpected floods due to snow melting. In the case of solar production, even in sunny days with sparse clouds, the cover of the solar panels during a few minutes can rapidly change the production. In meteorology environment and following the glossary of meteorology of American Meteorological Society, the skill forecasting (SF) or skill score, shown in equation (2.84) , is a measure which relates the forecast accuracy of a certain forecasting model in a comparison with a reference model. The result is a single value which gives an idea about the forecast accuracy of a certain forecast model comparing to some reference model [127]. Ref 1 Forecast SCORE SF SCORE = − (2.84) A perfect forecast results on a unitary SF, while performances which are equal to the reference model result on a null value of SF. If underperforming, the SF value is negative. There are, generically, two risk indices, the Meteo-Risk-Index (MRI) [71],[128] and Normalised Prediction Risk Index (NPRI) [71],[129]. From the risk indices it is possible to have a previous hint for the magnitude of the forecasting errors and to do a critical analysis concerning the level of accuracy of the expected forecasts. The MRI reflects the spread of the available NWP ensemble at a given time, which can be achieved by perturbing the initial conditions of the NWP model or by different NWP models Chapter 2 - Scheduling – General Overview 67 with equal initial conditions. The same spread can be created with lagged forecasts (same look-ahead time k but done at different t of t+k|t) applying different initial conditions to an unperturbed model. Following [129], the NPRI results from the NWP ensemble forecasts converted to wind power. The spread is computed with a weighted standard deviation of the ensemble members and can be analysed as the ability of each ensemble member to provide information on predictability. 2.11 Summary and main conclusions This overview gives a general idea about the met challenges in the scheduling with uncertainty. It is visible the interest that the scientific community still has concerning these issues. To fulfill the entire process there are several areas that must be investigated, each one conducting research areas addressed in this overview. The uncertainty associated to the power prediction gave a new interest to the scheduling, namely UC, ED and reliability assessment, introducing new models for the promotion of efficient analysis of the uncertainty. Follow-on from these problem there are some concerns about the need to increase the computational velocity. The reserves assessment created new challenges that have been investigated with good margin for future research. Concerning the load and renewable power forecasts, the models which produced point forecasts have reached good performances. However, researchers are now focused on uncertainty modelling. Consequently a wide number of researches concerning probabilistic forecasts, scenarios generation and risk indices are available. However, we were able to conclude that there is still space for research on complete scheduling tools in order to improve the robustness, velocity and computational burden reduction. Chapter 2 - Scheduling – General Overview 68 CHAPTER 3 Equation Chapter (Next) Section 1 Power forecast with uncertainty Contents In this chapter are presented some methodologies used in the forecasting of renewable power production, the load and the net load. Comparisons between probabilistic forecasting methods and evaluation criteria are also presented. Chapter 3 - Power forecast with uncertainty 70 3. Power forecast with uncertainty 3.1 Introduction In order to decrease the thermal generation fuel consumption and, consequently, decrease the power production costs in São Miguel Island, a large amount of renewable energy sources have been integrated. Nowadays, the generation capacity in São Miguel Island is ensured by 1 thermal power plant with 8 units (divided in 2 groups with different rated power), 2 geothermal power plants with 5 units, 7 small hydro power plants with 1 unit each and 1 wind power plant with 9 units. In table 3.1 the rated power to each power source is shown. Table 3.1 – Rated power of each power source Source (# units) total power Fuel (4) 28 MW (4) 64 MW Wind (9) 9,4 MW Small hydro (7) 5 MW Geothermal (5) 29,6 MW Other (2) < 1 MW To characterize the production and consumption in São Miguel Island, during 2012 were done hourly average measurements, and load variation between a minimum of 17,3 MW and a maximum of 70,2 MW was registered. The maximum renewable production reached was 43,8 MW. All this renewable production helps to decrease the thermal production during peak load periods but during off-peak periods the system is already saturated with RES power production. This situation leads to the necessity, most of the times, of limiting the wind power output in order to maintain the thermal generators operating above their technical minima. As an example, in figure 3.1 the load and production profile for a week in São Miguel Island is shown. 0 10 20 30 40 50 60 70 Power [MW] Geothermal Hydro Wind Thermal Load demand Figure 3.1 – Load/Production profile (one week) Chapter 3 - Power forecast with uncertainty 71 As depicted in figure 3.1 the load follows the traditional consumption profile with clear differences between peak and off-peak periods and between working days and weekends. In case of production, during almost all the time the geothermal power production is maintained constant and as the base of diagram, with only a slight decrease during a few hours on November 7th morning. As this happened in the peak load period, when more generation is required, it is expectable that this value is linked with some geothermal unit outage. Also, the hydro generation is generally constant. These hydro power plants have very low storage capacity, so the production is only dependent on the inflows and Azores Islands are characterized by frequent rains. On the other hand, wind generation (as would be expected) presents considerable variability along the week. This cannot be explained only with wind variability. Despite wind generation reach remarkable production during the peak load periods, during the off-peak periods it is reduced almost to zero. This is a clear sign of wind curtailment during off-peak periods (and this is not a particular case of São Miguel Island, in [130] several case studies concerning wind power curtailment are presented). The areas highlighted with red squares seem to be periods were there was wind curtailment, mainly because the reduction of wind power generation was exactly during off-peak periods. Notice that throughout the periods before and after the off-peak period, the wind generation was high. Seeing that in those off-peak periods the wind power generation remains without curtailment, the thermal production should strongly be reduced or even turned-off (which is not acceptable following the conservative approach done by the grid operator and the units’ operational limits). For operational security reasons it is mandatory that, at least, two thermal units must be on-line to avoid the complete loss of thermal production due to outages. It means that frequently, especially during off-peak periods, wind power production curtailments were necessary and, in some cases, to set the thermal units to work below their minimum technical limits. In figure 3.2, the hourly average thermal production, as well as the sum of minimum technical limits of 2 units, since 0h00 of December 1st up to 23h00 of December 31st is shown. 0 5 10 15 20 25 30 35 40 45 Power [MW] Figure 3.2 – Thermal production and minimum technical limits Chapter 3 - Power forecast with uncertainty 72 In these cases the thermal units were forced to work with poor efficiency and high fuel consumption while renewable sources are wasted. In 2012 those occurrences happened during 745 hours (17% of the year). As it is shown, during all over off-peak periods the thermal units worked below their minimum limits. To face this problem, it is clear the necessity of an efficient method to forecast load and renewable production in order to know the thermal production necessities. These forecasts, further than the spot values, must incorporate the uncertainty. This thermal production necessity is characterized as net load [31],[51] and [131]–[136], and it is calculated by the difference between forecasted load and the sum of forecasted renewable production. This way, the complete forecasting methodology involves the probabilistic forecasting of load and renewable production, in order to obtain the probabilistic forecasted values for the thermal production. It is intended to estimate the future conditional pdf of the random variable tk P + for each look-ahead time step t+k, given a learning set with N samples summarizing all historical information available up to the instant t. To develop this process, the data sets used in this work contain hourly average values since 0:00 of January 1st up to 23h00 of December 31st of 2012. The training/parameterization data set contain 7656 hourly average values since 0:00 of January 1st up to 23h00 of November 14th (87,2% of the total) and the test/validation set is composed by 1128 hourly averages values (12,8% of the total) since 0:h00 of November 15th up to 23h00 of December 31st. The NWP forecasts have 1 hour of temporal resolution and the forecast are available at 00h00 for 0h00 up to t+24. 3.2 Choice of forecast models As discussed in 2.10, in the literature there are several probabilistic estimation models. The choice of the model used in this work is based on analysis and comparison between several approaches available in the literature. In the case of point forecasts there are several standard error measures and evaluation criteria such as, Mean Error (resulting from systematic error), the Mean Square Error and Mean Absolute Error (resulting from systematic and random errors), Percentage Error, Mean Average Percentage Error (MAPE), and Standard Deviation Errors (SDE) (resulting from random errors). If calculated in function of installed capacity, it results on normalized values. Other error measures can be described as frequency distribution of the errors or coefficients of determination (R2) [72],[137]. However, evaluating probabilistic forecasts is harder than the evaluation of point forecasts. To evaluate probabilistic forecasts, there are four indicators that are widely used, namely, the reliability, sharpness, resolution and skill score [1],[8],[9],[11],[13] and [138]. Chapter 3 - Power forecast with uncertainty 73 The reliability is a measure of the agreement between nominal proportions (forecasted probabilities) and the computed from evaluation samples. The reliability diagrams give the empirical coverage versus the nominal coverages (proportions) for various nominal coverage rates. The nearer the diagrams are to a diagonal, the better the results. In alternative the diagrams can be drawn in function of the deviation from the “perfect reliability” in order to calculate the bias. In this case the deviation of empirical coverage in relation to the nominal one should be null. This is similar to the use of Probability Integral Transform (PIT) histograms, but in the reliability there is the added value to give the deviation magnitude to the “perfect reliability. The reliability is calculated with equations (3.1) up to (3.4). ˆ 1 0 tk t kt k if p q otherwise α α ξ ++ ≤  =   (3.1) { } ,1 , , 1 #1 N k ik ik i n αα α ξξ = = = = ∑ (3.2) { } ,0 , ,1 #0 k ik k n Nn αα α ξ = = = − (3.3) ,1 ,1 ,0 ˆk k kk n nn α α αα α =+ (3.4) In (3.1) pt+k is the realization of the variable and ˆ t kt q α + the forecasted quantile with the nominal coverage proportion α i. The result of this equation takes the value “1” if the realization of the variable pt+k hits the forecasted quantile with the nominal coverage α i and is “0” if it misses it. The number of “hits” and “misses” are counted by (3.2) and (3.3), ˆ k α α is calculated by equation (3.4) and gives the percentage of the hits . With the difference between empirical and nominal proportions it is possible to have an idea about the bias of the probabilistic forecasting methods and to measure the quality of the forecasts. Other indicator is the sharpness, which represents the capacity of the model to forecast extreme values. This criterion evaluates the prediction independently of the observations and gives an indication of the level of predictions usefulness. It measures the probability of forecasting values with probability falling near 0 or 1 instead on 0,5 (which is the dispersion around the 0,5). Following the equation (3.5), the quantiles are gathered by pairs in order to obtain intervals with different nominal coverage rates from narrow intervals up to wider intervals. The sharpness of the predictive intervals is measured by the average interval size, equation (3.6), where t kt α δ + is the size of the interval forecast with nominal coverage rate 1α estimated at time t for lead time t+k. Chapter 3 - Power forecast with uncertainty 74 22 1 , ˆˆ tk t kt t kt qq αα α δ − ++ = − (3.5) , 1 1N k tk t N αα δδ = = ∑ (3.6) With larger intervals it is possible to predict values with large dispersion around the median or even outliers but this does not mean that it is a good prediction. Larger intervals mean that the approach is very conservative and it is not designed to “take risks”. The reliability and sharpness have contradictory results: generally good values of reliability conduce to bad values of sharpness and vice-versa. For this, a trade-off between reliability and sharpness has to be accepted. The skill score, by equations (3.7) and (3.8), gives information about the model performance in a single measure for a set of m quantiles, ( ) ( ) ( ) 1 ˆˆ , ii m c tk i tk t kt t kt i Sf p p q αα ξα ++ ++ = =−− ∑ (3.7) ( ) ( ) 1 1 ˆˆ ,, N ck tk c tk t kt t kt t Sfp Sfp N ++ ++ = =∑ (3.8) where pt+k is the realized forecasted, α i is the quantile proportion, ˆ i tk q α + is the forecasted quantile and i α ξ is the indicator resulting from (3.1). The higher the value, the better the skill score with the maximum value of 0 for perfect probabilistic forecasts [6],[11]. Finally, the fourth indicator used for evaluating probabilistic forecasts is the resolution [6]. This indicator, used in case of intervals forecasting, can be calculated by the standard deviation of the size of the intervals. Similarly to other indicators, the resolution is dependent on the case being studied. A case with less uncertainty (small sharpness) generates less variability in intervals’ size. In opposition to the sharpness, as higher is the resolution, the better is the model [6]. In conclusion, sharpness is related to the average size of prediction intervals, whereas resolution is measured with the variability of their size. These indicators were used by several authors to study the performances of many probabilistic forecasts and will serve as base for the choice of the probabilistic forecasts used in this work. In [1] a comparison between Spline Quantile Regression (SQR), Quantile Regression Forest (QRF) and KDE is done, with the Linear Quantile Regression (LQR) serving as comparison base. Regarding the reliability, it is concluded that the QRF and KDE present equivalent performances whereas sharpness achieved nearly equal values for all approaches. In [6] the KDE and QRF were compared with B-splines Quantile Regression and it is concluded that KDE presents encouraging results of reliability and sharpness towards the remaining techniques. In [9] and [11] it is presented a time adaptive conditional KDE and its Chapter 3 - Power forecast with uncertainty 75 performances are compared with other probabilistic forecasts approaches, namely LQR and SQR in terms of reliability, sharpness and skill score. As a result, KDE showed to have a trend to present a better reliability while quantile regression presents the tendency of better sharpness. The skill score was quite similar to both approaches. Also in [118] it is done a comparison between LQR, Local Gaussian model and KDE, concerning reliability and sharpness. It is concluded that the provided examples did not give any clear preference between methods. It was also stated that regarding simplicity and implementation, the results suggest the KDE estimator as a good choice. Facing these results, it can be concluded that there is not an approach that stands out with clear improvement over the others. Given these conclusions, in this work, the probabilistic forecasting method to define the pdf as well as the expected values for each time-ahead forecasts is based on KDE, whose mathematic formulation was presented in point 2.10.3 and from equation (2.73) up to equation (2.77). 3.3 Choice of the Kernel function When using a KDE approach, the choice of kernel function is the first step. There are several possible functions that can be used for each variable. In the statistics literature are presented several kernel functions namely, Gaussian, Epanechnikov, Biweight, Triweight, Tricube, Cosine, Logistic, Silverman’s, among others. Generally, the most broadly used is the Gaussian kernel function, though in [122] it is stated that the Epanechnikov kernel has a slight improvement regarding the Gaussian kernel. In reference [6] it was used the Biweight kernel instead of the classical Gaussian kernel just to decrease the computational efforts. Following [14], and considering the wind power as a bounded variable, both Gaussian and truncated Gaussian kernels were compared, but without any practical benefits. As conclusion, in [1] and [6] it is stated that the kernel function when compared with the selection of bandwidth h has a minor impact on the estimation quality. On the other hand, and particularly to wind power forecast in [9] and [11], it is considered that the choice of the kernel function is a critical issue. The authors proposed different kernels to different forecasting variables, Beta kernel to the wind power (because it can be bounded between 0 and 1 p.u.), Gamma kernel to the wind speed (between 0 and + Inf.) and von Mises distribution to wind direction (circular between 0 and 2π). Subsequently to this analysis it was decided to choose as kernel the Gaussian function; it is widely used in the literature and it is easy to implement. In comparison, the Epanechnikov, Biweight, Triweight, Tricube and Cosine kernel functions demand the standardization of the values between -1 and 1. Moreover, there are no strong signs of the superiority of the remaining approaches. Chapter 3 - Power forecast with uncertainty 82 There are several techniques to optimize the values of h although many authors use a simple trial-and-error technique or some rules of thumb as the Silverman’s rule. More sophisticated approaches can be used, as Leave-One-Out Cross Validations (LOOCV) presented in equation (3.11), [110] and [122]. [ ] ( ) 2 , 1 1ˆ argmin N ii hi hi h Yf X N − =    = −     ∑ (3.11) Where (Yi,Xi) is the ith observation and [ ] , ˆhi f− is the estimated function, omitting the ith observation using bandwidth h. In this work, wind speed and direction smoothing parameters were calculated by LOOCV while wind power was calculated by trial-and-error. In table 3.2 the smoothing parameters obtained are shown. Table 3.2 – Bandwidth hj for wind power forecast Kernel’s Variable hj Wind speed 1,7 Direction 10 Wind power 0,8 As example, in figure 3.11 the conditional probability distribution of wind power forecast for forecasted wind speed of 12 ms-1 and direction of 247º is shown. 0.00 0.01 0.01 0.02 0.02 0.03 0.03 0.04 0.04 0.05 0.05 Probability Wind power [MW] Figure 3.11 – Conditional probability distribution of wind power forecast For each look-ahead time of the forecasting horizon, knowing the forecasted values of wind speed and direction, it is possible to compute a probability distribution for each time ahead. The point forecast value results from equation (2.77). In figure 3.12 the wind power forecasted and measured, as well as the 80% uncertainty interval is shown (following [13], a nominal coverage between 0,75 and 0,85 seems a good compromise). The forecasts were done at 0h00 up to 24 hours ahead and cover a time horizon of 168 hours. Chapter 3 - Power forecast with uncertainty 83 0 1 2 3 4 5 6 7 8 9 10 Wind power [MW] Uncertainty Forecasted Measured Figure 3.12 – Hourly measured and forecasted wind power generation Outside the off-peak periods, where there was no limitation in production, the measured values are reasonably covered by the uncertainty interval. A more accurate forecast is particularly hard to do, since there is only one power plant throughout the island consequently there is no smoothing effect resulting from some errors compensation. The forecasting performances indicators cannot be directly applied, since measures do not depend exclusively of the explanatory variables. 3.6 Hydro power forecast Despite the low rated power, compared with remain generation technologies implemented in São Miguel Island, the hydro power production was modelled and forecasted in order include its contribution to the generation mix. The total hydro production capacity is concentrated in 7 small hydro power plants (SHPP’s), with only one unit each. Except one, all present a rated power less than 1 MW. In table 3.3 the rated capacity of each facility is shown. Table 3.3 – Small hydro power capacity in São Miguel Island # units Power [kW] CHTN 1 1658 CHTB 1 94 CHFN 1 608 CHCN 1 400 CHFR 1 800 CHRP 1 800 CHSC 1 670 All the presented SHPP have low storage capacity and small watersheds, which means that the hydro production is strongly connected with rainfalls where the SHPP’s areas are located. In figure 3.13 the aggregated hourly average power production and hourly average precipitation in São Miguel Island during 2012 is shown. Chapter 3 - Power forecast with uncertainty 84 0 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Precipitation [mm] Power [MW] Hydro power Precipitation Figure 3.13 – Hourly average hydro power production and hourly forecasted precipitation Observation of figure 3.13 allows the verification that some variability on power production is inconsistent, which may be related with failures on data acquisition or outages. On the other hand, as discussed in point 2.6.1, the empirical regularities or periodicity are very often masked by noise. As such, it is comprehensible that it is not easy to establish a direct mathematical regression between rainfall and electric power production. These difficulties are related with several issues, such as: • Frequent phase-shift errors (temporal deviation between the forecasted and real moment of precipitation); • Delay between the instant when the precipitation occurs and the instant that the hydroresource arrives to the SHPP reservoirs (which depends of the watershed characteristics); • The electric energy production is mainly influenced by the operation strategy of the SHPP’s (the timing for generation is practically independent of the time of the precipitation). These reasons suggest the use of a moving average filter in the hydro power production and precipitation, leading to a better understanding of the relationship between both variables, as shown in figure 3.14. The data was filtered by a 24 values moving average centred in an hour h, which lead to a smoothing of high frequency noise and possible wrong measures. Although the aggregated production seems to be relatively constant throughout the year, there are some variations in hydro power production, with a fast increasing during periods with precipitation and decreases in production over dry periods. Chapter 3 - Power forecast with uncertainty 85 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 precipitation [mm] Power [MW] Hydro power Precipitation Figure 3.14 – Daily average hydro power precipitation and HPP In this work the hydro power forecasts are based on the methodology proposed by [96], named H4C, which is composed of three modules. The first is used to estimate the daily average power production of the SHPP, the second takes into account the operation strategy of the SHPP and the third is an assimilation module. The methodology used in this work, uses as inputs variables, past values of power production of SHPP and precipitation as along with forecasted values of precipitation (from an NWP tool). The first step, in order to create a relation between the forecasted values of precipitation and hydro power production, consists on the introduction of the concept of hydrological power Potential (HPP). The variable HPP represents the level of hydrological potentiality to produce electrical power and is defined by equation (3.12), ( ) 1 ˆˆ . hhh H B H AR − = + (3.12) where ˆ h H is the forecasted HPP to the hour h and depends on the previous forecasted HPP, Hh-1. The parameter A is related with the incremental response of the electric power generation to the precipitation (MW/mm/h), parameter B is dimensionless and related with the speed of decrease of such generation in dry days (its value must be positive and lower than 1). The variable ˆh R is the forecasted precipitation for hour h, in mm/h. In other words, parameter A is connected with the incremental power production due to the precipitation, while parameter B represents the decrease of the power production in dry days. As the area of the Island is relatively small, the precipitation is considered uniform all over the island so it was chosen a point near all SHPP and it is considered the average precipitation of the set of the points. Figure 3.15 shows the temporal evolution of the hydro power aggregated production of all SHPP’s, as well as its HPP. The parameters A and B were obtained by the least squares method where the fitness function was the absolute difference between the average power production and the corresponding HPP value. Chapter 3 - Power forecast with uncertainty 86 0 0.5 1 1.5 2 2.5 3 3.5 4 Power [MW] Hydro power HPP Figure 3.15 – Average hydro power production and HPP It can be seen that the HPP presents a good adherence to the average hydro power production. Nevertheless, some differences between the evolution of HPP and power production can occur due to nonlinearities associated with the generation limits. Thus, to model this nonlinearity, in [96], is proposed the function (3.13) to convert HPP in power production, max min min _ˆ 8 ˆ 1 hc s est h Hh h PP PP e  − −    − = +   +   (3.13) where _ ˆest h P is the forecasted power production for hour h, by the explanatory variable Hh, and Pmax and Pmin are the maximum and minimum power production, respectively. The parameters hc and hs are calculated by the equation (3.14), 2 om c s mo hh h hhh + =   = −  (3.14) where ho is the value when HPP reaches the minimum value of power production and hm is the value when it reaches the maximum. Following on from this process, the parameters used to calculate HPP and hydro power production are shown in table 3.4. These values were optimized by least squares method. Chapter 3 - Power forecast with uncertainty 87 Table 3.4 – Parameters used to hydro power forecast Parameter Value Pmin 2,18 Pmax 3,47 ho 1,855 hm 3,746 A 0,00855 B 0,99976 In figure 3.16 is depicted the measured hourly hydro power production during 2012 as well the calculated HPP and consequent hydro power point forecasts. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Hydro power, HPP [MW] Measured HPP Forecasted Figure 3.16 – Measures vs hydro power forecasted and HPP Using equation (3.13), in the presented case, the fitting of power forecast to the measured values is slightly better than with HPP, reaching low MAPE values. For a better perception about the results that this method is able to reach, a hydro power forecast was done for a period of 7 days (168 hours) divided in frames of 24 hours. To represent the uncertainty, the values of forecasted HPP were used as inputs of a KDE whose bandwidth values are depicted in table 3.5. Once again the bandwidth of HPP forecast was calculated by LOOCV method whereas for hydro power it was calculated by trial-and-error. Table 3.5 – Bandwidth for hydro power forecast Kernel’s Variable hj Hydro point forecast 0,08 Hydro power 0,1 The results are shown in the figure 3.17, namely the measured and hydro power forecasted with 80% uncertainty interval, as well as the measured precipitation. Chapter 3 - Power forecast with uncertainty 88 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Precipitation [mm] Hydro power [MW] Uncertainty Precipitation Measured Forecasted Figure 3.17 – Hourly measured and forecasted hydro power During this period, the forecasted power production presented small variations and the uncertainty interval remained almost constant. The measured values fall inside the uncertainty interval except on two occasions: On the first, despite the weak rainfall, the big production loss resulted from two power plants outages or measurement failures. Opposed, on the second, there were heavy rains, which substantially increased the power production. With the data from figure 3.17 it is possible to understand that there is a strong correlation between the rainfall and the hydro power. As in other methods, H4C can be affected by a BIAS error, due to possible errors in the forecasted precipitation which may imply significant deviations in the values of HPP. These deviations can be disseminated over the time, introducing significant forecast errors. The correction of HPP values allows reducing forecasting errors and can be done whenever new data is available. When the measured data of hydro power production is available for hour h (Ph), it can be used to adjust the HPP calculating Hnew_h by equation (3.15) which represents the new adjusted HPP for hour h[96]. ( ) max min _min 1 .ln . 1 8 s new h c h h H h PP PP  =− −−  −  (3.15) 3.7 Geothermal power forecast Geothermal power plants are considered as a renewable power source. Though, due to their specificities, they rarely are subject of study. Their implantation has to be done near areas with volcanic activity and as the resource is easy to control the variability is very low (when compared with other RES). In São Miguel Island there are two geothermal power plants namely: central geotérmica da Ribeira Grande (CGRG) and central geotérmica do Pico Vermelho (CGPV). The power capacity and number of units is depicted in table 3.6. Chapter 3 - Power forecast with uncertainty 89 Table 3.6 – Geothermal power capacity in São Miguel Island units Rated power/unit Total capacity CGRG 2 2,9 MW 16,6 MW 2 5,4 MW CGPV 1 13 MW 13 MW With 26,9 MW of power capacity it is the second source of energy in the island and, as depicted in figure 3.1, it ensures the base of load diagram. However, due to some over dimensioning, during the years, the production never reached the rated power, presenting an average production around 21 MW. In figure 3.18 the geothermal power production on Ribeira Grande and Pico Vermelho, during 2012 is shown. 0 2 4 6 8 10 12 14 16 18 20 22 24 Geothermal production [MW] Ribeira Grande Pico Vermelho Total Figure 3.18 – Hourly average geothermal power production all over 2012 Generally, the geothermal power production does not present great variability and it remains relatively constant around the set point power. Exceptionally, during most of the time of 2012 the CGRG remained unavailable and, when it was online, the apparent variability resulted from malfunctions of the power plant. It only restarted the production in August and a more constant production was only achieved by October. Since this date, the higher variability of total production can be explained by unit’s outages or errors in data acquisition. As the production is generally constant, the production forecast is based on a set point. In absence of information from the grid operator about the real set point, the spot forecast resulted from a moving average of the previous values. In order to model the uncertainty, the same proceedings as done in previous cases were implemented. The set point forecast is the unique explanatory variable and the bandwidths are shown in table 3.7. Table 3.7 – Bandwidth for geothermal power forecast Kernel’s Variable hj Point forecast 0,25 Geo power 0,5 In figure 3.19 the forecasted and measured geothermal power as well as the 80% uncertainty interval is shown. Chapter 3 - Power forecast with uncertainty 90 15 16 17 18 19 20 21 22 23 24 25 Power [MW] Uncertainty Forecasted Measured Figure 3.19 – Measured and forecasted geothermal with uncertainty It is visible that the production present a much lower variability compared with the remaining RES. 3.8 Load forecast Load forecast is one of the most important information for the scheduling process, being its knowledge fundamental to the UC, ED, security assessment and reserve capacity allocation, whether in deterministic or in stochastic formulations [4]. When the scheduling is deterministically formulated, the load forecasted values are generally provided as spot values (expected average demand) being the forecasted errors covered by a certain amount of predefined reserve [139]. In the case of stochastic formulation, due to the load (or production forecast uncertainty), the reserves can also be fixed or allocated dynamically in function of uncertainty sources [46],[59],[62],[63]. As in other forecasting process, the main concern is to choose the better explanatory variables which contain relevant information for better behaviour characterization of future load values. These explanatory variables are then used as inputs to a KDE in order to model the stochastic behavior of load forecast. Generally, the load has a typical behavior all over the 24h of a day, with off-peak periods during the night, beginning to grow early morning up to peak load period during the industry working period and decreasing again by the evening up to the off-peak. Moreover, the load behaviour presents differences between working days, weekends and holidays. In figure 3.20 the load profile along the 24 hours of all days of 2012 in São Miguel Island, is shown. Different values registered for each hour along the day and within each hour are clearly shown. On the other hand, in figure 3.21 the load behaviour for each hour all over the weekends of 2012 is depicted. Clearly the behaviour during weekends is different with regards to amplitude and shape. Chapter 3 - Power forecast with uncertainty 91 0 10 20 30 40 50 60 70 80 04812 16 20 24 Load demand [MW] time [h] Figure 3.20 – Behaviour of load for each hour all over 2012 (all days) 0 10 20 30 40 50 60 70 80 0 4 8 12 16 20 24 Load demand [MW] time [h] Figure 3.21 – Behaviour of load for each hour all over the weekends of 2012 As such, besides from the hour of the day, one of the explanatory variables should be the day of the week. In figures 3.22 a) to d), the load behaviour for different weeks of 2012 is depicted. It is shown that even for the same days of the week, the profile changes with different weeks along the year. 00 10 20 30 40 50 60 70 80 0246810 12 14 16 18 20 22 24 Load demand [MW] time [h] 1st week Wednesday Sunday Saturday 00 10 20 30 40 50 60 70 80 0 2 4 6 8 10 12 14 16 18 20 22 24 Load demand [MW] time [h] 13th week Wednesday Sunday Saturday a) b) 00 10 20 30 40 50 60 70 80 0 2 4 6 8 10 12 14 16 18 20 22 24 Load demand [MW] time [h] 25th week Wednesday Sunday Saturday 00 10 20 30 40 50 60 70 80 0246810 12 14 16 18 20 22 24 Load demand [MW] time [h] 37th week wednesday Sunday Saturday c) d) Figure 3.22 – Behaviour of load for different first weeks of 2012 This behaviour depends on many issues such as like weather conditions, summer or winter hours, holidays, school calendar, population growth due to tourism, among others. Following Chapter 3 - Power forecast with uncertainty 98 These parameters have to be calculated for every forecasted look-ahead time-step, since they are related with the expected value and variance of the distribution. The drawback is the time consuming, which can be reduced choosing a limited kernel as in [6]. In this work, in order to increase the speed, the model was programmed in a Graphic Processing Units (GPU) [142] by the Smartwatt team. 3.10 Summary and main conclusions The necessity of forecasting the load and production to optimize the thermal production in the scheduling process is a clear need. The classical point forecasts are not enough to characterize the uncertainties associated with the load, and furthermore the renewable power, which is strongly characterized by uncertainty. Analysis of challenges and the critical importance of the forecasts was accomplished and the awareness that quality of the forecasts is mandatory is one of the main conclusions. There are several probabilistic forecasting techniques, all with comparable performances, including the widely disseminated technique based on NadarayaWatson estimator, with conditional and marginal probability densities calculated with KDE. There are several kernel functions, in this work it was chosen the Gaussian function because it is one of the most used. In this work, wind generation reveals an extra challenge due to the wind curtailment, which is not an explanatory variable, forcing to take several approaches for dealing with it. The small area of the island and the reduced number of wind power plants and hydro power plants do not allow a real smoothing effect for decreasing the variability. It is also concluded that the choice of the explanatory variables is fundamental to the success of the forecasts and a careful choice of the smoothing factor can strongly influence the shape of the probability distributions. It is shown that convolution allows to aggregate probability distributions in order to calculate the net load forecasting. Finally, the probability distributions were approached to a Beta pdf in order to reduce the computational efforts and work with a very flexible and easy to implement parametric pdf. The complete study done in this chapter allowed the creation of a set of tools and proceedings focused in São Miguel Island case, which will support the remaining studies along the work. As conclusion, a good forecast of the net load could give trustable information to the grid operator in order to improve the grid operation. It is also notorious that the quality of the forecasts is essential to the success of the scheduling. Equation Chapter 4 Section 1 CHAPTER 4 Generation scheduling under uncertainty Contents In this chapter a full methodology for the generation scheduling under uncertainty, based on risk assessment is presented. The concept of equivalent generator is introduced and an original contribution of a metaheuristic in order to solve economic dispatch is addressed. Chapter 4 - Generation scheduling under uncertainty 100 4. Generation scheduling under uncertainty 4.1 Introduction Increasing introduction of electric energy production with RES, and mostly those with high variability, has created several challenges to the energy networks operators, especially in the scheduling. This problem is boosted in low power networks, particularly in islands without any connection to continental networks. Large variations on renewable production can introduce stability problems in the network, which can originate generation or load shed and, at limit, black-outs a strong possibility. With this into consideration and for security, scheduling is generally done by a conservative way, with low risk, although sometimes far away from an optimal operation. As such there is the necessity to introduce uncertainty of load/RES in scheduling for achieving a better management of the thermal unit’s commitment. 4.2 Generation scheduling under uncertainty When available, the RES production allows thermal production decrease, especially during the peak load periods. Optimizing the number and the power of the on-line thermal units lowers cost and emissions. On the other hand, an extreme reduction of the thermal committed capacity can lead to a situation where the spinning reserves are not sufficient to handle with great variations of load, renewable production or generation outages. Therefore, due to the uncertainty in load and renewable production forecast, it is sometimes hard to find a completely robust/economic scheduling solution. As analysed in chapter 2, the stochastic programming is an approach widely used to deal with the generation scheduling under uncertainty applying recourse problems, chance-constrained or robust optimization, with uncertainty described by scenarios [4],[15-16],[31-32] and [5767]. The scenario-based approach demands a great number of realizations in order to capture the temporal interdependence of the probabilistic behaviour of the uncertainty. One of the main problems of this approach is that it is time consuming to solve all scenarios, being necessary to appeal to some kind of scenario reduction. In all of the approaches presented in the scheduling overview, the UC and the ED are solved for each scheduling scenario, increasing the computational efforts with the crescent number of scenarios. To overcome this problem, this work develops a short-term scheduling approach to be used in insular power grids based on risk assessment, addressing the increase of variability and uncertainty created by RES. Chapter 4 - Generation scheduling under uncertainty 101 4.3 Description of the methodology The proposed generation scheduling is designed to minimize the sum of the estimated costs based on risk cost analysis. These costs result from the estimated normal operation cost plus the estimated cost of operating outside normal conditions. It is understood as “abnormal” conditions if there is the necessity of load shed due to the lack of available thermal production or RES curtailment caused by the lack of load. Following [143], decisions have to be made to accept a risk as long as it can be technically and financially justified. Contrary to widely used scenarios-based approach, in this work it is proposed the probabilistic estimation of costs based on estimation of risk, directly using the probability density function of the random variables. Stochastic programming allows the implicit determination of the reserves, by incorporating explicitly the stochastic nature of the uncertainty, with scenariobase. In the methodology proposed in this work it is also not considered a predefined value for the reserves, since reliability and operational risk minimization are expected to lead to solutions with enough levels of dynamic reserves. It is intended to evaluate the adequacy of each possible thermal GENeration mix SET (GENSET) in the UC for each hour under a probabilistic net load forecast. This way, the need of development of a large number of scenarios, by modelling explicitly the uncertainty of the forecast is avoided. For each hour h of the scheduling period (h = 1..H), the ability of each thermal generation mix set (GENSET) to meet the net load is verified. The risk of load shed or RES curtailment and thermal production below the technical minimums are used to define the objective function, as well as the probability of the thermal generators operating inside the appropriated ranges. This process is done in an independent way for each hour ahead (single period) being hereafter the probability of unexpected outages and start-up costs incorporated in the problem. An N-1 contingency of thermal units is taken into account and the start-up costs are integrated using a dynamic programming, based on the values obtained from the single period approach. The main objective is to choose the best configuration of thermal units in order to minimize the costs of unit commitment up to H hours ahead. The proposed methodology is divided into two stages; the first consists of a pre-processing which is done only once (offline) and the second, processed on-line. In figure 4.1a) a block diagram resumes the procedures of pre-processing stage. The proceeding starts with the determination of the fuel consumption characteristic of each thermal unit. Knowing the specific consumption in function of the power production it is possible to define a mathematical expression. Chapter 4 - Generation scheduling under uncertainty 102 Determination of fuel consumption curves for each thermal unit (with Specific Fuel Consumption) Definition of all thermal units combinations (GENSET) (Characterization of technical limits) Database with optimal production (each GENSET in function of net load) Resolution of economic dispatch (each GENSET in function of net load) Uncertainty aggregation (net load) Evaluation of each GENSET (Risk assessment and costs calculation) (Without start-up costs and contingencies) Contingency analysis (N -1 thermal unit) Dynamic programming (Startup costs) Scheduling Probabilistic forecasts Load demand Renewable production Single period Multi period a) b) Figure 4.1 – Proposed methodology Generally, these characteristics are given by the manufacturer and are evaluated by tests. With the set of thermal units it is created a dataset with all possible combinations of thermal mixes, defining a GENSET to each combination, with the respective maximum and minimum limits of generation. For each GENSET it is solved an ED for different values of net load allowing to define an equivalent generator. All results are stored in a database. Although exhaustive and very time consuming, this procedure is done only once, being updated when there is a change in the number or rated power of thermal units. The second stage procedure, shown in figure 4.1b), is always done whenever the scheduling performed. The load and RES production forecasts are received and aggregated as net load, which will be the procedure input. Knowing the probability density function (Beta) of the net load for each hour and all combinations of GENSET, it is done the risk assessment of each GENSET to be able to feed the net load, and the risk costs of the GENSET are computed. Next, is done a contingency analysis for (N-1). These procedures are recalculated to each hour of the scheduling period. At the end it is done a multi-period UC resulting in the solution with lower risk and lower cost with risk embedded. Chapter 4 - Generation scheduling under uncertainty 103 4.3.1 Thermal generation characterization In São Miguel Island the thermal power production is ensured by one power plant with eight generators divided in two sets, each one composed by four identical groups (prime mover plus generator) with the same rated power. The units are designated (1-4) for the set of machines with lower rated power while the designation (5-8) indicates the remaining units. The technical operation limits in steady-state of each set of groups are: 14 58 3848 kW 7200 kW 8410 kW 16500 kW GG GG P P − − ≤≤ ≤≤ With exception of start-up and shutdown periods (where the fuel is common diesel), the synchronous generators are powered by an internal combustion machines fed with heavy fuel oil. More information concerning the generation units is given in table AI.1 at annex I. In table 4.1 it is presented the specific consumption for each kind of machines for different values of power production. The production percentage in table 4.1 is calculated relatively to the maximum allowable power in steady-state of each generation unit. Table 4.1 – Specific consumption for thermal units % of rated power Rated power [kW] Specific consumption [g/kWh] Units 1 - 4 Units 5 - 8 Units 1 - 4 Units 5 - 8 50 3848 8410 222 218 75 5772 12615 213 207 100 7696 16820 212 205 With the values of specific consumption it is now possible to define a continuous function, generally modelled by a second order polynominal. In figure 4.2 is depicted the specific fuel consumption of each type of thermal units in g/kWh. 200 205 210 215 220 225 05000 10000 15000 20000 Specific fuel consumption [g/kWh] Power [kW] Polinomial (Units 1 - 4) Polinomial (Units 5 - 8) Figure 4.2 – Specific fuel consumption for each type of thermal units Chapter 4 - Generation scheduling under uncertainty 104 These curves represent the relationship between the specific fuel consumption and the power generation [144]. Examining this figure it is visible that specific consumption is remarkably higher when the units are operating at low power. Considering that the specific consumption continues to increase for even lower power values, this means that, beyond other technical issues, when the units operate below minimum limit, the efficiency is even lower. On the other hand the efficiency increases with the power growth, meaning that the maximum efficiency is near the maximum power. Given the specific fuel consumption (g/kWh) of each thermal unit and multiplying them by the power production (kW) it is possible to draw the fuel consumption curve (g/h). When the fuel price is known it is possible to create the cost function (€/h), characterizing the hourly cost of producing a certain amount of power. Generally, it is preferable to work with fuel consumption functions instead of cost functions because fuel costs are variable. Following [145], when the thermal units do not present valve-point effects, a second order trend line can be used to represent the fuel consumption values and consequently the generation costs. In figure 4.3 it is shown the fuel consumption for each type of thermal units. 0 500000 1000000 1500000 2000000 2500000 3000000 3500000 4000000 05000 10000 15000 20000 Fuel consumption [g/h] Power [kW] Polinomial (Unit 1 - 4) Polinomial (Unit 5 - 8) Figure 4.3 – Fuel consumption for each type of thermal units Applying a second order trend line to the fuel consumption points of figure 4.3 results in equations (4.1) and (4.2), where Pi is the power production of each thermal unit. These functions can easily be converted in currency multiplying them by the fuel cost. ( ) 32 14 3,89 10 160 172800 Units i i i FC P P P − − =× ++ (4.1) ( ) 32 58 1,70 10 150 445500 Unit i i i FC P P P − − =× ++ (4.2) 4.3.2 Thermal units combinations (GENSET) Once known the number of available thermal units in power system, together with minimum and maximum technical limits of each unit, it is possible to aggregate the units by defining all Chapter 4 - Generation scheduling under uncertainty 105 the available combinations. In this work, as there are 8 units, it is possible to define 255 combinations of thermal production mix and subsequent available power to feed the net load. However, as there are only two types of generators, in a restricted analysis of available power, the number of combinations can be reduced to 24 regardless of which thermal unit is on-line, this allows to create a much more tractable data set. For a matter of simplicity and to an easier understanding, the units (1-4) will be represented by GS (small power) and GB (big power). Table 4.2 summarizes the 24 possible unit’s combinations. Table 4.2 – Possible combinations of GENSET’s In figure 4.4 it is depicted the power production limits of each combination of thermal units. It is observable that, due to a finite set of GENSET’s, small variations in the production limits force several changes in the unit commitment. For instance, in the case of 2GS_3GB, the minimum limit is 32,9 MW and for decreasing this limit to 32,2 MW (a 700 kW decrease), it is necessary to adopt the configuration 4GS_2GB, which implies 3 manoeuvres. In this particular case two thermal type GS must by started and one type GB must be switch off. 3,8 7,7 8,4 11,5 12,3 15,4 16,1 16,8 20,0 20,7 23,8 24,5 25,2 28,4 29,1 32,2 32,9 33,6 36,8 37,5 40,6 41,3 45,249,0 7,2 14,4 16,5 21,6 23,7 28,8 30,9 33,0 38,1 40,2 45,3 47,4 49,5 54,6 56,7 61,8 63,9 66,0 71,1 73,2 78,3 80,4 87,6 94,8 0 10 20 30 40 50 60 70 80 90 100 1GS 0GB 2GS 0GB 0GS 1GB 3GS 0GB 1GS 1GB 4GS 0GB 2GS 1GB 0GS 2GB 3GS 1GB 1GS 2GB 4GS 1GB 2GS 2GB 0GS 3GB 3GS 2GB 1GS 3GB 4GS 2GB 2GS 3GB 0GS 4GB 3GS 3GB 1GS 4GB 4GS 3GB 2GS 4GB 3GS 4GB 4GS 4GB Figure 4.4 – Power production limits of the GENSET’s With the knowledge about of all possible GENSET’s makes conceivable the computation of optimal power production for each unit in function of the net load that each GENSET is able to produce. Therefore, in this work, it is proposed the establishment of the concept of equivalent optimal generation unit. As such, each GENSET is considered as an equivalent GENSET 1GS 2GS 3GS 4GS 1GS 1GS 1GS 1GS 2GS 2GS 2GS 2GS 3GS 3GS 3GS 3GS 4GS 4GS 4GS 4GS 0GS 0GS 0GS 0GS 0GB 0GB 0GB 0GB 1GB 2GB 3GB 4GB 1GB 2GB 3GB 4GB 1GB 2GB 3GB 4GB 1GB 2GB 3GB 4GB 1GB 2GB 3GB 4GB Max [kW] 7 200 14 400 21 600 28 800 23 700 40 200 56 700 73 200 30 900 47 400 63 900 80 400 38 100 54 600 71 100 87 600 45 300 61 800 78 300 94 800 16 500 33 000 49 500 66 000 Min [kW] 3 848 7 696 11 544 15 392 12 258 20 668 29 078 37 488 16 106 24 516 32 926 41 336 19 954 28 364 36 774 45 184 23 082 33 212 40 622 49 032 8 410 16 820 25 230 33 640 Power [MW] Chapter 4 - Generation scheduling under uncertainty 106 generator with an equivalent cost function. In figure AI.1 of annex I the specific consumption curves of all GENSET’s are shown. 4.3.3 Equivalent optimal generation unit As shown in figure 4.1a) calculation of equivalent generation unit is done offline for each GENSET merely once, and only if there are changes in the number of thermal power plants, units or in the rated power of the units recalculation is done. There are several methodologies to solve the ED, mainly when the cost functions are continuous and convex [145]. Following figure 4.3 and equations (4.1) and (4.2) it is visible that fuel consumption functions are convex, continuous and differentiable. In case of a GENSET with units of the same type (equal technical limits and cost functions), the power production is divided equitably for all units of the GENSET. If the units have different cost functions an ED must be solved in order to establish the production of each unit. In this work, in order to get a solution to the problem defined by equation (4.3), the ED was solved by the Lagrange multipliers. Later it is proposed a methodology to be used in case of non-convex cost functions. ( ) 1 min max 1 min .. G G N T ii i i ii N iN i F FP st P PP PL = = = ≤≤ = ∑ ∑ (4.3) Minimum and maximum production limits of each unit and the total production to be equal to the net load (LN) requirement are the unique restrictions taken into account in this formulation. For simplicity, the problem is formulated and solved without considering transmission losses, which may be introduced later, together with some changes in the resolution procedure. Equations (4.4) up to (4.7) specify the mathematical formulation for solving the ED of a GENSET, taking for instance the units 4 (GS) and 5 (GB), in function of the net load. The Lagrange function (4.4) is obtained when the equality restriction, equation (4.3), multiplied for an undetermined multiplier λ is added to the objective function [49]: ( ) ( ) 5 45 5 4 4 ,, ii iNi i PP L PFP λλ = =  = +−   ∑∑ L (4.4) where L is the Lagrangian operator and λ the Lagrange multiplier. Chapter 4 - Generation scheduling under uncertainty 107 The derivative of Lagrange function (4.5) gives the necessary conditions for the existence of a minimum cost operating condition for the GENSET. This happens when the incremental cost rates of all units are equal to some undetermined λ [49]. ( ) ( ) ,0 i ii ii P FP PP λλ ∂∂ = −= ∂∂ L (4.5) In order to calculate the most economic production of each unit for each value of net load. Equations systems (4.6) and (4.7), resulting from (4.5), can be solved by any linear programming solver. 4 44 55 5 45 2 01 02 1 1 10 N a Pb aP b L PP λ −−           −=−           − − −−      (4.6) 1 44 4 55 5 45 2 01 02 1 1 10 N Pa b Pa b L PP λ − −−         = −−         − − −−     (4.7) Considering the coefficients ai and bi from equations (4.1) and (4.2) with net load varying between 12258 kW and 23730 kW (limits of the GENSET 1GS_1GB), several values of fuel consumption and power production were obtained. Regarding the net load, in figure 4.5 is shown the power production profile of each unit of the GENSET as well as the equivalent fuel consumption. Fitting the resulting points makes possible the definition of the mathematical functions of each unit production as well as the fuel consumption for each net load and it delivers a linear function for the power production and a second order function to the fuel consumption. In figure 4.5 it can be clearly observed three distinct areas with different power production and fuel consumption functions. For smaller than 15594 kW net load values, the minimum production restriction of unit G4 is activated, remaining its production constant whereas the unit G5 continues to produce the lasting net load. With this, the fuel consumption function is also different. The opposite case happens when the net load is higher than 22425 kW. In this case G5 reaches the maximum production and the remaining load has to be fed by G4. Between the values of net load of 15594 kW and 22425 kW, both units are working inside their production limits. This procedure should be repeated to each GENSET combination and should be done only once and off-line. Chapter 4 - Generation scheduling under uncertainty 114 R 2 = 0,50 0 50 100 150 200 250 300 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 tp(i ) Xbest(i ) Figure 4.10 – Particle cloud behaviour for the variable x1 In case of variable x2 the determination coefficient has an important role, giving more weight to the best fitness particle than to the vertex of the function, making the new central point to move to a coordinate close to the minimum. R² = 0,39 -20 0 20 40 60 80 100 120 140 160 180 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 tp(i) Xbest(i) Figure 4.11 – Particle cloud behaviour for the variable x2 Depending of the function to be optimized and the search space covered by the cloud, the second-order fitness function could be concave indicating a maximum instead of a minimum. In this case it is necessary to calculate the trend point in order to ensure that the central point continues to move towards a minimum. This is done calculating the roots of the concave function subtracted by the best fitness, as indicated in equation (4.12) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 01 2 0 i ii ii bi a ax a x fx++− = (4.12) Then, the closest trend point xb(i) will be chosen using equation (4.13) ( ) ( ) ( ) { ( ) ( ) } 12 min , b i i i i bi p t xxxx=−− (4.13) To prevent a premature convergence into a local optimal and simultaneously to ensure a fast optimum search, the variance of the Gaussian distribution of the particles cloud must be dynamically adapted along the search space. Fitness value Fitness value Chapter 4 - Generation scheduling under uncertainty 115 Fs1 New cloud generation After the definition of a new central point position in iteration (k+1), the cloud is centred in that point, with new variance values for each dimension i calculated by equation (4.14). ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 21 12 2.. i k kk i si s i kFF σσ += (4.14) The changing in cloud distribution variance and consequently in the search space results from equations (4.15) and (4.17) represent both inverted sigmoid functions as shown in figures 4.12 and 4.13. ( ) ( ) ( ) ( ) 2 18 1 1 1 k i si R tc ts F e ϕ ϕ  −−     ∆− = + + (4.15) As seen in figure 4.12, if the value of ( ) 2 i R ϕ is near to 0, the value of Fs1(i) will remain near 1. If ( ) 2 i R ϕ is near to 1 the central point is near to an optimal solution, so the value of Fs1(i) will decrease depending from ∆ϕ(i). This value is dynamic and it is calculated by equation (4.16). ( ) ( ) 1 1. i i h ϕϕ ∆= +∇ (4.16) The slope calculated in (4.10) will define how narrow should be the next cloud, as depicted in figure 4.12. The h parameter of equation (4.16) can be used to speed up the search; however, if the value is too high, the search is faster, but not so accurate. 0.00 0.20 0.40 0.60 0.80 1.00 1.20 00.2 0.4 0.6 0.8 1 Determination coefficient ( ) 2 i R ϕ Figure 4.12 – Inverted sigmoid behaviour for different values of ϕ ∆ On the other hand, the value Fs2 in equation (4.14) intends to expand the cloud’s variance. If the coefficient of determination ( ) 2 fitness i R , calculated in the quadratic regression, has very low values, it means that the search space delimited by the cloud is quite flat and the cloud should ∆ϕ( i ) ∆ϕ(i) ∆ϕ( i ) Chapter 4 - Generation scheduling under uncertainty 116 be increased by a K value. Hence, if ( ) 2 fitness i R is approximately 1, low or no variation will be verified in Fs2, as shown in figure 4.13. ( ) ( ) ( ) ( ) ( ) 2 2 1 1 1 fitness aR b i k si K F e −+ − = − + (4.17) 0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 00,2 0,4 0,6 0,8 1 Determination coefficient of second order polynomial 2 ()fitness i R Figure 4.13 – Inverted sigmoid behaviour with K=1,5 The variance of the particles distribution to each dimension is given by (4.14), allowing the cloud to increase or decrease the space under search by a dynamic and automatic way. Stopping criteria In this algorithm several stopping criteria will be used. One stopping criterion is the definition of a fitness value to be reached, although there are no guarantees that the algorithm will reach the value, leading to a never-ending. On the other hand, the minimum value of the functions is sometimes unknown. For this, a variation is introduced for making the algorithm to run until the fitness value does not change in a certain amount during a determined number of iterations. The objective of this criterion is to decrease the processing time as long as the relation between processing time/results is considered acceptable. However, it can lead to premature stops in flat regions of the search spaces. In case of none of the above criteria could stop the algorithm, there is a defined maximum number of iterations. This criterion is not a guarantee to find the global or even a local optimum before reaching the maximum of iterations. In the proposed algorithm all these criteria can be used to stop the algorithm whenever one of them is met. Fs2 K Chapter 4 - Generation scheduling under uncertainty 117 4.4.2 Evaluation of SCO’s performance To demonstrate the performances and behavior of SCO, among several benchmark functions to test heuristics [146], the Rastrigin’s function (4.18) was used. Although being continuous, it is also multimodal, presenting several local minima which represent a big challenge to optimization algorithms. ( ) ( ) ( ) 22 1 10cos 2 10 ii i i fx x x π = =−+ ∑ (4.18) To an easier analysis, the function has only two dimensions in the range [-5,12 ; 5,12]. For the present example the chosen central point starting values were (5,0;5,0), close to the search space limits and far away from the global minimum f(0,0) = 0 (though in practice the starting values are randomly created). The initial value of variance for each cloud dimension was σ2 = 0,1. The algorithm parameters and related equations are presented in table 4.3. As stopping criteria, a maximum number of 500 iterations or minimum fitness value of 1e-10 were chosen (the remaining stopping criteria were not applied). Table 4.3 – Parameters of the SCO to solve Rastrigin’s function Parameter value Eq. nº particles 50 - h 10 (4.16) t c 0,5 (4.15) t s 1,0 (4.15) K 1,01 (4.17) a 50 (4.17) b 5 (4.17) With a small variance it is intended to create a small cloud far away from the global optimum to create particularly hard conditions and validate the search capacity of SCO. After the first iteration and as an example, the cloud dimensions were comprised between [4,1681;5,1098] for variable x1 and [4,251;5,1172] for variable x2. In figure 4.14 a) and b), it is shown the fitness value evolution together with cloud variance for each dimension. As seen in figure 4.14 b), the size of cloud greatly varies throughout the progress of the search. With an intentionally small initial cloud (σ2 = 0,1), the gathered search space information is limited. So, in early iterations, (less than 100 iterations), the cloud expanded, and the variance of variable x2 reached a value near 6. Chapter 4 - Generation scheduling under uncertainty 118 050 100 150 200 250 300 350 400 450 10 -20 10 -10 10 0 10 10 fitness a) 050 100 150 200 250 300 350 400 450 10 -10 10 0 b) variance σ (x1) σ(x2) a) fitness evaluation and b) cloud’s variance Figure 4.14 – SCO’s performance solving Rastrigin’s function This corresponds to an increase of about 60 times the initial variance, as it can be seen in figure 4.15, which shows the magnification of figure 4.14 b). 010 20 30 40 50 60 70 80 90 100 10-1 100 101 variance σ(x1) σ(x2) Figure 4.15 – Behaviour of cloud variance Successive increases in cloud size allowed the perception of the direction to be taken for reaching the near-optimal global. Thus, when the cloud was positioned above the global minimum, the search was refined contracting the cloud to perform a narrow search, with the objective function reaching the value 1e-9. To reach the goal of 1e-10, the cloud once again changed its size and shape. Once more, the variance took different values to each dimension separately, until it reached the fitness value of 4,93e-12 and stopped the search. In this case it is clear the behavior of SCO and its capacity to find the minimum of fairly difficult problems with wide search space and large number of local minima. 4.4.3 Evaluation of SCO’s in ED To evaluate the SCO’s performances for solving ED, several case studies with non-convex or non-continuous cost functions and non-continuous restriction, as presented in equation (2.14) and (2.15), were tested. To apply SCO to an ED problem, the following 12 steps must be done: Chapter 4 - Generation scheduling under uncertainty 119 Step 1 Create randomly a central particle Pq(i), with i=1..NG dimensions, (each dimension represents a generation unit) and according to its technical limits, as in (4.19). If there are starting values, then ( ) ( ) 0 qii PP= . ( ) ( ) ( ) max min min 0,1 . ii i qi P rand P P P= −+ (4.19) Step 2 Create the remaining cloud with [NG × NP] dimensions, where j = 1..NP particles and i = 1..NG dimensions normally distributed. The remaining particles are normally distributed, centered in the central particle and with standard deviation ( ) ( ) 1k i σ = , as in (4.20). ( ) ( ) ( ) ( ) ( ) ( ) ( ) , ~, k kk ij qi i P NP σ (4.20) Step 3 To each cloud’s particle ( ) ( ) k j P , calculate the transmission losses ( ) k L P by ( ) ( ) ( ) ( ) 11 1 GG G NN N k kk k L i iz z oi i oo iz i P P BP BP B = = = = ++ ∑∑ ∑ (4.21) Step 4 Evaluate each particle with (4.22) and retain the best fitness value ( ) ( ) k j BEST P position. ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ,, 1 1 G GN k kk j ij ij i Nk iL i f FP q P DP = = = + −− ∑ ∑ (4.22) The fitness function (4.22) has a penalty factor q to decrease the deviation between the power production and the sum of power demand and active losses as in (4.4). Step 5 Calculate the second-order regression coefficients to each dimension i, ( ) ( ) ( ) [ ] 12 ii i o βββ and the determination coefficient ( ) ( ) 2k i q R . Step 6 Verify the convexity of polynomials. If it is convex, then ( ) ( ) 2 ki pi i b ta = − . If not calculate ( ) ( ) k pi t by (4.12) and (4.13). Step 7 Generate new central particle by (4.23) and verify if it fulfils all constrains. ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 122 .1. kk k qi pi qi qi iBEST P tR R P + = +− (4.23) If the central particle fulfils all the constraints, then it is a feasible solution (as remain particles act as sensors and are not candidates to a solution, there is no need to satisfy all constrains). Chapter 4 - Generation scheduling under uncertainty 120 Step 8 Calculate the Euclidean distance for each particle j to the central particle, by (4.24) ( ) ( ) ( ) ( ) ( ) ( ) ,, k kk ij qi ij PP ϕ = − (4.24) Step 9 Calculate the linear regression coefficients [α0(i) α1(i)] and the determination coefficient ( ) ( ) 2k i R ϕ . Step 10 Calculate the new standard deviation for each dimension i of new cloud by (4.14). Step 11 If k = itmax go to step 12, otherwise, k=k+1 and go to step 2. Step 12 The central particle which generates the latest best fitness value represents the optimal power generation of each thermal unit and, consequently, the minimum total generation cost. To verify the feasibility of the proposed method, some cases were studied to demonstrate the capacity of the algorithm to reach the optimal values as well as the capacity to solve highly constrained problems with growing dimension. As any other heuristic method it may not converge to exactly the same solution at each run. Due to their stochastic behaviour, their performances could not be judged by the results of a single trial so all cases were performed 50 times keeping the average, maximum and minimum cost values. Mainly, in power systems literature the convergence tests in ED problems are mostly related with number of iterations or generations [33],[36],[41] and [147] or CPU time per iteration/generation [35]. However, this way does not give adequate information about the computational effort to perform a task in order to have the same base of comparison with other techniques [43]. Thus, in this work, to compare the computational efforts independently of the CPU or number of iterations, the number of objective function evaluations is used [148], [149]. The proposed algorithm was implemented in Matlab® (R2010a) and executed on a Core (i2) 1,59 GHz processor. First case study The first test system consists of six thermal units with prohibited operation zones and ramp limits, as shown in tables 4.4 and 4.5, feeding a load of 1263 MW. The network has 26 buses and 46 transmission lines characterized by the losses coefficients matrices Bij, B0j and B00, with 100 MVA capacity base. These matrices are shown in Annex II. This is a very common case study, largely used for comparison of performances between metaheuristics Chapter 4 - Generation scheduling under uncertainty 121 [34],[36],[43],[150] and [151]. The units have cost functions defined by second order continuous and convex functions and the initial values of each unit are defined by 0 i P . Table 4.4 – Generating unit’s data (first case study) Unit P i min (MW) P i max (MW) a i ($/MW2) b i ($/MW) c i ($) UR i (MW/h) DR i (MW/h) P i o (MW) 1 100 500 0.0070 7.0 240 80 120 440 2 50 200 0.0095 10.0 200 50 90 170 3 80 300 0.0090 8.5 220 65 100 200 4 50 150 0.0090 11.0 200 50 90 150 5 50 200 0.0080 10.5 220 50 90 190 6 50 120 0.0075 12.0 190 50 90 110 Table 4.5 – Generating unit’s prohibited zones (first case study) Unit Prohibited zones (MW) 1 [210 240] [350 380] 2 [90 110] [140 160] 3 [150 170] [210 240] 4 [80 90] [110 120] 5 [90 110] [140 150] 6 [75 85] [100 105] In this case simulation, the central particle, as well as each individual particle of the cloud will have 6 dimensions (P1...P6), one for each generation unit. Depending on the number of particles NP the dimension of cloud will be [6 x NP]. The number of evaluations was limited to 5000 and each attempt was ran 50 times as in [34]. As explained above, the particles act as “sensors” of the search space and allow the calculation of the first and second order polynomials, as depicted in figures 4.7 and 4.8. Therefore, there are a minimum number of particles necessary to describe the curve fitting. On the other hand, a large number and the need of evaluation of each particle will slow down the algorithm. After some trials, the number of particles was set to 50. In this case study, the value of h of (4.16) is 1 and the remaining values of the parameters are those presented in table 4.3. The penalty value q in (4.22) was set to 30. Figure 4.16 shows the best solution obtained by SCO for this case study and the results obtained by the PSO and GA, proposed by [34] and New-PSO (NPSO), PSO with Local Random Search (PSO-LRS) and New PSO-LRS (NPSO-LRS) all proposed by [36]. 15400 15420 15440 15460 15480 15500 15520 15540 SCO PSO [3] GA [3] PSO-LRS [10] NPSO [10] NPSOLRS [10] Figure 4.16 – Comparative minimum, maximum and average solutions (first case study) Cost ($/h) Chapter 4 - Generation scheduling under uncertainty 122 In figure 4.17 the convergence behaviour of SCO’s best solution is shown. 0500 1000 1500 2000 2500 3000 3500 4000 4500 5000 1.54 1.545 1.55 1.555 1.56 1.565 1.57 1.575 x 10 4 Number of evaluations Figure 4.17 – Convergence behaviour of SCO for a 6-units problem From figure 4.16, it can be observed that SCO reaches best solutions than the other 5 methods for the minimum and average costs. As in [36], with NPSO and NPSO-LRS, SCO it had a fast convergence to the cost value 15450$, (less than 400 evaluations), reaching a lower value few evaluations after. In addition, the losses obtained by SCO were fewer when compared to the remaining methods. The results are exposed in table 4.6. The algorithm demonstrated good velocity of convergence, reached a lower cost for the generation configuration and had lower losses. Table 4.6 – Results obtained (6-units 1263 MW) (Best individual) Power output (MW) Method SCO PSO [34] GA [34] PSO-LRS [36] NPSO [36] NPSO-LRS [36] PG1 445,3 447,5 474,8 447,4 447,5 447,0 PG2 176,8 173,3 178,6 173,3 173,1 173,4 PG3 265,0 263,5 262,2 263,4 262,7 262,3 PG4 135,3 139,1 134,3 139,1 139,4 139,5 PG5 167,7 165,5 151,9 165,5 165,3 164,7 PG6 85,3 87,1 74,2 87,2 88,0 89,0 PT (MW) 1275,50 1276,01 1276,03 1275,95 1275,95 1275,94 PLoss (MW) 12,50 12,96 13,02 12,96 12,95 12,94 Cost ($/h) 15443 15450 15459 15450 15450 15450 Second case study The second test is an extent of the first. Consisting of 15 thermal units system, whose characteristics are presented in tables 4.7 and 4.8, feeding a load of 2630 MW [33],[34],[43] and [54]. The thermal units are connected to a 30-bus network with active losses matrices shown in annex II. Cost ($/h) Chapter 4 - Generation scheduling under uncertainty 123 Table 4.7 – Generating unit’s prohibited zones (second case study) Unit Prohibited zones (MW) 2 [185 225] [305 335] [420 450] 5 [180 200] [305 335] [390 420] 6 [230 255] [365 395] [430 455] 12 [30 40] [55 65] Table 4.8 – Generating unit’s data (second case study) Unit Pi min (MW) Pi max (MW) ai ($/MW2) bi ($/MW) ci ($) URi (MW/h) DRi (MW/h) P i o (MW) [33] [34] 1 150 455 0.000299 10.1 671 80 120 395 400 2 150 455 0.000183 10.2 574 80 120 450 300 3 20 130 0.001126 8.8 374 130 130 50 105 4 20 130 0.001126 8.8 374 130 130 104 100 5 150 470 0.000205 10.4 461 80 120 426 90 6 135 460 0.000301 10.1 630 80 120 208 400 7 135 465 0.000364 9.8 548 80 120 286 350 8 60 300 0.000338 11.2 227 65 100 262 95 9 25 162 0.000807 11.2 173 60 100 95 105 10 25 160 0.001203 10.7 175 60 100 134 110 11 20 80 0.003586 10.2 186 80 80 67 60 12 20 80 0.005513 9.9 230 80 80 30 40 13 25 85 0.000371 13.1 225 80 80 46 30 14 15 55 0.001929 12.1 309 55 55 15 20 15 15 55 0.004447 12.4 323 55 55 52 20 The parameters of SCO were similar to the first case study. The results obtained were compared with Fast Evolutionary Programming (FEP), Improved Fast Evolutionary Programming (IFEP), Swarm Direction Fast Evolutionary Programming (SFEP) proposed by [33], Particle Swarm Optimization (PSO) and Genetic Algorithms (GA) proposed by [34]. Figure 4.18 shows the total costs and it is clear the improvement of SCO over all other methods; it presents the lowest values for all indicators as well as the lowest average value after 50 trials. 32400 32500 32600 32700 32800 32900 33000 33100 33200 33300 33400 SCO [59] FEP [59] IFEP [59] SFEP [59] SCO [3] PSO [3] GA [3] Minimum Average Maximum Figure 4.18 – Comparative minimum, maximum and average solution (second case study) Although in [33] and [34] the unit’s parameters were similar, the values of 0 i P were different, as shown in table 4.8. This fact is enough to produce fairly different results as it can be observed in table 4.9. These differences are due to the ramps limits associated with starting Cost ($/h) [33] [33] [33] [33] [34] [34] [34]