Impacto do Erro da Previsão Eólica nas Necessidades a Longo Prazo de Reserva Operacional

Faculdade de Engenharia da Universidade do Porto Impacto do Erro da Previsão Eólica nas Necessidades a Longo Prazo de Reserva Operacional João Teixeira VERSÃO FINAL Dissertação realizada no âmbito do Mestrado Integrado em Engenharia Eletrotécnica e de Computadores Energia Orientador: Professor Manuel António Cerqueira da Costa Matos, Ph.D. Coorientador: Leonel de Magalhães Carvalho, Ph.D. Porto, Portugal Outubro, 2015 ii © João Emanuel Leal Teixeira, 2015 iii iv v Earth is abundant with plentiful resources. Our practice of rationing resources through monetary control is no longer relevant and is counter-productive to our survival. Jacque Fresco vi vii Resumo A cada dia que passa, torna-se mais forte a evidência de que as centrais elétricas baseadas em combustíveis fósseis terão que ser substituídas por outras centrais mais sustentáveis, nomeadamente, centrais que usam fontes de energia renovável. Sendo o sistema elétrico um dos que mais afeta o meio ambiente, nomeadamente em termos de emissão de gases de efeito de estufa, tem havido um crescente incentivo por parte dos governos mundiais ao uso destas formas de energias alternativas. A penetração de fontes de energia renovável no SEE trás benefícios e desvantagens. De todas as fontes de energia renovável existentes, é a eólica a mais popular e a que mais energia injeta na rede. No entanto, o comportamento variável do vento trás incertezas adicionais não só para operação mas também para o planeamento do sistema de produção. Para manter o equilíbrio entre a geração e a carga durante a operação do sistema, são necessárias reservas de potência no sistema para fazer face às situações onde existe excesso ou défice de geração. Estas situações resultam não só do comportamento estocástico associado à energia renovável e à carga do sistema mas também da indisponibilidade de geradores que podem sair de serviço intempestivamente por avaria. Mediante esta incerteza crescente, torna-se, portanto, necessário planear e dimensionar adequadamente não só a reserva estática mas também a reserva operacional do ponto de vista estratégico, i.e., no planeamento dos sistemas de produção. O RESERVAS é um programa que avalia a disponibilidade a longo-prazo das reservas estática e operacional do sistema de produção. Baseando-se em Simulação Sequencial de Monte Carlo, consegue estimar os índices de fiabilidade para os dois tipos de reserva. Para se simular a incerteza associada à previsão da potência eólica durante a avaliação de reserva operacional, o programa RESERVAS utiliza o modelo Naïve que permite estimar erros de previsão usando o conceito de persistência. Este modelo é simples e amplamente usado para previsão a curto-prazo. No entanto, existem situações em que este modelo de previsão leva a estimação do erro de uma viii forma grosseira e, consequentemente, não permite aferir a verdadeira adequação da reserva operacional que existe no sistema de produção. De forma a ultrapassar este problema, é proposto um modelo ARIMA para melhorar o erro de previsão e substituir o modelo Naïve. O sistema teste utilizado é baseado no sistema didático IEEE RTS 79 onde se incluiu um parque eólico e variação mensal da capacidade de produção hídrica. O resultado destes testes são os índices de fiabilidade da reserva operacional para vários níveis de integração eólica, ao nível HL1. A este nível hierárquico, todos os geradores e cargas do sistema estão ligados a um barramento. Os índices de fiabilidade analisados são o LOLP, LOLE, EPNS, EENS, LOLF e LOLD. Com diferentes escalonamentos das unidades de produção do sistema, são avaliados vários cenários de capacidade máxima instalada de potência eólica, entre as quais: sem produção eólica, 800, 1600 e 2400 MW. Os escalonamentos das unidades produtoras são ordenados de duas formas: caso A - 1º parque eólico, 2º unidades térmicas, 3º unidades hídricas; caso B - 1º parque eólico, 2º unidades hídricas, 3º unidades térmicas. Em caso de cenário de produção eólica, utilizou-se os modelos de previsão Naïve e ARIMA. A análise destes dois modelos de previsão determinou qual é a influência dos erros de previsão a nível das necessidades de reserva operacional a longo-prazo. No final, e para o mesmo valor de risco, as avaliações obtidas pelo modelo de previsão eólica ARIMA foram melhores do que as obtidas pelo modelo Naïve. Com o modelo Naïve, a Reserva Operacional estava a ser sobrestimada, devido aos erros de previsão, e com esta melhoria no modelo de previsão eólica, conseguiu-se introduzir mais eólica no sistema, pois o erro de previsão diminuiu. O modelo ARIMA é um bom modelo para previsão eólica a curto-prazo, mas introduz erros de sazonalidade que podem vir a influenciar na avaliação dos índices de fiabilidade do sistema. No entanto, os valores obtidos na previsão da potência eólica por este modelo e as respetivas avaliações da Reserva Operacional, são mais satisfatórias do que as obtidas pelo modelo Naïve, dentro do mesmo nível de risco, levando a crer que este é um modelo de previsão que deve substituir o Naïve no programa RESERVAS. ix Abstract Every day becomes stronger the evidence of the fossil fuel-based power plants will have to be replaced by other more sustainable plants, namely, that who use renewable energy sources. Being the electric system one of the most that affects the environment, particularly in terms of emission of greenhouse gases, there has been a growing incentive for world governments to use these forms of alternative energy. The penetration of renewable energy sources in the power system has benefits and also disadvantages. Of all the renewable energy sources that exist, the most popular is the wind energy and that kind of energy is injecting more energy in the network. However, the behaviour of variable wind energy gives additional uncertainties, not only for operation but also for the production system’s management. To keep the balance between generation and load during system operation, power reserves are required in the system to deal with situations that have excess or deficit of generation. These situations are the result not only of the stochastic behaviour associated with the renewable energy and the load of the system, but also of the unavailability of generators that can leave the service releases for malfunction. Through this growing uncertainty, it is therefore necessary to plan and appropriately size not only static but also reserve operating reserve of the strategic point of view, i.e., in the planning of production systems. The RESERVE is a program that evaluates the long-term availability of static and operational reserves of the production system. Based on Sequential Monte Carlo Simulation can estimate the reliability indices for both types of reservation. o simulate the uncertainty associated with the prediction of wind power during the evaluation of operating reserve, the RESERVE program uses the Naïve model which allows to estimate prediction errors using the concept of persistence. This forecasting model is simple and widely used for the short-term forecast. However, there are situations in which this forecast model leads to the estimation of the error of a crudely and consequently does not allow assessing the true operational reserve adequacy that exists on the xvi Figura 3.9 – A. Autocorrelação dos resíduos; B. Correlação parcial dos resíduos.................... 49 Figura 3.10 – Histograma dos resíduos ARIMA(4,1,3) (laranja) e Naïve (azul) ...................... 51 Figura 3.11 – Gráficos Box para os modelos Naïve (esquerda) e ARIMA (direita) ................. 52 Figura 3.12 – Gráfico quantile-quantile dos modelos Naïve (esquerda) e ARIMA (direita) ..... 52 Figura 4.1 – Carga anual do IEEE RTS 79 HW ..................................................................... 56 Figura 4.2 – Produção hídrica mensal vs Carga média mensal ............................................... 57 Figura 4.3 – Parte negativa dos erros de previsão ................................................................. 58 Figura 4.4 – Valores do LOLP horários para potência eólica de 1600 MW............................. 64 Figura 4.5 – Reserva Operacional solicitada na simulação A-1600 MW................................. 64 Figura 4.6 – Valores do LOLP horários para potência eólica de 2400 MW............................. 65 Figura 4.7 – Reserva Operacional solicitada na simulação A-2400 MW................................. 66 Figura A.1 – Diagrama unifilar do sistema IEEE RTS 79 ...................................................... 76 xvii Lista de tabelas Tabela 3.1 – Tabela de valores estatísticos dos dados ............................................................41 Tabela 3.2 – Resultados do teste KPSS para a série eólica .....................................................43 Tabela 3.3 – Resultados do teste KPSS para a série eólica diferenciada ..................................45 Tabela 3.4 – Critérios de ordem e estimação dos parâmetros ARIMA(4,1,3)...........................46 Tabela 3.5 – Valores dos parâmetros para o modelo ARIMA(4,1,3) .......................................46 Tabela 3.6 – Avaliação dos resíduos para o ARIMA(4,1,3)....................................................47 Tabela 3.7 – Teste Ljung-Box para o modelo ARIMA(4,1,3).................................................48 Tabela 3.8 – Avaliação dos resíduos para o Naïve .................................................................49 Tabela 3.9 – Resultado do teste Ljung-Box para o modelo Naïve ...........................................50 Tabela 3.10 – Comparação da avaliação de cada modelo .......................................................51 Tabela 3.11 – Resultado do teste WRS .................................................................................53 Tabela 4.1 – Ordem de mérito das unidades produtoras .........................................................58 Tabela 4.2 – Sumário dos testes ...........................................................................................59 Tabela 4.3 – Parâmetros das simulações ...............................................................................59 Tabela 4.4 – Índices de fiabilidade sem produção eólica ........................................................60 Tabela 4.5 – Avaliação da ROPE para o modelo de previsão ARIMA.......................................61 Tabela 4.6 – Avaliação da ROPE para o modelo de previsão Naïve ..........................................61 Tabela 4.7 – Comparação dos índices de fiabilidade para a ordem de mérito A .......................63 Tabela 4.8 – Comparação dos índices de fiabilidade para a ordem de mérito B .......................66 Tabela A.1 – Informação geradores do sistema IEEE RTS 79 HW ..........................................75 xviii xix Abreviaturas e Símbolos Lista de Abreviaturas ACF Autocorrelation Function ADF Augmented Dickey–Fuller AIC Akaike Information Criteria ARIMA Autoregressive Integrated Moving Average BIC Bayesian Information Criteria EENS Expected Energy not Supplied EPNS Expected Power not Supplied FEUP Faculdade de Engenharia da Universidade do Porto HL0 Hierarchical Level Zero HL1 Hierarchical Level One HL2 Hierarchical Level Two HL3 Hierarchical Level Three INESC TEC Instituto de Engenharia de Sistemas e Computadores - Tecnologia e Ciência xx KPSS Kwiatkowski–Phillips–Schmidt–Shin LOLD Loss of Load Duration LOLE Loss of Load Expectation LOLF Loss of Load Frequency LOLP Loss of Load Probability MAE Mean Absolute Error MAPE Mean Absolute Percentage Error MCS Monte Carlo Simulation ME Mean Error MLE Maximum Likelihood Estimation MSE Mean Square Error MTTF Mean Time to Fail MTTR Mean Time to Repair PACF Parcial Autocorrelation Function PRE Produção em Regime Especial REE Red Eléctrica de España REN Rede Elétrica Nacional RMSE Root Mean Square Error SEE Sistema Elétrico de Energia SMAPE Symmetric Mean Absolute Percentage Error UE União Europeia VAR Vector Autoregression xxi WECS Wind Energy Conversion Systems WRS Wilcoxon rank-sum WTG Wind Turbine Generator Lista de Símbolos p Ordem do modelo AR d Ordem de diferenciação da série temporal q Ordem do modelo MA c Constante do modelo ARIMA(p,d,q) 𝜙𝑝 Parâmetro do modelo AR(p) yt Valor da série temporal a prever no instante t 𝜃𝑞 Parâmetro do modelo MA(q) 𝑒𝑡 Erro da previsão no instante t T Dimensão da série temporal eólica L Valor likelihood Q* Valor do teste Ljung-Box rk Valor da autocorrelação no instante k α Nível de confiança pt Percentagem do erro 𝑒𝑡 no instante t σ Desvio padrão xxii σ2 Variância 𝐸[𝐻(𝑋)] Valores estimados de 𝐸[𝐻(𝑋)] 𝐻(𝑋) Função teste na simulação Monte Carlo 𝑁 Número de estados do sistema 𝑉(𝐸[𝐻(𝑋)]) Variância estimada 𝛽 Coeficiente de variação (critério de paragem no MCS) 𝑆𝑋𝑓 Subconjuntos de estado de insucesso 𝑆𝑋s Subconjuntos de estado de sucesso 𝑁𝑌 Número de períodos simulados 𝜆 Taxa de avaria 𝜇 Taxa de reparação 𝑇𝑈𝑝/Down Tempo de residência em cada estado (Up/Down) 𝐶𝑘 Capacidade de produção da unidade k Δ𝑃𝑗𝑖 Perda de carga do estado i na hora j. 𝜎𝐶𝑃 Desvio standard distribuição gaussiana contínua de incertezas a curto-prazo 𝜎𝐿𝑃 Desvio standard distribuição gaussiana contínua de incertezas a longo-prazo STU Incerteza de curto-prazo LTU Incerteza de longo-prazo PMáx Valor horário da carga 𝑅STA Reserva Estática 𝑅𝑂𝑃𝐸 Reserva Operacional 𝑅S Reserva Secundária xxiii 𝑅T Reserva Terciária 𝛥𝐿 Variação da carga em curto-prazo 𝛥𝑃𝑊 Variação da produção eólica 𝛥𝐺 Variação da capacidade de produção p-value Valor p resultante de testes estatísticos χ2α,h-K Valor Chi-square com nível de confiança α e h-K graus de liberdade z-value Valor z resultante do teste Wilcoxon rank-sum xxiv Introdução 1.1 Motivação e Contexto da Dissertação Ao longo dos últimos anos, o Sistema Elétrico de Energia (SEE) está a mudar a sua matriz energética substituindo centrais que usam combustíveis fósseis por alternativas renováveis, sendo a Alemanha um bom exemplo internacional [1]. Sendo o sector energético um dos mais críticos no que toca ao impacto ambiental, esforços estão a ser feitos para que este problema tenha uma solução. O objetivo “20-20-20”, diretiva da UE (União Europeia) é um exemplo de mudança e preocupação ambiental [2]. Dentro desta diretiva estão os seguintes objetivos a serem atingidos até ao ano 2020, dentro do tema da alteração ambiental e energia sustentável:  Reduzir 20% a emissão dos gases de efeito de estufa de 1990;  20% da energia ser proveniente de energias renováveis;  Aumentar 20% a eficiência energética. Recentemente, os líderes industriais G7 acordaram em reduzir os gases de efeito de estufa por fase, cortando o uso de combustíveis fósseis no final deste século [3], reduzindo de 40% a 70% até 2050 a emissão dos gases de efeito de estufa, tendo como referência os valores de 2010. Um dos países que tem atraído as atenções dentro deste contexto é a Dinamarca. Este país tem sido pioneiro na integração de energia proveniente de fontes renováveis, maioritariamente energia eólica. E recentemente conseguiu um valor de produção de energia proveniente de fonte renovável que atingiu os 140% do valor do pico de carga [4], mostrando que é possivel efetuar a mudança do uso de combustíveis fósseis por soluções renováveis. A par deste país estão o Reino Unido e 8 Estado da Arte aplicaram várias vertentes do modelo ARIMA. Este trabalho centra-se na aplicação e comparação destes dois modelos, onde serão falados mais detalhadamente nas subsecções 2.1.4 e 2.1.5. Um dos métodos de previsão são as Redes Neuronais Artificiais. Inspirada na biologia dos neurónios do sistema nervoso, este modelo tem como elemento base o neurónio, a unidade de processamento elementar. Nas Redes Neuronais Artificiais também existem ”neurónios” e são eles que permitem complexas relações não lineares entre a variável de resposta e as suas entradas (inputs) [10], [11]. A arquitetura das Redes Neuronais Artificiais baseia-se na organização dos neurónios por camadas, com ligações unidirecionais, em sentido progressivo, dos inputs para as respostas (outputs). As redes mais simples contêm camadas escondidas e são equivalentes a uma regressão linear, obtendo resultado através de uma combinação linear dos inputs. As ligações entre neurónios estão associadas a pesos que são selecionados usando um algoritmo de aprendizagem que minimiza uma dada função de custo, por exemplo, MSE (Mean Square Error) [11]. Se quisermos assumir que eventos futuros influenciam as previsões, podemos usar um modelo de previsão designado modelo de regressão dinâmico. O efeito das férias, mudança de leis, a economia geral, atividade concorrente ou outra variável externa pode explicar uma variação histórica e, com isso, permitir uma previsão mais exata. Este modelo tem como base o modelo ARIMA, apesar de ter informações sazonais que influenciam as variáveis de entrada, e consequentemente, os valores de previsão [11]. O modelo de Vetores Autorregressivos (VAR) tem na sua base influências entre variáveis do vetor de regressão. Por exemplo, variáveis como as despesas de consumo pessoal e rendimento pessoal são interdependentes, e para englobar as dependências lineares entre várias séries temporais, é criado um modelo autorregressivo multivariado. As alterações de despesas de consumo pessoal com base no rendimento pessoal é um exemplo simples de aplicação deste VAR. As despesas são afetadas com o rendimento pessoal e vice-versa, portanto, uma variação numa variável influencia a outra [11]. Existem outros tipos de modelos de previsão, tal como a Suavização Exponencial (Exponential Smoothing) [12] e a previsão de séries temporais hierárquica ou agrupadas [13]. A primeira baseia-se, em grande parte, num modelo ARIMA mais primitivo, o que é interessante no ponto de vista de previsão para séries temporais. Já a segunda cria uma desintegração da série temporal por níveis ou hierarquias. Introdução 9 O modelo Hierárquico pode-se aplicar à previsão de vendas de bicicletas de uma dada loja. As bicicletas podem ser divididas em várias camadas: montanha, estrada, híbrida ou para crianças; a das crianças pode ainda ser dividida por idades. Também se podem agrupar as bicicletas por género: masculinas, femininas ou unissexo. A partir desta desintegração, é possível fazer a previsão hierárquica utilizando um modelo adequado para cada hierarquia [11]. De entre estes modelos, os mais populares e aplicados a séries temporais eólicas são o Naïve, o ARIMA e as Redes Neuronais Artificiais, que serão referenciadas nos seguintes parágrafos pela sua aplicação em alguns trabalhos. 2.1.2 Modelos de Previsão Aplicados à Produção Eólica Tal como dito na subsecção anterior, os modelos mais utilizados na previsão de séries temporais são o modelo ARIMA e o modelo de Redes Neuronais Artificiais. Nos seguintes parágrafos, estarão alguns exemplos de aplicação divididos pelos dois modelos de previsão, e outros modelos de previsão eólica. 2.1.2.1 ARIMA Um dos trabalhos que incluía o uso do ARIMA para previsão da produção eólica, foi proposto por R. Billinton, et al, em [14] e depois completado pelos mesmos, em [15], onde utilizaram um modelo ARIMA para simular a velocidade do vento e analisar os índices de fiabilidade da reserva estática com os dados simulados, para estudar o efeito que os aerogeradores teriam no sistema produtor. Outros modelos foram aplicados para previsão do vento, principalmente, a sua velocidade. Kavasseri and Seetharaman, em [16], usaram um ARIMA fracionário (f-ARIMA) para prever as velocidades do vento para o dia seguinte, que mostram bastantes melhorias na previsão, comparando com o método Naïve. Palomares-Salas, et al, em [17], compararam a performance das Redes Neuronais Artificiais com ARIMA e concluíram que o modelo ARIMA é melhor para intervalos de curto-prazo (10 minutos, 1 hora, 2 horas e 4 horas). Os autores Cadenas e Rivera, em [18], também compararam Redes Neuronais Artificiais com ARIMA sazonal (s-ARIMA). O resultado foi que o modelo sARIMA é melhor para prever velocidades do vento e tem melhor resposta a mudanças estruturais do regime de vento. Brown, Katz e Murphy, [19], apresentaram um modelo Autorregressivo com base nas propriedades dos dados da velocidade do vento: autocorrelação, distribuição não-Gaussiana e 10 Estado da Arte série do vento diurna não estacionária. O seu objetivo foi tornar os dados estacionários no período diurno e aproximar a uma distribuição Gaussiana melhorando, deste modo, a previsão da velocidade do vento. Já Chen, et al, em [20], propuseram um modelo ARIMA tendo em conta os dados da velocidade do vento não estacionários bem como os limites estocásticos da geração eólica. O modelo introduz limites e carateriza a geração eólica estocástica por nível, correlação temporal e condução do ruído. Na literatura existem alguns estudos que incorporam outros parâmetros relacionados com a previsão da velocidade do vento. Um desses estudos desenvolvido por Erdem e Shi, [21], comparam vários modelos ARIMA que consideram a direção e velocidade do vento, com o objetivo de prever esses mesmos parâmetros. Erdem, Shi e Peng, em [22], voltaram a propor um novo método, desta vez usando um modelo ARMA que incorpora direção e velocidade para previsão da potência do vento. Com este modelo obtêm-se ligeiras melhorias na estimação para a velocidade do vento a curto-prazo. 2.1.2.2 Redes Neuronais Artificiais Em [23], Kariniotakis, et al, desenvolveram um modelo baseado em Redes Neuronais Artificiais para a previsão da potência gerada dum parque eólico, mostrando que é melhor do que o clássico Naïve. Com base nesses resultados, foi implementado num controlo de sistema avançado para operação e gestão de um SEE real. Li e Shi, [24], demonstraram que a performance dos três tipos de Redes Neuronais Artificiais na previsão da velocidade de vento na seguinte hora pode variar consoante a métrica de avaliação. Além disso, a seleção do tipo de Redes Neuronais Artificiais para melhorar desempenho na previsão depende da qualidade dos dados. Também usando Redes Neuronais Artificiais, Barbounis e Theocharis, estudaram o problema de previsão a longo-prazo da velocidade e potência do vento. Comparando três tipos de Redes Neuronais Artificiais: infinite impulse response multilayer perceptron (IIR-MLP), local activation feedback multilayer network (LAFMLN), e diagonal recurrent neural network (RNN), chegaram a resultados que demonstram que os modelos treinados pelos métodos sugeridos exibiam melhorias face ao método de persistência [25]. 2.1.2.3 Outros Modelos de Previsão Eólica Monfared, et al, em [26], propuseram uma nova estratégia na previsão de velocidade de vento baseado na lógica fuzzy e redes neurais artificiais. Este método, não só fornece significativamente Introdução 11 menos regras de base, mas também aumenta a precisão de velocidade estimada de vento quando comparado com o tradicional. Já Slougther, et al, [27], usaram a abordagem de modelagem Bayesiana para fornecer os conjuntos na abordagem de previsão probabilística para a velocidade do vento. Ao contrário da previsão da velocidade média do vento, estes desenvolveram uma metodologia para estimar a velocidade do vento máxima analisando a natureza da topologia do sítio e a época do ano. Concluíram que a velocidade pode ter diferentes propriedades, que exigem considerações adicionais. 2.1.3 Influência da Previsão Eólica nos Índices de Fiabilidade Ao longo dos anos surgiram vários estudos que conjugam modelos de previsão eólica com os cálculos dos seus índices de fiabilidade nos respetivos sistemas inseridos no estudo. Billinton, et al, em [15], usando o método de Simulação de Monte Carlo para a velocidade do vento, desenvolveram um estudo capaz de incluir WECS (Wind Energy Conversion Systems) na avaliação dos índices de fiabilidade. Usando um modelo ARMA para simulação a longo-prazo da velocidade do vento, considerando a sua localização, Karki e Hu em [28], não só analisaram os índices de fiabilidade como também os custos da inclusão da energia eólica nos tradicionais sistemas de geração. Mais tarde, e juntando-se a Billinton, [29], desenvolveram um modelo simplificado de geração eólico para avaliação dos índices de fiabilidade. À semelhança dos dois últimos, analisaram os índices de fiabilidade com a inclusão de energia eólica no sistema produtor, preocupando-se com a exatidão da previsão da velocidade do vento. Normalmente utiliza-se um método de Simulação Sequencial de Monte Carlo MCS ou uma representação multiestados do parque eólico, no entanto, o modelo apresentado simplifica o método de avaliação. O método é simplificado através da determinação da mínima representação multiestado da geração, na avaliação dos índices de fiabilidade. Este é aplicável a múltiplas localizações geográficas e adequado para avaliação dos índices de sistemas de potência com significativa penetração eólica. Huang e Billinton [30], estudaram o efeito na probabilidade, frequência e duração média das três fases de estados de operação quando adicionados a potência eólica a um sistema, ou substituindo algumas capacidades de produção convencionais existentes por energia eólica. Para a integração de grandes partes da geração eólica em sistemas de energia, há a necessidade do desenvolvimento de novos algoritmos e ferramentas de previsão para a tomada de decisões no domínio operacional, tendo em conta das incertezas de previsão eólica. Uma dessas decisões 12 Estado da Arte considera os requisitos de reserva operacional disponíveis para atender às variações do vento e da carga. Matos e Bessa em [31], apontaram para esse risco de decisões, mostrando que é possível descrever as consequências de cada nível de reserva, através de um conjunto de índices de risco útil para tomada de decisão durante a operação do sistema. 2.1.4 Autoregressive Integrated Moving Average Este foi o modelo escolhido para o estudo de previsão da produção de um parque eólico, quando aplicado a sua respetiva série eólica temporal [9], [11], [32], com o objetivo de substituir o modelo Naïve no programa RESERVAS. O modelo é na realidade a junção de dois modelos de previsão, decompondo-se em : AR – autoregressive; MA – moving average 1 2.1.4.1 Autoregressive Models O modelo Autorregressivo (AR) faz uma regressão linear dos valores passados da série temporal, com o objetivo de prever o seu valor futuro. Esta regressão é feita numa ordem p, AR(p), em que p é a ordem do atraso na regressão à série temporal. Este modelo pode ser descrito do seguinte modo 𝑦𝑡=𝑐+𝜙1𝑦𝑡−1+𝜙2𝑦𝑡−2+⋯+𝜙𝑝𝑦𝑡−𝑝+𝑒𝑡, (2.1) onde 𝑦𝑡 representa o valor a prever, 𝑦𝑡−𝑝 são os valores passados da série até à ordem de atraso p, e 𝜙𝑝 os seus pesos respetivos. A constante é representada por 𝑐 e 𝑒𝑡 é o erro de previsão para o instante t [9], [11], [32]. Idealmente, este erro de previsão deverá ser ruído branco. O modelo é bastante flexível para lidar com uma certa gama de diferentes padrões de séries temporais, isto quer dizer que, o modelo pode ser ajustado para várias séries temporais, com a mesma ordem de regressão p, mas com pesos 𝜙𝑝 diferentes. 1 Apesar de não estar referenciado, o Intregrated não é um modelo de previsão, mas sim a diferenciação da série a que é aplicado. Na subsecção Série Estacionária e não Estacionária será falado em maior detalhe. Introdução 13 Figura 2.1 – Dois diferentes exemplos de modelos autorregressivos [11] Na Figura 2.1 ilustram-se dois gráficos de modelos autorregressivos diferentes, tanto em ordem p, como em peso 𝜙𝑝, onde ambos apresentam erros de previsão 𝑒𝑡 com distribuição normal, média zero e variância um (ruído branco). 2.1.4.2 Moving Average Models Ao contrário do modelo autorregressivo, o modelo de média móvel usa erros de previsão passados num modelo de regressão tipo. A ordem de regressão é representada por q e genericamente presenta-se MA(q) da seguinte forma 𝑦𝑡=𝑐+𝑒𝑡+𝜃1𝑒𝑡−1+𝜃2𝑒𝑡−2+⋯+𝜃𝑞𝑒𝑡−𝑞, (2.2) onde 𝑐 é a constante e 𝑒𝑡 é o ruído branco, ou erro da previsão, 𝑒𝑡−𝑞 são erros de previsão no instante t-q, sendo q a ordem de regressão de MA [9], [11], [32]. Na realidade, não é possivel observar os valores de 𝑒𝑡, portanto, este tipo de modelo não é um modelo regressivo. Figura 2.2 – Dois diferentes exemplos de modelo de média móvel [11] 14 Estado da Arte Na Figura 2.2 verificam-se dois modelos de média móvel de ordem q e diferentes pesos 𝜃𝑞. A aplicação destes modelos para a mesma série temporal resulta em dois resultados distintos em valor. 2.1.4.3 Série Estacionária e não Estacionária Uma série pode ser estacionária se as suas propriedades não dependerem do tempo em que a série é observada. No caso de uma série temporal, como os dados de velocidade do vento, não se pode afirmar que a série seja estacionária, visto que, a velocidade do vento não estaciona para nenhum valor. Se a série é não estacionária, recorre-se à diferenciação na tentativa de a estabilizar. A diferenciação consiste na diferença entre dois valores consecutivos da série, tal como se pode visualizar na seguinte equação 𝑦𝑡′=𝑦𝑡−𝑦𝑡−1 . (2.3) Na equação (2.3), [32], a ordem de diferenciação d é primeira (d=1). Na prática, raramente é necessário recorrer além da segunda ordem de diferenciação (d=2), que se expressa do seguinte modo 𝑦𝑡′′=𝑦𝑡′−𝑦𝑡−1 ′=(𝑦𝑡−𝑦𝑡−1)−(𝑦𝑡−1−𝑦𝑡−2)=𝑦𝑡−2𝑦𝑡−1+𝑦𝑡−2 . (2.4) Em alguns casos, existe dificuldade em perceber se uma dada série é ou não estacionária. Uma das formas para se verificar se a série é estacionária, consiste na análise das suas autocorrelações. O valor das autocorrelações entre yt e yt+k é obtido de acordo com a expressão [9] 𝑟𝑘=𝑐𝑘 𝑐0 , (2.5) onde 𝑟𝑘 é o valor da autocorrelação para o atraso k, sendo que este atraso varia de 0 e o último atraso (K). A variável 𝑐0 é o valor da variância da série temporal, e 𝑐𝑘 é obtido por 𝑐𝑘=1 𝑇−1∑(𝑦𝑡−𝜇)(𝑦𝑡−𝑘−𝜇) 𝑇−𝑘 𝑡=1 , (2.6) onde T é o número de elementos da série, 𝑦𝑡 e 𝑦𝑡−𝑘 são respetivamente os valores da série no instante t e t-1, e 𝜇 é a média da série. Introdução 15 Figura 2.3 – Autocorrelação de uma série temporal não estacionária 2 A Figura 2.3 mostra um gráfico de autocorrelação, onde se procura encontrar as seguintes caraterísticas referentes as séries não estacionárias:  O gráfico ACF (Autocorrelation function) da série decai lentamente para zero;  O primeiro valor da autocorrelação r1 é normalmente alto (≈1) e positivo. Na Figura 2.3, a linha tracejada representa o valor crítico da autocorrelação dado por ±1,96 √𝑇 ⁄ (2.7) sendo 𝑇 o período de observações. Pela Figura 2.3 identifica-se de que a série é não estacionária, devido a respeitar as últimas caraterísticas. No entanto, existem testes mais eficazes para se determinar se uma série é ou não estacionária. Estes testes são baseados em testes de hipóteses estatísticos sendo o teste Augmented Dickey-Fuller (ADF) um dos mais populares [11]. O teste ADF baseia-se na hipótese nula de que a serie é não estacionária. Se o valor de p for acima do valor crítico (normalmente 0,05), indica que a série é não estacionária sendo que a diferenciação é necessária. Outro teste bastante popular é o de Kwiatkowski-Phillips-Schmidt-Shin (KPSS). Este assume a hipótese nula de a série ser estacionária. Neste caso para valores de p mais baixos que 0,05 sugere que a série não é estacionária, logo a sua diferenciação é necessária. Ambos os testes também admitem que as séries são não sazonais [11]. 2 Gráfico obtido pela autocorrelação de uma série eólica no programa R. 16 Estado da Arte 2.1.4.4 ARIMA(p,d,q) O modelo ARIMA é a junção do modelo AR e MA, com as respetivas ordens p e q, e de grau de diferenciação d, como dito anteriormente. O modelo pode ser interpretado de acordo com as seguintes caraterísticas ARIMA(0,0,0) White noise 𝑦𝑡=𝑒𝑡 (2.8) ARIMA(0,1,0) sem constante Random Walk 𝑦𝑡=𝑦𝑡−1+𝑒𝑡 (2.9) ARIMA(0,1,0) com constante Random Walk with drift 𝑦𝑡=𝑐+𝑦𝑡−1+𝑒𝑡 (2.10) ARIMA(1,0,0) AR(1) 𝑦𝑡=𝑐+𝜙1𝑦𝑡−1+𝑒𝑡 (2.11) ARIMA(0,0,1) MA(1) 𝑦𝑡=𝑐+𝑒𝑡+𝜃1𝑒𝑡−1 (2.12) ARIMA(1,1,1) 𝑦𝑡=𝑐+𝑦𝑡−1+𝜙1𝑦𝑡−1+𝑒𝑡+𝜃1𝑒𝑡−1 (2.13) Consoante a série temporal que se analisa, são ajustadas as suas respetivas ordens p e q, e os seus coeficientes 𝜙 e 𝜃 (a ordem d é determinada pelos testes referidos em 2.1.4.3). Para se entender melhor os modelos ARIMA temos as seguintes hipóteses:  Se 𝑐=0 e 𝑑=0, a previsão a longo-prazo decresce para zero;  Se 𝑐=0 e 𝑑=1, a previsão a longo-prazo tende uma constante diferente de zero;  Se 𝑐=0 e 𝑑=2, a previsão a longo-prazo segue uma linha reta;  Se 𝑐≠0 e 𝑑=0, a previsão a longo-prazo tende para a média da série;  Se 𝑐≠0 e 𝑑=1, a previsão a longo-prazo segue uma linha reta;  Se 𝑐≠0 e 𝑑=2, a previsão a longo-prazo segue uma tendência quadrática. Usando gráficos da autocorrelação e correlação parcial, pode-se ter uma ideia de quais os valores de p e q mais adequadas para se ajustar aos dados temporais. A autocorrelação mede a relação entre 𝑦𝑡 e 𝑦𝑡−𝑘 para diferentes valores de k. E segue a lógica que se 𝑦𝑡 e 𝑦𝑡−1 têm correlação, então 𝑦𝑡−1 e 𝑦𝑡−2 também devem ter correlação. Mas 𝑦𝑡 e 𝑦𝑡−2 podem estar correlacionados, simplesmente por estarem conectados a 𝑦𝑡−1, em vez de qualquer nova informação contida em 𝑦𝑡−2 que pode influenciar a previsão de 𝑦𝑡 [32]. Introdução 17 Usando as correlações parciais, são retirados os efeitos dos atrasos 1 até k-1 das medições das relações entre 𝑦𝑡 e 𝑦𝑡−𝑘. Figura 2.4 – Correlação parcial de uma série temporal não estacionária A Figura 2.4, que representa a correlação parcial da mesma série temporal da Figura 2.3, que mostra uma clara diferença de correlações parciais entre os diferentes atrasos. Neste exemplo, apenas quatro atrasos são significativos, ou seja, estão fora da banda do limite crítico da série. Esta banda está representada por duas linhas azuis tracejadas, e que respeitam a equação (2.7). Devido a T ser elevado, a banda dos limites críticos fica muito próxima da linha horizontal. Pelo gráfico das correlações parciais, pode-se estimar a ordem do modelo Autorregressivo, ou por outras palavras, ter a ideia de quais são os atrasos mais significativos. Apesar de haver algoritmos mais eficazes para deteção dos valores dos coeficientes 𝜙𝑝, pelas correlações parciais percebe-se a influência dos atrasos k no modelo Autorregressivo. 2.1.4.5 Estimação e Ordem do Modelo ARIMA Um dos meios para determinar a ordem de um modelo ARIMA é o critério da informação Akaike (AIC - Akaike Information Criterion), [11], [32], e é definido pela equação 𝐴𝐼𝐶 = −2log(𝐿)+2(𝑝+𝑞+𝑘+1), (2.14) 24 Estado da Arte  H0: a distribuição de ambas as populações é igual;  HA: os valores das duas populações tendem a exceder os outros: 𝑃(𝑋>𝑌)+0.05∙ 𝑃(𝑋=𝑌)>0.05. Para construir o ranking, todas as amostras são juntas e ordenadas por valor. Ao serem ordenadas, é atribuído a cada uma, uma posição (rank). São então separadas e cada valor terá a sua posição no rank. Calcula-se estatisticamente U para cada valor das populações. em que 𝑚 𝑒 𝑛 são os respetivos tamanhos das populações, 𝑟𝑚,𝑛 é o somatório do respetivo ranking. Este teste segue uma distribuição Z, quando a quantidade das amostras é superior a 20. Depois de calculados os U, escolhe-se qua o que teve menor resultado de acordo com O valor de z pode ser calculado por em que ou seja, 𝜇𝑚 e 𝜎𝑚 são respetivamente a média e o desvio padrão da população m (caso seja a população com menor U). Quando determinado o valor z e para um nível de confiança de 95%, o teste indica a possibilidade de rejeição, ou não, da hipótese nula de relação entre os resíduos de cada modelo, se aplicado para dois modelos de previsão, como no caso do modelo de previsão ARIMA e Naïve. 𝑈𝑚=𝑚∙𝑛+𝑚(𝑚+1) 2−𝑟𝑚, (2.28) 𝑈𝑛=𝑚∙𝑛+𝑛(𝑛+1) 2−𝑟𝑛, (2.29) 𝑈=min(𝑈𝑚,𝑈𝑛). (2.30) 𝑧=𝑈−𝜇𝑚 𝜎𝑚, (2.31) 𝜇𝑚=𝑚(𝑚+𝑛+1) 2, (2.32) 𝜎𝑚=√𝑚∙𝑛(𝑚+𝑛+1) 12 , (2.33) Métodos de Simulação 25 2.2 Métodos de Simulação O método de Simulação de Monte Carlo (MCS – Monte Carlo Simulation) pode ser uma ferramenta para resolver tanto problemas estocásticos como determinísticos e as técnicas deste método são aplicados em diversos campos, tais como: cálculos matemáticos complexos, processo de simulação estocásticos, estatísticas médicas e análise de sistemas de engenharia [40]. O interesse ampliou-se devido a fatores como o incremento do poder computacional ao dispor dos utilizadores a preço acessível e, com o desenvolvimento de técnicas de aceleração de convergência, que em algumas centenas de simulações, se pode atingir resultados com a mesma precisão teórica que sem essas mesmas técnicas [41]. O método MCS pode ser visto como uma agregação de estados, incluindo o estado de carga e, por este motivo, é vantajoso face à dificuldade de encontrar um modelo analítico que englobe todas as dependências de funcionamento entre os diferentes elementos do sistema. Devido à sua flexibilidade, têm sido amplamente utilizados em estudos de sistemas. Dependendo das diferentes representações do método de Simulação de Monte Carlo, o mesmo é classificado por: sequencial, quando os estados do sistema são amostrados tendo em conta a cronologia dos acontecimentos; não-sequencial, se o sistema tiver uma representação espacial, ou seja, ser tirada uma “fotografia” ao estado do sistema; sistemas híbridos, tais como o pseudo sequencial [42], pseudo-cronológico [43] e quasi-sequencial [44], sendo que estes não adotam uma representação espacial pura do estado do sistema, nem uma representação cronológica pura [45]. No método de Simulação não Sequencial de Monte Carlo, estima-se 𝐸[𝐻(𝑋)] a partir de um sorteio aleatório de 𝑁 estados do sistema. 𝐸[𝐻(𝑋)]=1 𝑁∑𝐻(𝑥𝑖) 𝑁 𝑖=1 , (2.34) em que 𝐸[𝐻(𝑋)] são valores estimados de 𝐸[𝐻(𝑋)]; 𝐻(𝑋) é a função teste, por exemplo, um índice de fiabilidade (LOLP – Loss of Load Probability); 𝑁 é o número de amostras; 𝑥𝑖 o i-ésimo estado de sistema sorteado [41]. Como em qualquer processo de amostragem, o seu valor médio distribui-se à volta do valor “verdadeiro”, de tal forma que se pode representar a incerteza na estimativa de acordo com a variância do estimador 26 Estado da Arte 𝑉(𝐸[𝐻(𝑋)])=𝑉(𝐻(𝑋)) 𝑁. (2.35) No entanto, a verdadeira variância 𝑉(𝐻(𝑋)) não é conhecida, pelo que deverá também ser estimada através de 𝑉(𝐸[𝐻(𝑋)])=1 𝑁−1∑(𝐻(𝑥𝑖)−𝐸[𝐻(𝑋)]) 𝑁 𝑖=1 , (2.36) sendo 𝑉(𝐸[𝐻(𝑋)]) a variância estimada, e onde 𝑥𝑖 é um estado do sistema de N exemplos. A convergência ou critério de paragem no método de Simulção de Monte Carlo é baseada num coeficiente de variação 𝛽, definido matematicamente por 𝛽=√𝑉(𝐸[𝐻(𝑋)]) 𝐸[𝐻(𝑋)] . (2.37) Figura 2.6 – Exemplo de Simulação não Sequencial de Monte Carlo A Figura 2.6 mostra um exemplo de Simulação não Sequencial de Monte Carlo com o objetivo de estimar o valor de π, onde os pontos azuis são as variáveis de sucesso, e as vermelhadas de insucesso. Neste exemplo, foram simulados pontos com coordenadas compreendidas entre zero e um, onde se pretendia descobrir quais os que se encontravam dentro (sucesso) ou fora (insucesso) de um quarto de área de uma circunferência de raio 1. Métodos de Simulação 27 2.2.1 Método de Simulação não Sequencial de Monte Carlo O método MCS não sequencial, também conhecido por não cronológico, é caraterizado por retirar amostras aos estados do sistema através de “fotografias” do seu comportamento estocástico. Ou seja, as amostras analisadas não têm qualquer dependência temporal dos estados anteriores do sistema. Tal como na Figura 2.6, cada estado está localizada no espaço de estados, não sendo o tempo um fator a ser considerado. O índice de confiança será dado matematicamente pela equação (2.34), sendo 𝐻(𝑥𝑖) é o resultado da função de teste 𝐻 do estado do sistema 𝑥𝑖. Um exemplo desta função aplicada ao índice LOLP (Loss Of Load Probability) [46] 𝐻𝐿𝑂𝐿𝑃(𝑥𝑖)={1 𝑠𝑒 𝑥𝑖∈𝑆𝑋𝑓 0 𝑠𝑒 𝑥𝑖∈𝑆𝑋𝑠, (2.38) em que 𝑆𝑋=𝑆𝑋𝑓∪𝑆𝑋𝑠 é o conjunto de todos os estados do sistema, dividido em dois subconjuntos de estado de sucesso, 𝑆𝑋𝑠, e de insucesso, 𝑆𝑋𝑓 [45]. Devido à sua dificuldade de modelização temporal, a Simulação não Sequencial de Monte Carlo tem dificuldades em modelizar sistemas como os de carga ou variações eólicas. Mesmo assim, estão a ser feitos esforços para ultrapassar esta limitação. 2.2.2 Método de Simulação Sequencial de Monte Carlo O MCS sequencial, ou cronológico, é baseado em modelos estocásticos dos componentes e carga do sistema. Neste método são simulados os estados dos diferentes componentes, num dado período de tempo [47]. É “traçada” a história de vida de um só componente do sistema, garantindo que dois estados consecutivos diferem entre si [45]. Como a sua representação cronológica pode reproduzir a operação do sistema torna-se, por isso, um método muito adequado para avaliação de processos cronológicos, tais como energias renováveis e incertezas naturais (hidrológicas e eólicas), bem como variação de cargas. 28 Estado da Arte Figura 2.7 – Exemplo de uma representação sequencial durante o período T [45] Na Figura 2.7 observa-se uma ilustração do MCS sequencial para a avaliação da capacidade de produção disponível. O estado da carga altera-se entre dois estados consecutivos, estando estes separados por intervalos de tempo. A área sombreada representa os estados onde a capacidade de produção disponível é inferior à carga, ou seja, estados onde existe corte de energia devido a falhas de geradores. Adequando a equação (2.34) ao MCS sequencial, a estimação dos índices de confiança é obtida de acordo com 𝐸[𝐻(𝑋)]=1 𝑁𝑌∑𝐻(𝑦𝑖) 𝑁𝑌 𝑖=1 , (2.39) sendo neste caso 𝑁𝑌 o número de períodos simulados, 𝑦𝑖 a sequência de estados simulados para o período 𝑖 definido por 𝑦𝑖={𝑥𝑛}𝑛=1 𝑆𝑖= {𝑥1,…,𝑥𝑆𝑖} ,𝑆𝑖∈ℕ. (2.40) Aplicando a função de teste ao índice LOLP, define-se como 𝐻𝐿𝑂𝐿𝑃(𝑦𝑖)=1 𝑇∑[𝑑( 𝑥𝑛)×𝐻𝐿𝑂𝐿𝑃(𝑥𝑛)] 𝑆𝑖 𝑛=1 . (2.41) À semelhança do que foi definido anteriormente, 𝑥𝑛 é o estado 𝑛 da sequência, 𝑇 a duração do período simulado (normalmente um ano – 8760h), 𝑑( 𝑥𝑛) é a duração em cada estado e 𝐻𝐿𝑂𝐿𝑃(𝑥𝑛) é o resultado da função de teste para o estado 𝑥𝑛, tal como definido na equação (2.38) em [8], [45], [48]. Métodos de Simulação 29 2.2.3 Modelização dos Componentes do Sistema A capacidade dos componentes do sistema depende de dois modelos distintos:  O modelo do ciclo estocástico falha/restauro;  O modelo de dependência temporal da capacidade do componente. O primeiro modelo está relacionado com o comportamento estocástico das falhas furtuitas dos componentes. O modelo define a disponibilidade do componente, isto é, a transição entre o estado de funcionamento (up state) e de não funcionamento (down state). A disponibilidade do componente pode ser vista como a capacidade máxima que este tem num determinado estado. O modelo do ciclo estocásticos falha/restauro são modelos de Markov de dois estados e de multiestados [7]. Um modelo de Markov de dois estados é utilizado para modelizar os estados dos componentes do sistema. Estes estados variam entre o estado de funcionamento de não funcionamento, onde as taxas de avaria λ e de reparação 𝜇=1 𝑟, são as transições entre os estados up/down e down/up respetivamente, em que r é o tempo médio de reparação, em horas. Figura 2.8 – Modelo de Markov de dois estados [45] A Figura 2.9 mostra um ciclo estocástico falha/restauro (modelo de Markov de dois estados) tem a sua disponibilidade máxima quando está no estado de funcionamento (up state). No estado em de não funcionamento (down state), a capacidade do componente é zero. Assumindo que a duração em cada estado segue uma distribuição exponencial, o tempo de residência em cada estado é dado de acordo com 𝑇𝑈𝑝=−1 𝜆ln𝑈1, (2.42) 𝑇𝐷𝑜𝑤𝑛=−1 𝜇ln𝑈2, (2.43) sendo 𝑇𝑈𝑝/𝐷𝑜𝑤𝑛 o tempo residência em cada estado de funcionamento; 𝑈1 e 𝑈2 são amostras uniformemente distribuídas entre um intervalo {0,1}. Up Down λ μ 30 Estado da Arte B. Figura 2.9 – Modelo de Markov de multiestados [45] A Figura 2.9 exemplifica um modelo de Markov de multiestados, que pode ser interpretado como a agregação do ciclo falha/restauro de N componentes. Nesta situação, a capacidade máxima de cada estado, 𝐶𝑘 é obtida de acordo com 𝐶𝑘=(𝑁−𝑘)×𝐶, (2.44) sendo 𝑘=0, 1, …, N e C a capacidade de cada componente. Se a duração em cada estado é exponencialmente distribuída, o tempo de residência em cada estado 𝐶0 e 𝐶𝑁 é 𝑇𝐶0=− 1 𝑁𝜆ln𝑈1, (2.45) 𝑇𝐶𝑁=− 1 (𝑘+1)𝜇ln𝑈2, (2.46) onde 𝑇𝐶0 é o tempo de residência no estado 𝐶0; 𝑇𝐶𝑁 é o tempo de residência no estado 𝐶𝑁; 𝜆 e 𝜇 são respetivamente taxas de avaria e de reparação de cada componente agregado. O tempo de residência em cada estado é calculo por [7] 𝑇𝐶𝑘=𝑚𝑖𝑛{−1 (𝑁−𝑘)𝜆ln𝑈1,−1 𝑘𝜇ln𝑈2}. (2.47) Normalmente são definidos parâmetro dos componentes MTTF (Mean Time To Fail) e MTTR (Mean Time To Repair). Estes podem ser obtidos por 𝑀𝑇𝑇𝐹=1 𝜆, (2.48) 𝑀𝑇𝑇𝑅=1 𝜇. (2.49). Por outro lado, o segundo modelo enunciado preocupa-se com a representação da dependência temporal da capacidade de cada componente. Sendo isto de extrema importância devido à capacidade de alguns componentes seguirem uma variação horária, mensal e até mesmo anual C0C1 μ CkCN (N-k)λ (k+1)μ Nλ ... (N-1)λ(N-k+1)λ kμ2μ Métodos de Simulação 31 [45], [47], [48]. Esta caraterística de dependência temporal é dada por séries horárias e/ou mensais com probabilidades associadas. Estas séries são obtidas através de anos de observações. 2.2.3.1 Unidades Térmicas As unidades térmicas são as que convertem energia térmica contida nos combustíveis fosseis ou nos núcleos atómicos por via de um ciclo térmico-dinâmico. O modelo de Markov de dois estados é usado para representar o ciclo de falha/restauro deste tipo de unidade de geração. As transições entre estados seguem uma distribuição probabilística exponencial. Quando a central é capaz de produzir à sua capacidade máxima, então está no estado de funcionamento (up state). Caso contrário está no estado de não funcionamento (down state). A capacidade de produção deste tipo de unidade são pré-especificadas e fixas a longo-prazo. 2.2.3.2 Unidades Hídricas Estas unidades de produção convertem a energia potencial da água em eletricidade. E tal como as unidades térmicas, o ciclo de falha/restauro é representado por um modelo de Markov de dois estados, em que as transições entre estados seguem uma distribuição probabilística exponencial. O modelo da dependência temporal da capacidade da unidade hídrica é complexo devido á dependência nos níveis de água disponíveis, tanto em albufeira como em fio de água. Com uma simples e robusta aproximação, é comum usar-se algumas séries temporais hidrológicas. Baseadas em eventos históricos, estabelece-se uma relação proporcional entre o total de água armazenada e a produção de energia através das unidades hídricas, em cada mês do ano [45], [47], [48]. O uso deste modelo simples é justificado devido à elevada complexidade em estabelecer condições de despacho a longo-prazo [48]. A capacidade disponível em cada mês é calculada através da multiplicação da capacidade máxima, que é obtida pelo modelo estocástico falha/restauro, com o correspondente valor mensal da série temporal hidrológica. As séries temporais hidrológicas são construídas baseadas em registos de vários anos. Apesar da natureza probabilística das séries hidrológicas, estas são um parâmetro determinístico do método de Simulação Sequencial de Monte Carlo. [47]. 2.2.3.3 Parques Eólicos Os aerogeradores (WTG – Wind Turbine Generator) são os responsáveis por converter energia do vento em eletricidade. Tradicionalmente, os aerogeradores de um parque eólico são 32 Estado da Arte iguais. Para tal, é usado um modelo de Markov multiestado para representar o ciclo de falha/restauro de um parque eólico. Neste tipo de modelo, as transições entre estados seguem uma distribuição probabilística exponencial. Tal como o modelo de dependência temporal da capacidade hídrica, a capacidade máxima dos estados dos aerogeradores é multiplicada pelo correspondente valor horário da série temporal eólica. Estas séries são compostas por valores horários das produções dos parques eólicos, em percentagem da capacidade total. Apesar de haver uma natureza probabilística, e tal como acontece nas unidades hídricas, também nos parques eólicos são usados os mesmos parâmetros no método de Simulação Sequencial de Monte Carlo. A vantagem do uso deste método é a ausência de um modelo para correlações espácio-temporais para parques eólicos de diferentes pontos geográficos, desde que, estas correlações sejam implicitamente incluídas nas séries e, naturalmente, contabilizadas para a simulação sequencial através do tempo. As incertezas da previsão eólica a curto-prazo são também incluídas no processo de simulação. 2.2.3.4 Carga A carga é modelada através do uso de uma representação sequencial que contém o nível da carga para cada hora do ano. O modelo de Simulação Sequencial de Monte Carlo segue cronologicamente estes passos da carga com o progresso da simulação. Cada barramento que suporta a carga tem o seu próprio perfil de carga, em valores percentuais em relação ao respetivo valor de pico. Cada percentagem é obtida dividindo o pico de carga da respetiva hora pelo valor de pico anual. Devido à flexibilidade do método de Simulação Sequencial de Monte Carlo, não há necessidade do uso discreto de representações de distribuições gaussianas. A inclusão de incertezas a curto-prazo nas simulações sequenciais consistem em adicionar ao valor horário de carga, um valor aleatório da distribuição gaussiana contínua de média zero e desvio standard, 𝜎𝐶𝑃 (Figura 2.10). De igual forma, às incertezas a longo-prazo são contabilizadas através da adição ao valor de pico anual da carga, um valor aleatório de uma distribuição gaussiana contínua de média zero e desvio standard, 𝜎𝐿𝑃 (Figura 2.10). Estes desvios controlam a amplitude da incerteza e são pré-especificados no início das simulações [45]. Métodos de Simulação 33 Figura 2.10 – Função de densidade de probabilidade da distribuição normal [8] A incerteza a curto-prazo da carga horária pode ser obtida através do uso da seguinte expressão [8] 𝑃(𝑡)=𝑃(𝑡)×(1+𝜎𝐶𝑃×𝑆𝑇𝑈), (2.50) onde 𝑃(𝑡) é o valor da carga para a hora t; 𝜎𝐶𝑃 é o desvio standard distribuição gaussiana contínua de média zero para as incertezas a curto-prazo; e STU é a incerteza de curto-prazo, em valor percentual. A incerteza a longo-prazo é obtida de forma semelhante, mas neste caso, o valor de pico de carga é alterado para o valor de pico de carga anual 𝑃𝑀á𝑥=𝑃𝑀á𝑥×(1+𝜎𝐿𝑃×𝐿𝑇𝑈), (2.51) onde 𝑃𝑀á𝑥 é o valor horário da carga; 𝜎𝐿𝑃 é o desvio standard distribuição gaussiana contínua de média zero para as incertezas a longo-prazo; e LTU é a incerteza de longo-prazo, em valor percentual [8]. 40 Estado da Arte Figura 2.15 – Resultado de uma simulação parcial [8] A Figura 2.15 mostra um exemplo do resultado de uma simulação parcial, em que se pode ver os índices de fiabilidade para cada tipo de reserva: LOLE, EENS, LOLF. Os restantes índices são de bem-estar (well-being) e não serão considerados neste trabalho: FREQ(S), FREQ(M), PROB(M). Este programa oferece uma panóplia de análise de resultados que se poderia introduzir aqui. Todo o seu conteúdo informativo está disposto em [8]. Previsão Eólica: Caso de Estudo 3.1 Introdução Neste Capítulo 3 é analisada uma série de 3 anos da produção de um parque eólico, adaptado de um concurso de previsão eólica [50]. Esta série eólica é a agregação de três anos de sete parques eólicos. Com estes dados, parametrizou-se um modelo ARIMA que melhor se ajusta à série, de forma a obter-se melhores resultados de previsão do que os obtidos pelo modelo Naïve. Tabela 3.1 – Tabela de valores estatísticos dos dados Ano 1 Ano 2 Ano 3 Total Média 0.2632 0.3048 0.3017 0.2899 Mediana 0.1973 0.2277 0.2201 0.2151 σ 0.2175 0.2467 0.2467 0.2381 Mínimo 0.0000 0.0000 0.0000 0.0000 Máximo 0.9354 0.9519 0.9349 0.9519 Como todo o conjunto de valores torna-se de pesada visualização, pela Tabela 3.1 conseguese ter uma ideia estatística de como variam os valores, em que todos os valores estão em relativos à capacidade máxima instalada. 42 Previsão Eólica: Caso de Estudo 3.2 Modelo de Previsão ARIMA 3.2.1 Estudo da Série Cada ano representa a medição horária da potência gerada pelo parque eólico, ou seja, 8760 elementos de cada ano. O objetivo da análise destes dados é a modelização de um ARIMA(p,d,q), para futuras previsões e comparar com os resultados de previsão do modelo Naïve. Figura 3.1 – Potência gerada pelo parque eólico ao longo de 3 anos Figura 3.2 – A. Autocorrelação da série; B. Correlação parcial da série A Figura 3.1 mostra o gráfico da agregação da produção dos sete parques ao longo dos três anos, verificando que a potência produzida nunca atinge a potência nominal. Existe uma periocidade na produção, isto é, existe alturas do ano em que a produção eólica é maior comparativamente a outras épocas do ano. Facilmente se observam três períodos de maior produção, intercalados com “vales”. Poder-se-á dizer que existe uma sazonalidade, no entanto, a Modelo de Previsão ARIMA 43 modelização do ARIMA é para muito curto-prazo, e não é de todo interessante saber os padrões mensais ou anuais, mas sim saber o que se passou nas horas anteriores mais próximas. Na Figura 3.2.A está representado o gráfico da autocorrelação da série, observa-se uma descida lenta dos valores para zero, mostrando de que a série é não estacionária. A linha vermelha tracejada representa o valor crítico da autocorrelação da série, dada pela equação (2.7), em que neste caso 𝑇=3×8760 ℎ. Já na Figura 3.2.B, gráfico das correlações parciais, o caso muda de figura e percebe-se claramente que os atrasos 1, 2, 3 e 4 estão acima do valor crítico, concluindo que estes atrasos influenciam significativamente o valor yt, sobretudo os atrasos 1 e 2. Pela correlação parcial percebe-se qual será, possivelmente, a ordem de regressão p. Apesar de ser evidente que a série é não estacionária, é necessário efetuar o teste KPSS para o confirmar. Este teste é feito até 20 atrasos, sob o critério do p-value ser maior que 0,05, rejeitando a hipótese nula da série ser não estacionária. Tabela 3.2 – Resultados do teste KPSS para a série eólica Atraso Resultado p-value 0 Verdadeiro 0.01 1 Verdadeiro 0.01 2 Verdadeiro 0.01 3 Verdadeiro 0.01 4 Verdadeiro 0.01 5 Verdadeiro 0.01 6 Verdadeiro 0.01 7 Verdadeiro 0.01 8 Verdadeiro 0.01 9 Verdadeiro 0.01 10 Verdadeiro 0.01 11 Verdadeiro 0.01 12 Verdadeiro 0.01 13 Verdadeiro 0.01 14 Verdadeiro 0.01 15 Verdadeiro 0.01 16 Verdadeiro 0.01 17 Verdadeiro 0.01 18 Verdadeiro 0.01 19 Verdadeiro 0.01 20 Verdadeiro 0.01 A Tabela 3.2 dos resultados do teste KPSS mostra que a hipótese de que a série ser estacionária é rejeitada. Desta forma, a série necessita de ser diferenciada, atribuindo d=1. 44 Previsão Eólica: Caso de Estudo Figura 3.3 – Série eólica diferenciada Figura 3.4 – A. Autocorrelação da série diferenciada; B. Correlação parcial da série diferenciada A Figura 3.3 mostra o gráfico da série eólica diferenciada, o que na prática representam as variações da produção de uma hora para a outra. Existe uma concentração dos valores em torno de zero, mas pode-se ver picos de variações. Na Figura 3.4, os primeiros valores dos gráficos foram retirados devido a serem mais elevados que os restantes (≈1), sendo possível assim, uma melhor interpretação dos valores correlacionais, bem como uma melhor visualização da banda de limite crítico. A Figura 3.4.A, onde é representado o gráfico das autocorrelações, mostra que os valores têm um comportamento diferente dos valores da Figura 3.2.A, pois aproximam-se bastante de zero, o que leva a crer que esta série diferenciada pode ser estacionária. Pode-se ver que os atrasos de 2 a 4 estão fora da banda do limite crítico, especialmente o segundo e o terceiro que têm valores mais elevados. Também outros valores estão fora da banda, no entanto, não são significativos em Modelo de Previsão ARIMA 45 comparação com os três primeiros. Na Figura 3.4.B, também estes atrasos estão fora da banda do limite crítico, estando os restantes em torno da banda do limite crítico. Com estes resultados estima-se que a ordem p de autorregressão não será maior que 4, pois a partir deste atraso, os valores já quase ou nada se correlacionam com o valor em yt. Para confirmar se esta série diferenciada também é estacionária, realizou-se novamente um teste KPSS. Tabela 3.3 – Resultados do teste KPSS para a série eólica diferenciada Atraso h p-value 0 Falso 0.1 1 Falso 0.1 2 Falso 0.1 3 Falso 0.1 4 Falso 0.1 5 Falso 0.1 6 Falso 0.1 7 Falso 0.1 8 Falso 0.1 9 Falso 0.1 10 Falso 0.1 11 Falso 0.1 12 Falso 0.1 13 Falso 0.1 14 Falso 0.1 15 Falso 0.1 16 Falso 0.1 17 Falso 0.1 18 Falso 0.1 19 Falso 0.1 20 Falso 0.1 A Tabela 3.1 mostra pelo teste KPSS que a série diferenciada é estacionária. Caso contrário, dever-se-ia proceder novamente à diferenciação da série. 3.2.2 Parametrização do Modelo ARIMA Depois de efetuado o teste estacionário e tendo sido diferenciada, introduz-se a série num ciclo para que seja possivel determinar todas as combinações possiveis das ordens p e q entre 1 e o último atraso fora da banda do limite crítico, que para o caso seria 4. Para cada combinação, é calculado o seu respetivo AIC e BIC e a seleção do modelo baseia-se na minimização do critério BIC. 46 Previsão Eólica: Caso de Estudo Figura 3.5 – Resultado do ciclo com o critério de minimização do BIC A Figura 3.5 mostra um ARIMA(4,1,3), resultado do critério de seleção pela minimização do BIC. Como seria de esperar, a ordem de autorregressão afirmou-se ser efetivamente 4. Tabela 3.4 – Critérios de ordem e estimação dos parâmetros ARIMA(4,1,3) AIC BIC -96451 -96430 A Tabela 3.4 mostra qual o resultado do AIC e BIC do modelo ARIMA(4,1,3). Segundo o critério de seleção, este é o modelo que melhor se ajusta aos dados da série. Para chegar a este resultado e estimar os valores dos parâmetros do modelo, foi utilizada a função automática estimate, do programa MatLab [51]. Esta função utiliza como critério o MLE. Tabela 3.5 – Valores dos parâmetros para o modelo ARIMA(4,1,3) d 𝛟𝟏 𝛟𝟐 𝛟𝟑 𝛟𝟒 𝛉𝟏 𝛉𝟐 𝛉𝟑 σ2 c 1 1.8367 -1.0622 0.2696 -0.0601 -1.3998 0.3522 0.0486 0.0015 0.0000 Matematicamente, o modelo é escrito do seguinte modo 3.2.3 Diagnóstico dos Resíduos As previsões foram realizadas através da função forecast do MatLab [52], que inserindo o modelo ARIMA correspondente, determinava as previsões para a hora t+1, para um intervalo de confiança de 95%. Através da expressão dos erros de previsão facilmente se determinavam os resíduos 𝑦𝑡=𝑦𝑡−1+1.8367∙𝑦𝑡−1−1.0622∙𝑦𝑡−2+0.2696∙𝑦𝑡−3−0.0601∙𝑦𝑡−4+𝑒𝑡 −1.3998∙𝑒𝑡−1+0.3522∙𝑒𝑡−2+0.0486∙𝑒𝑡−3 (3.1) Modelo de Previsão ARIMA 47 sendo 𝑦𝑡 o valor observado e 𝑦𝑡 o valor previsto. No entanto, através função infer do MatLab [53] e fornecendo o modelo ARIMA(p,d,q) bem como a série, determinaram-se os resíduos do respetivo modelo. Com os resíduos avalia-se o modelo usando MAE, RMSE e MAPE. Tabela 3.6 – Avaliação dos resíduos para o ARIMA(4,1,3) ME MAE RMSE SMAPE 0.00% 2.73% 3.86% 7.92% A Tabela 3.6 mostra os diferentes tipos de avaliação dos resíduos. Posteriormente neste documento, estas avaliações serão comparadas com as do Naïve em 3.4. Um dos requisitos destes resíduos seria a sua média (ME) ser zero, o que se confirma na Tabela 3.6. Outra passava pela análise da sua autocorrelação (Figura 3.7). Figura 3.6 – Resíduos do modelo de previsão ARIMA(4,1,3) Figura 3.7 – A. Autocorrelação dos resíduos; B. Correlação parcial dos resíduos 𝑒𝑡=𝑦𝑡−𝑦𝑡 , (3.2) 48 Previsão Eólica: Caso de Estudo Para melhor visualização, foi retirado o primeiro valor de correlação, tal como nos gráficos da Figura 3.4, devido a ser muito maior que os restantes. Por isso, facilmente se visualiza nas Figura 3.7.A e Figura 3.7.B de que os resíduos são não correlacionados. Em ambos, vê-se que todos os valores se encontram dentro dos limites críticos, com a exceção dos atrasos 18, 23 e 24 que saem ligeiramente dos limites. Contudo, esta banda é muito estreita, {-0.0123, 0.0123}, e de valores muito baixos em comparação ao primeiro valor de autocorrelação (≈ 1), desprezando-se assim estes atrasos que saem dos limites. Tabela 3.7 – Teste Ljung-Box para o modelo ARIMA(4,1,3) h h-K χ2 α,h-K Q* p-value Resultado 10 3 7.81 14.07 0.0028 Verdadeiro 11 4 9.49 17.67 0.0014 Verdadeiro 12 5 11.07 17.76 0.0033 Verdadeiro 13 6 12.59 18.29 0.0055 Verdadeiro 14 7 14.07 19.14 0.0078 Verdadeiro 15 8 15.51 19.70 0.0115 Verdadeiro 16 9 16.92 19.71 0.0198 Verdadeiro 17 10 18.31 22.72 0.0118 Verdadeiro 18 11 19.68 28.63 0.0026 Verdadeiro 19 12 21.03 31.20 0.0018 Verdadeiro 20 13 22.36 32.73 0.0019 Verdadeiro Aplicando o teste Ljung-Box no MatLab, sendo nível de confiança 95%, o resultado “Verdadeiro” indica rejeição da hipótese nula da não autocorrelação residual, para o respetivo atraso h e grau de liberdade h-K. Verifica-se que, em todos eles, a hipótese é rejeitada. Entra-se aqui numa discordância: pela média zero, visualização dos gráficos de autocorrelação e correlação parcial, estes resíduos identificam-se como ruído branco. No entanto, o teste Ljung-Box afirma que a hipótese de não autocorrelação residual é rejeitada. Esta contradição torna os resultados confusos e suspeitos. Existem algumas suposições que podem justificar estes resultados:  O teste pode não estar a ser bem aplicado e/ou não adequado para o modelo ARIMA, uma vez que nas referências [9], [11] e [34], apenas é apontado como recomendado para modelos ARMA;  O conjunto de dados pode não ser adequado para ser aplicado a este tipo de teste. Devido ao tamanho elevado da série, o teste pode perder força com esta dimensão, [11]. Modelo de Previsão Naïve 49 3.3 Modelo de Previsão Naïve Sendo este o modelo de previsão aplicado no programa RESERVAS, é o mais básico de todos os modelos de previsão. Como não há qualquer métrica de cálculo e/ou determinação de qualquer tipo de ajustamento do modelo à série onde é aplicado, rapidamente se obtém resultados. Figura 3.8 – Resíduos do modelo de previsão Naïve Figura 3.9 – A. Autocorrelação dos resíduos; B. Correlação parcial dos resíduos Quando aplicando este modelo à série em estudo, verifica-se que não existe diferença entre os resíduos do modelo de previsão Naïve (Figura 3.8) e a série diferenciada (Figura 3.3), tal como seria de esperar. Tabela 3.8 – Avaliação dos resíduos para o Naïve ME MAE RMSE SMAPE 0.00% 3.08% 4.34% 8.72% 56 Avaliação da Reserva Operacional 4.2 IEEE RTS 79 HW O sistema usado nas simulações foi o IEEE RTS 79 HW, como referido anteriormente, é uma versão alterada do original IEEE RTS 79 [54] (Anexos, Tabela A.1) que inclui unidades hídricas e eólicas no sistema produtor. O perfil da carga varia entre uma potência mínima de 965.25 MW, e 2850 MW de pico de carga. A capacidade de produção hídrica é a média mensal de uma série de cinco anos da referência [48]. Figura 4.1 – Carga anual do IEEE RTS 79 HW A Figura 4.1 mostra o perfil da carga do sistema, sendo que o valor de pico é atingido no final da série, ou mais precisamente, no final da série os valores são os mais elevados. Este é um ponto bastante importante a ter em conta, pois é nesta altura do ano que a carga exige mais do sistema produtor. 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 01000 2000 3000 4000 5000 6000 7000 8000 9000 Carga (% Carga Máx.) Tempo (h) IEEE RTS 79 HW 57 Figura 4.2 – Produção hídrica mensal vs Carga média mensal O sistema produtor tem seis geradores hídricos com 50 MW cada um. Na Figura 4.2, está a comparação da capacidade de produção hídrica mensal (azul), em percentagem da capacidade máxima instalada hídrica, e da carga mensal (laranja), estando esta em relação ao valor de pico de carga. Apesar dos valores estarem em valores percentuais dos respetivos valores máximos (capacidade máxima de produção e respetivo pico de carga) existem alturas do ano em que a carga é maior que a capacidade de produção hídrica, em termos de valores relativos. Tanto para as primeiras 744 horas, bem como para os últimos dois meses, onde a carga atinge valores de pico, a produção hídrica atinge os valores mais baixos da sua produção. As unidades térmicas representam a grande maioria da capacidade de produção do sistema, com uma capacidade de geração de 3105MW, entre as 26 unidades térmicas. Quando se agrega a potência eólica, todas as unidades eólicas agregam-se num único parque eólico, aumentando o número de geradores do sistema de 32 para 33. Como o trabalho incide na busca do impacto do erro de previsão eólica, a capacidade de geração eólica varia de 800, 1600 e 2400 MW de capacidade instalada. É também avaliada a reserva operacional sem produção eólica. No planeamento a curto-prazo, neste sistema teste didático, definiram-se duas ordens de mérito, apresentadas na Tabela 4.1, de forma a simular diferentes casos. 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 01000 2000 3000 4000 5000 6000 7000 8000 9000 Potência (% Cap. Máx) Tempo (h) Hídrica Carga 58 Avaliação da Reserva Operacional Tabela 4.1 – Ordem de mérito das unidades produtoras Prioridade de Escalonamento Ordem de Mérito A B 1ª Eólica Eólica 2ª Térmica Hídrica 3ª Hídrica Térmica Tal como dito anteriormente, pretende-se fazer o estudo da reserva operacional quando esta está a subir, pretendendo-se com isso, evitar que a reserva operacional entre em nível excessivo por uma má previsão. Neste caso, só é interessante retirar os valores negativos dos erros de previsão, isto quer dizer que 𝛥𝑃𝑊𝑡<0⇔𝑃𝑊𝑡−𝑃𝑊𝑡<0⇔𝑃𝑊𝑡<𝑃𝑊𝑡, (4.1) sendo 𝑃𝑊𝑡 o valor observado e 𝑃𝑊𝑡 o valor previsto. Por esta razão, a equação de risco da reserva operacional (2.52) terá de ser reescrita da seguinte forma 𝑅𝑂𝑃𝐸<𝛥𝐿−𝛥𝑃𝑊+𝛥𝐺. (4.2) Figura 4.3 – Parte negativa dos erros de previsão Na Figura 4.3 estão representados os erros médios horários da previsão de três anos, obtidos pelos dois modelos de previsão eólica. Apesar do modelo ARIMA ser melhor do que o Naïve em termos de valores médios, como visto na secção ARIMA vs Naïve, existem alturas da série onde isso não se verifica. Entre a hora 7296 e o final da série, os erros de previsão são superiores, em valores absolutos. -0,14 -0,12 -0,1 -0,08 -0,06 -0,04 -0,02 001000 2000 3000 4000 5000 6000 7000 8000 9000 Erros previsão (% Cap. Instalada) Tempo (h) Naïve ARIMA Resultados das Simulações 59 4.3 Resultados das Simulações As simulações realizadas diferem nas ordens de mérito dos geradores, produção eólica e respetivo modelo de previsão. A Tabela 4.2 mostra o sumário dos testes realizados. Tabela 4.2 – Sumário dos testes Número aerogeradores Cap. unitária eólica Cap. instalada eólica Cap. Máx. Produtora MW MW MW 0 2 0 3405 400 800 4205 800 1600 5005 1200 2400 5805 Pico Carga Valley Load Número Unidades Térmicas Número Unidades Hídricas Número Parques Eólicos MW MW 2850 965.62 26 6 1 Com diferentes quantidades de aerogeradores, existem diferentes capacidades máximas do sistema produtor. Considera-se que cada gerador eólico tem capacidade de 2MW, como pode ser visualizado nos Anexos, Tabela A.1. O perfil da carga pode ser visualizado na Figura 4.1. Com pico de carga de 2850 MW e o mínimo de 965,62 MW, verifica-se que a capacidade máxima produtora sem produção eólica é superior ao pico de carga (3405MW sem produção eólica). Um reforço de capacidade eólica no sistema aumenta esta diferença entre o pico de carga e a capacidade máxima instalada de geração. Foi admitido que todos os aerogeradores foram aglomerados no mesmo parque eólico. Com isto, admite-se que todos os aerogeradores são afetados pelo mesmo comportamento de vento, ou seja, as previsões de produção eólica em cada gerador seguirão a mesma série, bem como os respetivos valores observados. Em todas as simulações, foram mantidos os seguintes parâmetros Tabela 4.3 – Parâmetros das simulações Parâmetro Convergência β 2.50% Incerteza de carga a curto-prazo STU 2.50% Reserva Primária RP 50 MW Reserva Secundária RS 200 MW 60 Avaliação da Reserva Operacional 4.3.1 Caso sem Produção Eólica As simulações do sistema sem unidades eólicas incorporadas no sistema resultam numa avaliação da reserva operacional quando não há reforço de geração por via eólica. Tabela 4.4 – Índices de fiabilidade sem produção eólica LOLP LOLE EPNS EENS LOLF LOLD prob. h/ano MW MWh/ano oc./ano h/oc. Ordem de mérito A 0.00056 4.90 0.0299 262.04 4.89 1.002 Ordem de mérito B 0.00277 24.28 0.2372 2078.17 24.23 1.002 Claramente que quando as unidades hídricas são as últimas a entrar em funcionamento, há uma maior flexibilidade da reserva terciária em satisfazer as variações da carga e de geração. Já no caso da ordem de mérito B, em que as unidades hídricas entram antes das unidades térmicas, os índices de Reserva Operacional pioram, devido à pouca flexibilidade de produção das unidades térmicas de fazer face a estas variações. Para algumas destas unidades térmicas, uma variação na produção, ou entrada em funcionamento, leva mais do que uma hora, não sendo as ideais para atuarem a curto-prazo. Claramente que as unidades hídricas são bem mais flexíveis e melhores para serem usadas na reserva terciária. 4.3.2 Caso com Produção Eólica Nesta subsecção introduz-se produção eólica no sistema, aumentando a capacidade máxima de geração do sistema produtor. De forma a simular vários cenários de eólica, foram considerados três valores de potência: 800, 1600 e 2400 MW. E foram testados os dois modelos de previsão, ARIMA e Naïve. 4.3.2.1 ARIMA Tal como se observou no Capítulo 3, este é o melhor modelo de previsão comparativamente ao Naïve, que quando aplicado no programa RESERVAS, obtiveram-se os seguintes resultados da avaliação da reserva operacional. Resultados das Simulações 61 Tabela 4.5 – Avaliação da ROPE para o modelo de previsão ARIMA Nº geradores eólicos Cap. unitária eólica Cap. instalada eólica LOLP LOLE EPNS EENS LOLF LOLD MW MW prob. h/ano MW MWh/ano oc./ano h/oc. Ordem de Mérito A 400 2 800 0.00043 3.76 0.0240 210.11 3.76 1.001 800 1600 0.00055 4.81 0.0349 305.98 4.80 1.002 1200 2400 0.00186 16.28 0.1443 1263.92 16.20 1.005 Ordem de Mérito B 400 2 800 0.00305 26.75 0.2119 1856.13 26.65 1.004 800 1600 0.00936 82.00 0.4819 4221.03 80.00 1.025 1200 2400 0.02409 211.00 1.6141 14139.23 202.50 1.042 À semelhança da Tabela 4.4, também neste caso se pode concluir que a ordem de mérito influencia bastante os índices de fiabilidade da reserva operacional. A Tabela 4.5 mostra que quando maior é a potência eólica instalada, piores são os seus índices. Devido aos erros de previsão eólica, quanto maior fora a potência instalada no parque, maior irá ser o 𝛥𝑃𝑊 em MW. Para um 𝛥𝑃𝑊 (% 𝐶𝑎𝑝.𝐼𝑛𝑠𝑡𝑎𝑙𝑎𝑑𝑎)=0,01, e um 𝑃𝑊𝑀á𝑥(𝑀𝑊)=1600 𝑀𝑊 (para quando não se consideram as saídas por falhas furtuitas dos aerogeradores), o 𝛥𝑃𝑊(𝑀𝑊)=16 𝑀𝑊. Mas se a potência eólica máxima aumentar para 2400 MW, então o erro de previsão leva a um erro de 24 MW. Isto leva a uma maior probabilidade de insucesso de reserva, devido a limitações das unidades hídricas em termos de potência disponível, e de unidades térmicas, em termos de resposta às variações a curto-prazo. Com o aumento da potência eólica, os erros atingem magnitudes que a potência hídrica tem dificuldade de compensar. 4.3.2.2 Naïve Usando o modelo Naïve para a previsão da produção eólica, realizou-se a avaliação dos índices de fiabilidade da reserva operacional com os seguintes resultados. Tabela 4.6 – Avaliação da ROPE para o modelo de previsão Naïve Nº geradores eólicos Cap. unitária eólica Cap. instalada eólica LOLP LOLE EPNS EENS LOLF LOLD MW MW prob. h/ano MW MWh/ano oc./ano h/oc. Ordem de Mérito A 400 2 800 0.00043 3.74 0.0241 211.04 3.74 1.001 800 1600 0.00053 4.67 0.0300 262.91 4.66 1.002 1200 2400 0.00263 23.04 0.1853 1623.15 21.45 1.074 Ordem de Mérito B 400 2 800 0.00325 28.46 0.2166 1897.58 28.24 1.008 800 1600 0.01373 120.28 0.7112 6230.37 107.22 1.122 1200 2400 0.03288 288.00 2.2237 19479.35 238.50 1.208 62 Avaliação da Reserva Operacional Tal como aconteceu com o modelo de previsão ARIMA, o efeito da ordem de mérito é igual: quando a produção hídrica está fora da base, melhores são os seus índices. Neste caso, as unidades hídricas fazem a parte da reserva terciária, atendendo assim, às eventuais variação de carga, erro de previsão eólica e variação na produção. A Tabela 4.6 mostra esses resultados. Quanto maior a potência instalada eólica, maior é o índice LOLP, bem como o índice EENS, pelas mesmas razões dadas no modelo ARIMA. Discussão de Resultados 63 4.4 Discussão de Resultados Nas subsecções anteriores analisaram-se os resultados obtidos da avaliação do sistema produtor do IEEE RTS 79 HW, tendo este capacidade máxima produtora, sem produção eólica, de 3405 MW. Apesar de haver capacidade produtora para a carga, devido ao erro de previsão eólica e variação da carga, há um risco associado de não existir flexibilidade para atender às variações de curto-prazo. Tabela 4.7 – Comparação dos índices de fiabilidade para a ordem de mérito A Cap. instalada eólica Modelo previsão LOLP LOLE EPNS EENS LOLF LOLD MW prob. h/ano MW MWh/ano oc./ano h/oc. Ordem de Mérito A - - 0.00056 4.90 0.0299 262.04 4.89 1.002 800 ARIMA 0.00043 3.76 0.0240 210.11 3.76 1.001 Naïve 0.00043 3.74 0.0241 211.04 3.74 1.001 1600 ARIMA 0.00055 4.81 0.0349 305.98 4.80 1.002 Naïve 0.00053 4.67 0.0300 262.91 4.66 1.002 2400 ARIMA 0.00186 16.28 0.1443 1263.92 16.20 1.005 Naïve 0.00263 23.04 0.1853 1623.15 21.45 1.074 A Tabela 4.7 mostra a comparação da avaliação da reserva operacional, quando utilizados os dois modelos de previsão, bem como, quando não se consideram unidades eólicas na geração. Como visto no Capítulo 3 e dito na secção 4.2, o modelo ARIMA é um melhor modelo de previsão, em termos de valores médios, do que o Naïve. Por isso, esperava-se que a avaliação da reserva fosse melhor do que o Naïve. Analisando as simulações de 800 e 1600 MW, verifica-se de que os índices são bastante semelhantes, tendo o Naïve melhores índices. Mas isso é diferente para a simulação de 2400 MW. 64 Avaliação da Reserva Operacional Figura 4.4 – Valores do LOLP horários para potência eólica de 1600 MW Como visualizado na Figura 4.4, os erros maiores do modelo de previsão ARIMA acontecem entre {7296, 8760}. Nesse mesmo intervalo de tempo, a carga atinge o seu pico bem como os valores das cargas médias mensais dos últimos dois meses são as maiores. E nesses mesmos últimos meses, a capacidade de geração por fontes hídricas é a menor (Figura 4.2). Visualizando a Figura 4.4 verifica-se que há maior perda de carga nos instantes referidos. Uma vez que os erros do ARIMA são maiores do que os do Naïve, em valores absolutos, entre 7296 e 8760h há uma afirmação maior da perda de carga deste modelo. Para este caso de simulação, os erros do Naïve não estão tão afirmados ao longo da série, fazendo com que o ARIMA se torne um mau modelo apenas porque estar a errar mais nos últimos instantes. Figura 4.5 – Reserva Operacional solicitada na simulação A-1600 MW A Figura 4.5 representa a reserva operacional solicitada para uma potência eólica de 1600 MW, e ordem de mérito A. Apesar de não ser muito percetível, estão representados graficamente 0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 LOLP (prob.) Tempo (h) Naïve Arima 280 305 330 355 380 405 430 455 01000 2000 3000 4000 5000 6000 7000 8000 9000 Potência (MW) Tempo (h) Reserva Operacional Arima Naïve Discussão de Resultados 65 os resultados das simulações de reserva operacional dos dois modelos, que se apresentam resultados muito semelhantes. Pelo desenho do gráfico, poder-se-á dizer que as unidades hídricas estão a assegurar a reserva terciária do sistema produtor. As variações de reserva são provocadas pela variação de carga. Estas são maiores nos últimos dois meses, meses onde a carga atinge os valores mais altos. Comparativamente à avaliação entre os sistemas com e sem produção eólica, chega-se à conclusão de que para uma capacidade instalada inferior a 1600 MW, há um menor LOLE, ou seja, o risco é menor. O facto de haver muita produção eólica, faz aumentar a incerteza das unidades de reserva operacional satisfazerem os erros de previsão. Isto quer dizer que a probabilidade da reserva operacional causar falha no sistema é maior. Figura 4.6 – Valores do LOLP horários para potência eólica de 2400 MW Quando a potência instalada eólica é de 2400 MW, a avaliação torna-se melhor usando o modelo ARIMA do que o Naïve. A Figura 4.6 mostra que os erros do Naïve provocam maior impacto na reserva operacional, isto devido à variação ∆𝑃𝑊 ser maior em valores de MW neste caso, do que no caso de 1600 MW. Aqui, o fator da ponta dos erros ARIMA não é tão significativo comparativamente aos erros do Naïve, e os erros deste modelo sobressaem-se mais ao longo da série. 0,00 0,05 0,10 0,15 0,20 0,25 0,30 0,35 0,40 01000 2000 3000 4000 5000 6000 7000 8000 9000 LOLP (prob.) Tempo (h) Naïve Arima 72 Conclusão Anexo 1. Glossário Falhas furtuitas (forced outage) [7]– Uma interrupção que resulta de condições de emergência diretamente associadas a um componente que deve ser retirado de serviço imediatamente, automaticamente ou com comutação de operações que podem ser executadas, ou uma interrupção causada por operação indevida de equipamentos ou erro humano. Likelihood [55] – em estatística, é uma função de parâmetros de um modelo estatístico que permite inferir sobre o seu valor a partir de um conjunto de observações. Este desempenha um papel fundamental na inferência estatística, especialmente em métodos de estimação de parâmetros. Reserva Estática [45] – os estudos de reserva estática pontam definir para um certo nível de risco, a capacidade necessária para encontrar a carga esperada. A incerteza associada a este tipo de reserva, que é uma variável estocástica, é causada pela intermitência dos recursos energéticos primários, as saídas planeadas e furtuitas, e a incerteza das cargas. Os eventos de reserva estática insuficiente acontecem quando 𝑅𝑆𝑇𝐴=𝐺−𝐿<0, (A.1) onde G representa a capacidade de produção e L a carga do sistema. 74 Glossário Anexo 2. IEEE RTS 79 HW Tabela A.1 – Informação geradores do sistema IEEE RTS 79 HW ID Número de geradores Nome Nome tecnologia Potência Mínima Potência Máxima MTTF MTTR Tipo Modelo Markov Estado Mobilizável MW MW (h) (h) 1 1 U20 Oil/CT 0 20 450 50 Dois estados 0 false 2 1 U20 Oil/CT 0 20 450 50 Dois estados 0 false 3 1 U76 Coal/Steam 0 76 1960 40 Dois estados 0 false 4 1 U76 Coal/Steam 0 76 1960 40 Dois estados 0 false 5 1 U20 Oil/CT 0 20 450 50 Dois estados 0 false 6 1 U20 Oil/CT 0 20 450 50 Dois estados 0 false 7 1 U76 Coal/Steam 0 76 1960 40 Dois estados 0 false 8 1 U76 Coal/Steam 0 76 1960 40 Dois estados 0 false 9 1 U100 Oil/Steam 0 100 1200 50 Dois estados 0 false 10 1 U100 Oil/Steam 0 100 1200 50 Dois estados 0 false 11 1 U100 Oil/Steam 0 100 1200 50 Dois estados 0 false 12 1 U197 Oil/Steam 0 197 950 50 Dois estados 0 false 13 1 U197 Oil/Steam 0 197 950 50 Dois estados 0 false 14 1 U197 Oil/Steam 0 197 950 50 Dois estados 0 false 15 1 U12 Oil/Steam 0 12 2940 60 Dois estados 0 false 16 1 U12 Oil/Steam 0 12 2940 60 Dois estados 0 false 17 1 U12 Oil/Steam 0 12 2940 60 Dois estados 0 false 18 1 U12 Oil/Steam 0 12 2940 60 Dois estados 0 false 19 1 U12 Oil/Steam 0 12 2940 60 Dois estados 0 false 20 1 U155 Coal/Steam 0 155 960 40 Dois estados 0 false 21 1 U155 Coal/Steam 0 155 960 40 Dois estados 0 false 22 1 U400 Nuclear 0 400 1100 150 Dois estados 0 false 23 1 U400 Nuclear 0 400 1100 150 Dois estados 0 false 24 1 U155 Coal/Steam 0 155 960 40 Dois estados 0 false 25 1 U155 Coal/Steam 0 155 960 40 Dois estados 0 false 26 1 U350 Coal/Steam 0 350 1150 100 Dois estados 0 false 27 1 U50 Hydro 0 50 1980 20 Dois estados 0 true 28 1 U50 Hydro 0 50 1980 20 Dois estados 0 true 29 1 U50 Hydro 0 50 1980 20 Dois estados 0 true 30 1 U50 Hydro 0 50 1980 20 Dois estados 0 true 31 1 U50 Hydro 0 50 1980 20 Dois estados 0 true 32 1 U50 Hydro 0 50 1980 20 Dois estados 0 true 33 100 W1 Wind 0 2 1914.74 80 Multiestado 0 false 76 IEEE RTS 79 HW Figura A.1 – Diagrama unifilar do sistema IEEE RTS 79 Referências Bibliográficas [1] Tim Smedley. "Goodbye nuclear power: Germany's renewable energy revolution," 18 August 2015; http://www.theguardian.com/sustainable-business/nuclear-powergermany-renewable-energy. [2] "Europe 2020 targets," 07 August 2015; http://ec.europa.eu/europe2020/europe-2020-ina-nutshell/targets/index_en.htm. [3] Kate Connolly. "G7 leaders agree to phase out fossil fuel use by end of century," 8 June 2015; http://www.theguardian.com/world/2015/jun/08/g7-leaders-agree-phase-outfossil-fuel-use-end-of-century. [4] Arthur Neslen. "Wind power generates 140% of Denmark's electricity demand," 07 Agosto 2015; http://www.theguardian.com/environment/2015/jul/10/denmark-windwindfarm-power-exceed-electricity-demand. [5] John Vidal. "UK and Germany break solar power records," 18 August 2015; http://www.theguardian.com/environment/2014/jun/23/uk-and-germany-break-solarpower-records. [6] REN. "Estatística Diária - SEN," 8 Agosto 2015 2015; http://www.centrodeinformacao.ren.pt/PT/InformacaoExploracao/Paginas/EstatisticaDi aria.aspx. [7] R. Billinton, and R. N. Allan, Realiability Evaluation of Power Systems, 2nd ed.: Plenum Press, 1994. [8] M. A.da Rosa, L. M. Carvalho, and R. Ferreira, RESERVE - Program Manual, 1.0 ed., Porto: INESC, 2011. [9] George E. P. Box, Gwilym M. Jenkins, and Gregory C. Reinsel, Time Series Analysis: Forecasting and Control 3rd ed., Englewood Cliffs, New Jersey, 1994. [10] Vladimiro Miranda, "Redes Neuronais – treino por retropropagação," Faculdade de Engenharia da Universidade do Porto, 2007, p. 10. 78 Referências Bibliográficas [11] R. J. Hyndman, and G. Athanasopoulos. "Forecasting: principles and practice," 25 June 2015. [12] Everette S. Gardner, “Exponential Smoothing: The State of the Art,” Journal of Forecasting, vol. 4, pp. 1-28, 1985. [13] R. J. Hyndman, R. A. Ahmed, G. Athanasopoulos, and H. L. Shang, "Optimal combination forecasts for hierarchical time series," 2010, p. 24. [14] R. Billinton, H. Chen, and R. Ghajar, "Time-series models for Reliability Evaluation of Power Systems Including Wind Energy," Microelectron. Reliab., 36, 1996. [15] R. Billinton, H. Chen, and R. Ghajar, "A Sequential Simulation Techinque for Adequacy Evaluation of Generating Systems Including Wind Energy," IEEE Transactions on Energy Conversion, 11, 1996. [16] R. G. Kavasseri, and K. Seetharaman, "Day-ahead wind speed forecasting using fARIMA models," 34, 2009. [17] J.C. Palomares-Salas, J.J. G. de la Rosa, J.G. Ramiro, J. Melgar, A. Agüera, and A. Moreno, "ARIMA vs. Neural Networks for Wind Speed Forecasting," 2009. [18] E. Cadenas, and W. Rivera, "Wind speed forecasting in the south coast of Oaxaca, Mexico," Renewable Energy 32, 2007. [19] B. G. Brown, R. W. Katz, and A. H. Murphy, “Time Series Models to Simulate and Forecast Wind Speed and Wind Power,” Journal of Climate and Applied Meteorology, vol. 23, pp. 1184-1195, 1984. [20] P. Chen, T. Pedersen, B. Bak-Jensen, and Z. Chen, "ARIMA-Based Time Series Model of Stochastic Wind Power Generation," IEEE Transactions on Power Systems, 25, 2010. [21] E. Erdem, and J. Shi, "ARMA based approaches for forecasting the tuple of wind speed and direction," Applied Energy, 88, 2011. [22] E. Erdem, J. Shi, and Y. Peng, "Short-Term Forecasting of Wind Speed and Power - A Clustering Approach," Industrial and Systems Engineering Research Conference, 2014. [23] G. N. Kariniotakis, G. S. Stavrakakis, and E. F. Nogaret, "Wind Power Forecasting Using Advanced Neural Networks Models," IEEE Transactions on Energy Conversion, 11, 1996. [24] G. Li, and J. Shi, "On comparing three artificial neural networks for wind speed forecasting," Applied Energy, 87, 2010. [25] T. G. Barbounis, J. B. Theocharis, M. C. Alexiadis, and P. S. Dokopoulos, "LongTermWind Speed and Power Forecasting Using Local Recurrent Neural Network Models," IEEE Transactions on Energy Conversion, 21, 2006. [26] M. Monfared, H. Rastegar, and H. Madadi Kojabadi, "A new strategy for wind speed forecasting using artificial intelligent methods," Renewable Energy, 34, 2009. Referências Bibliográficas 79 [27] J. M. Sloughther, T. Gneiting, and A. Raftery, “Probabilistic Wind Speed Forecasting Using Ensembles and Bayesian Model Averaging,” American Statistical Association, vol. 105, pp. 25-35, 2010. [28] R. Karki, and P. Hu, "Wind Power Simulation Model For Reliability Evaluation," 2005. [29] R. Karki, P. Hu, and R. Billinton, "A Simplified Wind Power Generation Model for Reliability Evaluation," IEEE Transactions on Energy Conversion, 21, 2006. [30] D. Huang, and R. Billinton, "Effects of Wind Power on Bulk System Adequacy Evaluation Using the Well-Being Analysis Framework," IEEE Transactions on Power Systems, 24, 2009. [31] M. A. Matos, and R. Bessa, “Operating Reserve Adequacy Evaluation using Uncertainties of Wind Power Forecast,” in IEEE Bucharest Power Tech Conference, Bucharest, Romania, 2009, pp. 1 - 8. [32] Rob J. Hyndman, Forecasting: Principles & practice: University of Western Australia, 2014. [33] G. E. P. Box, and D. A. Pierce, “Distribution of Residual Autocorrelations in Autoregressive-Integrated Moving Average Time Series Models,” Journal of the American Statistical Association, vol. 65, pp. 1509–1526, 1970. [34] G. M. Ljung, and G. E. P. Box, “On a Measure of a Lack of Fit in Time Series Models,” Biometrika, vol. 65, pp. 297-303, 1978. [35] Benito E. Flores, “A Pragmatic View of Accuracy Measurement in Forecasting,” Journal of Mgmt Sci., vol. 14, pp. 93-98, 1986. [36] Z. Chen, and Y. Yang, "Assessing Forecast Accuracy Measures," Department of Economics, Iowa State University, 2004, pp. 1-26. [37] Frank Wilcoxon, "Individual Comparisons by Ranking Methods," Biometrics Bulletin, 6, International Biometric Society, 1945, pp. 80-83. [38] H. B. Mann, and D. R. Whitney, “On a test of whether one of two random variables is stochastically larger than the other,” pp. 1-11, 1947. [39] V. S. Srinivasan, and K.V. Ramachandra, “Decision Making in Reliability,” IEEE Transactions on Reliability, vol. R-26, pp. 261-265, October 1977. [40] Roy Billinton, and Wenyuan Li, Reliability Assessment of Electrical Power Systems Using Monte Carlo Methods, New York: Plenum Press, 1994. [41] Vladimiro Miranda, Fiabilidade em Sistemas de Potência – Uma Introdução, 2.0 ed., Porto: INESC, 2006. [42] J. C. O. Mello, M. V. F. Pereira, and A. M. Leite da Silva, "Evaluation of reliability worth in composite systems based on pseudo-sequential monte carlo simulation," IEEE Transactions on Power Systems, 9, 1994. 80 Referências Bibliográficas [43] A. M. Leite da Silva, L. A. F. Manso, J. C. O. Mello, and R. Billinton, "PseudoChronological Simulation for Composite Reliability Analysis with Time Varying Loads," IEEE Transactions on Power Systems, 15, 2000. [44] A. M. Leite da Silva, R. A. González-Fernández, W. S. Sales, and L. A. F. Manso, "Reliability Assessment of Time-Dependent Systems via Quasi-Sequential Monte Carlo Simulation," 2010. [45] Leonel M. Carvalho, “Advances on the Sequential Monte Carlo Reliability Assessment of Generation-Transmission Systems using Cross-Entropy and Population-based Methods,” PhD, Universidade do Porto, Porto, 2013. [46] A. C. G. Melo, M. V. F. Pereira, and A. M. Leite da Silva, "Frequency and duration calculations in composite generation and transmission reliability evaluation," IEEE Transactions on Power Systems, 7, 1992. [47] M. A. Matos, J. P. Lopes, M. A. da Rosa, R. Ferreira, A. Leite da Silva, W. Sales, L. Resende, L.Manso, P. Cabral, M. Ferreira, N. Martins, C. Artaiz, F. Soto, and R. López, "Probabilistic evaluation of reserve requirements of generating systems with renewable power sources: The Portuguese and Spanish cases," Electrical Power and Energy Systems, 31, 2009. [48] A. M. Leite da Silva, W. S. Sales, L. A. F. Manso, and R. Billinton, “Long-Term Probabilistic Evaluation of Operating Reserve Requirements With Renewable Sources,” IEEE Transactions on Power Systems, vol. 25, no. 1, pp. 106-116 February 2010. [49] A. M. Leite da Silva, A. M. Cassula, R. Billinton, and L. A. F. Manso, "Integrated reliability evaluation of generation, transmission and distribution systems," IEE Proc.- Gener. Transm. Distrib, 149, 2002. [50] T. Hong, P. Pinson, and S. Fan, “Global Energy Forecasting Competition 2012,” International Journal of Forecasting, vol. 30, pp. 357-363, April - June 2014. [51] The MathWorks. "estimate," 29 July 2015; http://www.mathworks.com/help/econ/arima.estimate.html. [52] The MathWorks. "Forecast," 29 July 2015; http://www.mathworks.com/help/econ/arima.forecast.html. [53] The MathWorks. "Infer," 29 July 2015; http://www.mathworks.com/help/econ/arima.infer.html#outputarg_E. [54] Probability Methods Subcommittee, "IEEE Reliability Test System," IEEE Transactions on Power Apparatus and Systems, PAS-98, 1979. [55] A. Raue, C. Kreutz, T. Maiwald, J. Bachmann, M. Schilling, U. Klingmüller, and J. Timmer, “Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood,” BIOINFORMATICS, vol. 25, pp. 19231929, June 2009.