ANÁLISE CLÁSSICA DE SÉRIES TEMPORAIS E REDES NEURAIS ARTIFICIAIS PARA PREVISÃO DE CARGA A CURTO PRAZO: UM ESTUDO COMPARATIVO

UNIVERSIDADE FEDERAL FLUMINENSE ESCOLA DE ENGENHARIA DEPARTAMENTO DE ENGENHARIA ELÉTRICA

EDUARDO LUTTERBACH SAPORETTI AZEVEDO

ANÁLISE CLÁSSICA DE SÉRIES TEMPORAIS E REDES NEURAIS ARTIFICIAIS PARA PREVISÃO DE CARGA A CURTO PRAZO: UM ESTUDO COMPARATIVO

Niterói 2011

ii

EDUARDO LUTTERBACH SAPORETTI AZEVEDO

ANÁLISE CLÁSSICA DE SÉRIES TEMPORAIS E REDES NEURAIS ARTIFICIAIS PARA PREVISÃO DE CARGA A CURTO PRAZO: UM ESTUDO COMPARATIVO

Trabalho de Conclusão de Curso apresentado ao Curso de Graduação em Engenharia Elétrica da Universidade Federal Fluminense, como requisito parcial para obtenção do Grau de Engenheiro Eletricista. Área de Concentração: Sistemas de Potência.

Orientador: Prof. VITOR HUGO FERREIRA, D. Sc.

Niterói 2011

iii

iv

EDUARDO LUTTERBACH SAPORETTI AZEVEDO

ANÁLISE CLÁSSICA DE SÉRIES TEMPORAIS E REDES NEURAIS ARTIFICIAIS PARA PREVISÃO DE CARGA A CURTO PRAZO: UM ESTUDO COMPARATIVO

Trabalho de Conclusão de Curso apresentado ao Curso de Graduação em Engenharia Elétrica da Universidade Federal Fluminense, como requisito para obtenção do Grau de Engenheiro Eletricista.

Aprovada em julho de 2011.

BANCA EXAMINADORA

______________________________________________________ Prof. Vitor Hugo Ferreira, D. Sc. UFF

______________________________________________________ Prof. Júlio César Stacchini de Souza, D. Sc. UFF

______________________________________________________ Prof. Márcio Zamboti Fortes, D. Sc. UFF

Niterói 2011

v

SUMÁRIO 1

INTRODUÇÃO ....................................................................................................................1 1.1

2

ORGANIZAÇÃO DO TRABALHO .............................................................................3

PREVISÃO DE CARGA ......................................................................................................5 2.1

SÉRIES TEMPORAIS ...................................................................................................6

2.1.1

3

Processos Estocásticos ..................................................................................................................... 8

2.2

TRANSFORMAÇÕES ..................................................................................................9

2.3

TENDÊNCIA ...............................................................................................................10

2.4

SAZONALIDADE .......................................................................................................11

MODELO SARIMA ...........................................................................................................12 3.1

MODELOS LINEARES ESTACIONÁRIOS .............................................................14

3.2

MODELOS AUTO-REGRESSIVOS (AR) .................................................................16

3.3

MODELOS DE MÉDIAS MÓVEIS (MA) .................................................................17

3.4

MODELOS AUTO-REGRESSIVOS E DE MÉDIAS MÓVEIS (ARMA) ................19

3.5

FUNÇÃO DE AUTOCORRELAÇÃO PARCIAL ......................................................20

3.6

MODELOS NÃO-ESTACIONÁRIOS ........................................................................21

3.7

MODELOS ARIMA ....................................................................................................22

3.8

IDENTIFICAÇÃO DE MODELOS ARIMA ..............................................................23

3.9

MÉTODOS BASEADOS EM UMA FUNÇÃO PENALIZADORA ..........................24

3.9.1

Critério de informação de Akaike .................................................................................................. 25

3.9.2

Critério de informação Bayesiano ................................................................................................. 26

3.10 ESTIMAÇÃO DE MODELOS ARIMA......................................................................26 4

REDES NEURAIS ARTIFICIAIS .....................................................................................28 4.1

MODELO DE UM NEURÔNIO .................................................................................31

4.2

TIPOS DE FUNÇÃO DE ATIVAÇÃO .......................................................................33

4.2.1

Função de Limiar ........................................................................................................................... 33

4.2.2

Função Linear por Partes (Piecewise-linear) ................................................................................. 34

vi

4.2.3

4.3

ARQUITETURAS DE REDE .....................................................................................36

4.3.1

Redes Alimentadas Adiante (feed forward) com Camada Única .................................................. 36

4.3.2

Redes Feed Forward com Múltiplas Camadas .............................................................................. 37

4.3.3

Redes Recorrentes ......................................................................................................................... 39

4.4

PROCESSOS DE APRENDIZAGEM.........................................................................40

4.4.1

Aprendizagem com um Professor .................................................................................................. 41

4.4.2

Aprendizagem sem um Professor .................................................................................................. 43

4.5

5

Função Sigmóide ........................................................................................................................... 35

ETAPAS GERAIS DE PROJETOS DE REDES NEURAIS ......................................44

4.5.1

Definição do espaço de entrada ..................................................................................................... 45

4.5.2

Pré-processamento ......................................................................................................................... 45

4.5.3

Escolha da Estrutura ...................................................................................................................... 47

4.5.4

Treinamento do Modelo ................................................................................................................ 48

4.5.5

Avaliação do Modelo..................................................................................................................... 48

RESULTADOS ..................................................................................................................50 5.1

DADOS DO HISTÓRICO DA CARGA .....................................................................50

5.2

DEFINIÇÃO DOS MODELOS ...................................................................................52

5.2.1

Modelos SARIMA ......................................................................................................................... 52

5.2.2

Modelos Neurais ............................................................................................................................ 54

5.3

COMPARAÇÃO DOS RESULTADOS OBTIDOS ...................................................58

6

CONCLUSÃO ....................................................................................................................63

7

REFERÊNCIAS BIBLIOGRÁFICAS ...............................................................................65

vii

LISTA DE ILUSTRAÇÕES Figura 1 – Sistema Dinâmico ..................................................................................................... 6 Figura 2 – Um processo estocástico interpretado como uma família de variáveis aleatórias. (7) .................................................................................................................................................... 9 Figura 3 – Série não-estacionária quanto ao nível e inclinação. .............................................. 12 Figura 4 – Filtro linear, com entrada at , saída Z t e função de transferência  ( B) . .............. 15 Figura 5 – Um Neurônio Artificial ........................................................................................... 29 Figura 6 – Uma Rede Neural Alimentada Adiante com dupla camada ................................... 30 Figura 7 – Modelo não-linear de um neurônio. (retirado de (22)) ........................................... 31 Figura 8 – Transformação afim produzida pela presença de um bias (24)............................... 33 Figura 9 – Função de Limiar .................................................................................................... 34 Figura 10 – Função linear por partes (Piecewise-linear) ......................................................... 35 Figura 11 – Função Sigmóide para parâmetro de inclinação a variável ................................... 36 Figura 12 – Rede Feed Forward ou acíclica com única camada de neurônios ........................ 37 Figura 13 – Rede Feed Forward ou acíclica com uma única camada de neurônios ................ 38 Figura 14 – Rede recorrente sem laços de auto-realimentação e sem neurônios ocultos ......... 39 Figura 15 – Rede recorrente com neurônios ocultos ................................................................ 40 Figura 16 – Diagrama em blocos da aprendizagem com um professor.................................... 42 Figura 17 – Diagrama de blocos da aprendizagem não supervisionada ................................... 43 Figura 18 – Diagrama em Blocos das Etapas Gerais de um Projeto de Redes Neurais ........... 44 Figura 19 – Análise de Componentes Principais ...................................................................... 47 Figura 20 – Demandas Máximas Diárias de 1997 e 1998 (Dias do Ano x Demanda [MW]) .. 51 Figura 21 – Mês de Janeiro e Julho de 1997 ............................................................................ 51 Figura 22 – Rotina de Construção de um Modelo .................................................................... 52 Figura 23 – Rotina da Previsão da Rede Neural sem Sazonalidade ......................................... 55 Figura 24 – Rotina da Previsão com Análise Espectral para modelagem da Sazonalidade ..... 56 Figura 25 – Rotina da Previsão por meio de Variáveis Binárias .............................................. 57 Figura 26 – Previsões de Janeiro .............................................................................................. 62 Figura 27 – Distribuição dos Competidores de Cada País ....................................................... 62

viii

LISTA DE TABELAS

Tabela 1 – Previsão dos Modelos em Estudo ........................................................................... 59 Tabela 2 – MAPE, MAE e Tempo para Previsão..................................................................... 60 Tabela 3 – Resultado dos 10 Primeiros Colocados na competição da EUNITE. ..................... 61

ix

RESUMO A impossibilidade de armazenamento, em escala conveniente, da energia elétrica acarreta a necessidade de expansão do sistema a frente da demanda. Investimentos antecipados onerariam desnecessariamente ou até inviabilizariam a expansão do sistema. Ser preciso na previsão da carga possibilita o planejamento justo à necessidade. A tomada de decisões referentes à operação também sofre expressiva influência do comportamento futuro da carga e a escolha de modelos que dispõem de segurança na previsão trazem maior estabilidade e confiança. Neste trabalho é comparado o desempenho da análise clássica de séries temporais e de modelo neural à previsão de carga a curto prazo. SARIMA (AutoRegressivos Integrados de Médias Móveis Sazonal) é o modelo utilizado na análise clássica, em que são feitas duas abordagens, uma com a escolha dos parâmetros feita graficamente pelo usuário e a outra autônoma através de métodos baseados em uma função penalizadora, pelo critério de informação Akaike (AIC) e o Bayesiano (BIC). Para os modelos neurais são comparadas diversas formas de processamento dos dados, incluindo análise espectral via séries de Fourier e ainda a sazonalidade processada por rede neural com entradas binárias. É comparada também a possibilidade de seleção automática das entradas dos modelos neurais utilizando teoria do caos contra a seleção dos atrasos feita pelo usuário. Além da forma de especificação dos padrões de entrada, um método de agrupamento automático dos dados é testado, visando produzir modelos locais de previsão e sua comparação contra abordagem global, ou seja, sem o uso de métodos de agrupamento para definição do conjunto de treinamento da rede neural. Os resultados obtidos pelas previsões com essas características foram comparados entre si e com os de uma competição internacional para fins de avaliação da eficácia dos modelos.

x

ABSTRACT The impossibility of electricity energy storage in appropriate scale entails the need for expansion of the system ahead of demand. Anticipated investments unnecessarily onerous or even made it impossible to expand the system. Be accurate in predicting the load enables the planning just to the need. The decisions concerning the operation also suffers significant influence of future behavior of the load and the choice of models that have security in predicting bring stability and confidence. This paper compared the performance of classical time series analysis and neural model to forecast short-term load. SARIMA (Auto-Regressive Integrated Moving Average seasonal) is the model used in classical analysis, in which two approaches are made, with a choice of parameters made graphically by the user and by other independent methods based on a penalty function, the criterion Akaike information (AIC) and Bayesian (BIC). For the neural models are compared several ways of data processing, including spectral analysis via Fourier series and also seasonality processed by neural network with binary inputs. It also compared the possibility of automatic selection of the inputs of neural models using chaos theory against the delays selection made by the user. Addition to the specifying form of input patterns, a method of automatic data clustering is tested, to produce local models for prediction and its comparison with the global approach, i.e. without the use of clustering methods for defining the set of training neural network. The results obtained by the predictions with these characteristics were compared with other and with an international competition for the purpose of evaluating the effectiveness of the models.

1

INTRODUÇÃO A vontade de antever acontecimentos é inerente ao ser humano. A possibilidade de

previsão de comportamentos desperta interesse ou, no mínimo, curiosidade. Em muitos campos de estudo essa vontade não passa a ser uma necessidade. Em outros há a necessidade de se estimar o porvir, é utilizada a palavra estimar, pois, exatidão perfeita, não é passível de ser atingida. Em sistemas de potência a estimativa da carga futura, imprescindível ao bom funcionamento do sistema, é conhecida como previsão de carga. Na previsão de carga os intervalos de tempo para previsão são usualmente [1]: 

curto prazo, as previsões são feitas de alguns minutos ou até um mês, e podem ser dividas em minutos, horas ou dias;



médio prazo, as previsões são de uma semana até cinco anos, podendo ser divididas em dias, semanas ou meses;



longo prazo, as previsões são de meses até trinta anos, dividas mensalmente ou anualmente.

Dependendo do estudo a ser feito, estes horizontes podem variar (o que é considerado médio para uns pode ser longo para outros), portanto essa definição não é imutável. Segundo [1], tais horizontes de previsão são úteis, necessários e utilizados de muitas maneiras e em diversos setores como pode ser visto a seguir. As previsões em longo prazo são utilizadas em funções relacionadas com o planejamento técnico e financeiro das empresas de energia, tais como planejamento da expansão dos sistemas de transmissão e distribuição e do parque gerador, programação anual da manutenção de unidades geradoras, gerenciamento energético de longo prazo, desenvolvimento de estratégias operacionais, estudos de viabilidade econômica, planejamento dos investimentos e do orçamento, e pesquisa de mercado. Em mercados regulamentados, tais previsões podem ser utilizadas para o desenvolvimento de políticas tarifárias.

2

As previsões em médio prazo fornecem subsídios para diversas atividades relacionadas ao planejamento da expansão e da operação de sistemas de potência, podendo ser citados: programação da compra de combustíveis; planejamento da manutenção de equipamentos, do intercâmbio entre áreas, das transações energéticas e do orçamento; otimização da programação das unidades geradoras; e desenvolvimento de estratégias de gerenciamento energético. Tais previsões também podem ser utilizadas também para desenvolvimento de políticas tarifárias. Esse trabalho lidará com a previsão em curto prazo. As previsões para este horizonte são de suma importância para a operação e o controle em tempo real de sistemas de potência. Dentre as funções inerentes ao planejamento da operação, de acordo com [1], a previsão de carga em curto prazo fornece subsídios para análise de segurança, incluindo estudo de contingências e elaboração de estratégias de gerenciamento da carga; programação da geração, abrangendo coordenação hidrotérmica, programação da compra e alocação de combustível, comissionamento de unidades térmicas e despacho econômico; estudos de fluxo de potência, como fluxo de potência ótimo e programação do intercâmbio entre áreas; programação da alocação de reserva girante; programação e avaliação das transações de compra e venda de energia; e programação da manutenção. Dentre as atividades relacionadas ao controle em tempo real de sistemas de potência, estas previsões fornecem informações importantes para controle automático da geração e controle do fluxo de potência reativa. Para o caso específico de empresas de distribuição de energia, o conhecimento do comportamento futuro da carga, particularmente do pico de carga, nas diversas barras do sistema é um dos requisitos mais importantes para o aumento da eficiência da operação. Estas informações constituem a base para a estimação do estado do sistema e para cálculos técnicos e econômicos, possibilitando assim melhorias na manutenção dos equipamentos e no planejamento da operação dos sistemas de distribuição. Tais melhorias podem ser obtidas através de instalações de equipamentos de emergência, desligamento de circuitos, transferências de carga, aumento da refrigeração de equipamentos críticos e ajuste dos tap’s dos transformadores das subestações. Exposta a importância sob o ponto de vista técnico, ainda deve ser levado em consideração, o grande interesse suscitado na área econômica. A exigência de que a oferta de energia elétrica deva estar sempre a frete da demanda [2], e o lapso do tempo que decorre entre a decisão de construir uma instalação de geração, transmissão ou distribuição elétrica e o

3

momento de sua entrada em operação, fazem que seja necessário projetar a demanda com uma antecipação superior à de outros setores da economia. Este problema não existiria se fosse possível simplesmente construir novas instalações a um ritmo tão elevado, que o avanço da demanda estivesse ultrapassado. Logicamente este método significa, ter, na prática, uma margem excessiva de reserva a um volume de investimentos que a maioria dos países não suportaria, pois aquelas instalações estariam funcionando permanentemente com um grau de utilização inferior ao ótimo desejável. Nos países que privatizaram seu mercado de energia, onde, além disso, há a competição por melhores preços, os agentes de mercado procuram atingir elevados níveis de eficiência com a redução nas margens da operação. Segundo a estimativa apresentada em [3], para empresas de energia que apresentem gastos com combustível da ordem de centenas de milhões de dólares anuais, melhorias da ordem de 1 (um) % na precisão das previsões em curto prazo podem resultar em reduções nos custos operacionais da ordem de centenas de milhares de dólares por ano. Outros estudos mostrando o impacto da precisão das previsões nos custos das empresas de energia podem ser encontrados em [4][5][6]. Notada a inserção da previsão de carga em curto prazo nas diversas atividades relacionadas com a operação de sistemas de potência, a precisão de tais previsões está intimamente ligada à redução dos custos operacionais das empresas de energia. Já que uma previsão de carga é crucial sob o ponto de vista técnico e econômico, muitos estudos vêem sendo feitos neste campo nas últimas décadas. Várias metodologias foram criadas, porém, quase todas necessitam de intervenção humana especializada em suas variadas fases de modelagem. Quando se está de frente à pequenos sistemas com uma ou algumas barras a intervenção, ainda que se gaste algum tempo, é plausível. Porém é inviável (até impossível), quando o estudo toma uma proporção maior, como no sistema elétrico brasileiro que possui mais de cinco mil barras. Um modelo de previsões autônomas se fez necessário. 1.1

ORGANIZAÇÃO DO TRABALHO Este trabalho está organizado em cinco Capítulos, resumidos a seguir. O Capítulo II irá tratar das características inerentes à energia elétrica que a difere de

outros tipos de energia, da inserção da previsão de carga no contexto do sistema elétrico brasileiro (SIN) e mostrará como esta previsão pode ser analisada como uma série temporal.

4

Ainda será vista a definição de processo estocástico e a importância da análise no tratamento da série. O Capítulo III trará uma metodologia (SARIMA - Auto-Regressivos Integrados de Médias Móveis Sazonal ou em inglês Seasonal Autoregressive integrated moving average) para parametrização na análise de modelo de séries temporais e irá mostrar como esta lida com a previsão de carga. O Capítulo IV descreve a inspiração para a criação de um modelo neural e a motivação para o seu desenvolvimento, detalha as etapas para sua implantação e os cuidados a serem tomados para que se obtenha sucesso no tratamento do problema. O Capítulo V mostra os resultados obtidos com as técnicas de previsão de carga utilizadas incluindo o seu tempo de resposta, compara a eficácia entre elas e ainda faz um paralelo com uma competição internacional. Será feito o uso de um histórico da carga de dois anos anteriores ao mês de previsão. No Capítulo VI serão expostas as conclusões retiradas sobre o estudo.

2

PREVISÃO DE CARGA Etimologicamente (prae e videre) [7], a palavra previsão sugere que se quer ver uma

coisa antes que ela exista. Alguns autores preferem a palavra predição, para indicar algo que deverá existir no futuro. Ainda outros utilizam o termo projeção. Neste trabalho será usada a palavra previsão com o sentido indicado acima. É importante salientar que a previsão não constitui um fim em si, mas apenas um meio de fornecer informações para uma consequente tomada de decisões, visando determinados objetivos. A principal característica da energia elétrica que a difere de outras formas de energia é a incapacidade de ser armazenada em quantidades significantes. Pequenas quantidades podem ser armazenadas através, por exemplo, de baterias, porém, estas seriam ínfimas perante a energia utilizada no sistema. O Operador Nacional do Sistema (ONS) traz no módulo 5 do Procedimentos de Rede a “Consolidação da previsão de carga” que estabelece a sistemática do fornecimento da previsão da carga para os estudos realizados pelo mesmo, que abrangem desde as ampliações e reforços na rede básica até os que tratam do planejamento da operação elétrica e energética, em seus diversos horizontes [8]. O processo de consolidação da previsão de carga está inserido nas atividades de planejamento da operação atribuídas ao ONS. Devido ao fato de não poder ser armazenada, a energia elétrica deve ser gerada concomitantemente ao seu consumo. Deste fazem parte, tanto a energia requerida pelos consumidores ligados ao sistema, quanto a energia nas perdas inerentes ao seu transporte. Portanto, como era de se esperar, a geração seguirá um padrão lógico, já que não existe a possibilidade de saber, precisamente, de antemão o comportamento da carga. Este comportamento é caracterizado como uma série temporal.

6

2.1

SÉRIES TEMPORAIS Segundo [7], uma série temporal é um conjunto de observações ordenadas, yt , cada

uma observada em um instante de tempo. As séries são divididas em dois tipos: contínuas, em que o tempo varia de forma contínua, e discretas, em que o tempo varia discretamente. Há, basicamente, dois enfoques usados na análise de séries temporais. Em ambos, o objetivo é construir modelos para as séries, com propósitos determinados. No primeiro enfoque, a análise é feita no domínio temporal (interesse na magnitude dos eventos que ocorrem em determinado instante do tempo) e os modelos propostos são modelos paramétricos (com um número finito de parâmetros). No segundo, a análise é conduzida no domínio de frequência (interesse na frequência com que certos eventos ocorrem em um período de tempo) e os modelos propostos são modelos não-paramétricos. Dentre os modelos paramétricos tem-se, por exemplo, os modelos ARIMA (Autoregressive integrated moving average – modelos auto-regressivos integrados e de médias móveis), que serão analisados nesse trabalho levando em conta a sazonalidade dos eventos (SARIMA, Seasonal Autoregressive integrated moving average – modelos auto-regressivos integrados e de médias móveis sazonais). No domínio da frequência tem-se a análise espectral. O objetivo das análises de séries temporais é, segundo [7] analisar um conjunto de dados (ao longo do tempo), procurar periodicidades relevantes nos dados e selecionar e estimar um modelo estocástico que possivelmente tenha gerado o conjunto de dados. E são dois os principais motivos por utilizar estes modelos: entender o mecanismo gerador da série temporal e predizer o comportamento futuro do sistema. Independente do modelo a ser escolhido, este deve ser simples e parcimonioso, ou seja, o número de parâmetros envolvidos deve ser o menor possível e sua utilização, dentro do possível, não deve suscitar dificuldades às pessoas que o irão manipular. Em muitas situações existe uma relação entre entrada e saída de dados. Relação esta que pode ter inúmeras características. É chamado sistema dinâmico o processo caracterizado por um série de entrada X  t  , uma série de saída Z  t  e uma função de transferência   t  (Figura 1).

X (t )

 (t )

Z (t )

Figura 1 – Sistema Dinâmico

7

De particular importância são os sistemas lineares, onde a saída é relacionada com a entrada através de um funcional linear envolvendo   t  . Um exemplo típico é 

Z  t     X  t   

(2.1)

 0

Problemas de interesse aqui são: a) estimar a função de transferência   t  , conhecendo-se as séries de entrada e saída; b) fazer previsões da série Z  t  , com o conhecimento de observações da série de entrada X  t  ; c) estudar o comportamento do sistema, simulando-se a série de entrada; d) controlar a série de saída Z  t  , de modo a trazê-la o mais próximo possível de um valor desejado, ajustando-se convenientemente a série de entrada X  t  ; este controle é necessário devido a perturbações que normalmente afetam um sistema dinâmico. O comportamento mais esperado de uma série temporal é a estacionaridade, ou seja, que ela apresente, ao longo do tempo, mudanças aleatórias sem que se distancie de um valor médio constante, permanecendo estável quando o tempo é grande o suficiente para tornar suas oscilações imperceptíveis. Entretanto, a maioria das séries observadas, na prática, são nãoestacionárias e apresentam tendências, sendo a forma mais simples, destas, quando as oscilações flutuam em torno de uma reta com inclinação positiva ou negativa com tendência linear. Outra forma não-estacionária é a explosiva, um exemplo é o crescimento ou decrescimento exponencial. Os modelos utilizados para descrever séries temporais são processos estocásticos, isto é, processos controlados por leis probabilísticas. Sendo

Zt , t  1, 2,

, N  observações de uma série temporal. Um modelo de

decomposição consiste em escrever Z  t  como uma soma de três componentes nãooberváveis,

Zt  Tt  St  at

(2.2)

8

onde Tt e St representam a tendência e sazonalidade, respectivamente, enquanto at é uma componente aleatória de média e variância constante  a2 . Supondo at  como um processo estacionário então Z t  será uma série não estacionária. O interesse principal em considerar um modelo do tipo (2.2) será o de estimar St e construir a série livre de sazonalidade ou sazonalmente ajustada. Isto é, se Sˆt for uma estimativa de St ,

ZtSA  Zt  Sˆt

(2.3)

será a série sazonalmente ajustada. Há várias razões para considerar este procedimento de ajustamento sazonal. As componentes Tt e St são, em geral, bastante relacionadas e a influência da tendência sobre a componente sazonal pode ser muito forte, segundo [9], por duas razões: a) método de estimação de St podem ser bastante afetados se não se levar em conta a tendência; b) a especificação de St depende da especificação de Tt . Por isso, não poderá se isolar uma das componentes sem tentar isolar a outra. Estimando-se Tt e St e subtraindo de Z  t  obtém-se uma estimativa da componente aleatória at . 2.1.1

Processos Estocásticos

Definição de [7]: Seja T um conjunto arbitrário. Um processo estocástico é uma família Z  Z (t ), t  T  , tal que, para cada t  T , Z  t  é uma variável aleatória. Nestas condições, um processo estocástico é uma família de variáveis aleatórias (v.a.), que serão supostamente definidas um mesmo espaço de probabilidades  , ,   . O conjunto T é normalmente tomado como o conjunto dos inteiros Z  0, 1, 2,... ou o conjunto dos reais ℝ. Também, para cada t  T , Z  t  será uma v.a. real. Como, para t  T , Z  t  é uma v.a. definida sobre  , na realidade Z  t  é uma função de dois argumentos, Z  t ,   , t  T ,   . A Figura 2 ilustra esta interpretação de processo estocástico.

9

Figura 2 – Um processo estocástico interpretado como uma família de variáveis aleatórias. [7] Para cada t  T , tem-se uma v.a. Z  t ,   , com uma distribuição de probabilidades; é possível que a função densidade de probabilidade (fdp) no instante t1 seja diferente da fdp no instante t2 , para dois instantes t1 e t2 quaisquer, porém, o mais comum é aquela em que a fdp de Z  t ,   é a mesma, para todo t  T . Para cada   fixado, obterem-se uma série temporal, ou seja, uma realização ou trajetória do processo que será uma função de t . Os métodos de previsão podem ser divididos em duas categorias: 

automáticos: que são aplicados diretamente com a utilização de um computador;



não-automáticos: exigem a intervenção de pessoal especializado, para serem aplicados.

2.2

TRANSFORMAÇÕES Existem dois motivos essenciais que levam a transformar os dados originais:

estabilizar a variância e tornar o efeito sazonal aditivo. Para tornar uma série estacionária serão tomadas diferenças. Para se ter uma idéia do tipo de transformação que é adequada, [7] mostra que

pode ser utilizado um gráfico que traz no eixo das abscissas médias de

subconjuntos de observações da série original e no eixo das ordenadas a amplitude de cada um desses subconjuntos; se Z1 ,

, Z k for um tal subconjunto com k observações, calcula-se

10

Z

1 k  Zt , k i 1 i

  max  Zt   min  Zt  , i

i

que são medidas de posição e variabilidade, respectivamente; o par  Z ,   será um ponto do gráfico. O número de elementos em cada sub-série pode ser igual ao período, no caso de séries sazonais. Se  independer de Z , obteria-se pontos espalhados ao redor de uma reta paralela ao eixo das abscissas, e neste caso não haverá necessidade de transformação. Se, no entanto,  for diretamente proporcional a Z , a transformação logarítmica é apropriada. Uma outra razão para efetuar transformações é obter uma distribuição para os dados mais simétrica e próxima da normal. 2.3

TENDÊNCIA Inicialmente irá se supor que a componente sazonal St não esteja presente. O modelo

considerado será:

Zt  Tt  at

(2.4)

onde at é ruído branco, com variância  a2 . Há vários métodos para estimar Tt . De acordo com [7], os mais utilizados consistem em: (i).

ajustar uma função do tempo, como um polinômio, uma exponencial ou outra função suave de t ;

(ii).

suavizar (ou filtrar) os valores da série ao redor de um ponto, para estimar a tendência naquele ponto;

(iii).

suavizar os valores da série através de sucessivos ajustes de retas de mínimos quadrados ponderados (“lowess”).

Estimando-se a tendência através de Tˆt , pode-se obter a série ajustada para tendência ou livre de tendência,

Yt  Zt  Tˆt . Um procedimento que é também utilizado para eliminar a tendência de uma série é aquele de tomar diferenças. Normalmente a primeira diferença

11

Zt  Zt  Zt 1 já é estacionária. 2.4

SAZONALIDADE [10] interpreta sazonalidade como sendo um movimento regular de uma série dentro

de um período, ou ainda pode ser interpretada como a representação de movimentos sistemáticos. O objetivo será ajustar uma série para que esta fique livre de sazonalidade. Deverá se estimar St e subtrair a série estimada de Z t no modelo (2.2) ,

Zt  Tt  St  at ,

t  1, 2,

,N .

Desta maneira, um procedimento de ajustamento sazonal consiste em: a) obter estimativas Sˆt de St ; b) calcular

ZtSA  Zt  Sˆt .

É necessário estimar tanto a tendência quanto a sazonalidade que, como estão intimamente associadas, não se conseguirá isolar uma da outra. Estimando-se Tt e St e subtraindo de Z t se obterá uma estimativa da componente aleatória at .

3

MODELO SARIMA Uma metodologia foi proposta por [11] para parametrização na análise de modelo de

séries temporais [7]. Essa metodologia é chamada ARIMA (p,d,q), e consiste em modelos ajustados auto-regressivos integrados de médias móveis, a um conjunto de dados. A classe dos modelos ARIMA conseguirá descrever séries estacionárias e não estacionárias de maneira satisfatória desde que não apresentem característica explosiva. Quando a série apresenta esse comportamento é dita homogênea, mesclando o comportamento estacionário com o não-estacionário com tendência. A Figura 3 mostra esta forma de não-estacionaridade.

Figura 3 – Série não-estacionária quanto ao nível e inclinação. Devido ao fato de a maioria dos procedimentos de análise estatística de séries temporais intuir estacionaridade das mesmas, os dados de origem devem ser transformados de modo a serem vistos como estacionários se estes ainda não o forem. A transformação mais corriqueira consiste em tomar diferenças sucessivas da série original, até o ponto em que estas diferenças sejam estacionárias. A primeira diferença de Z  t  é definida por

13

Z (t )  Z (t )  Z (t 1)

A segunda diferença será

2 Z (t )    Z (t )    Z (t )  Z (t  1) , Ou seja,

2 Z (t )  Z (t )  2Z (t 1)  Z (t  2) . De modo geral, a n -ésima diferença de Z  t  é  n Z (t )    n1Z (t ) 

Normalmente, tomar uma ou duas diferenças será suficiente para tornar a série estacionária. A classe de modelos SARIMA, avaliada nesse trabalho, é capaz de descrever séries estacionárias e não estacionárias, e ainda as que possuem sazonalidade, de maneira satisfatória, mas com o pressuposto de não possuir comportamento explosivo. Uma metodologia bastante utilizada na análise de modelos paramétricos é conhecida como abordagem de [11]. Tal metodologia consiste em ajustar modelos auto-regressivos integrados de médias móveis, ARIMA (p, d, q), a um conjunto de dados. Três classes de processos podem ser descritos pelos modelos ARIMA: 

processos lineares estacionários;



processos lineares não-estacionários homogêneos;



processos de memória longa.

Processos descritos de maneira adequada pelos chamados modelos auto-regressivos e de médias móveis de ordens p, d e q: ARIMA (p, d e q), e podem ser generalizados pela inclusão de um operador sazonal. A estratégia para a construção do modelo será baseada em um ciclo iterativo, no qual a escolha da estrutura do modelo é baseada nos próprios dados. Os estágios do ciclo iterativo são: 

uma classe geral dos modelos é considerada para análise (especificação);



há identificação de um modelo, com base na análise de autocorrelações, autocorrelações parciais e outros critérios;



a seguir vem a fase de estimação, na qual os parâmetros do modelo identificado são estimados;

14



finalmente, há a verificação ou diagnóstico do modelo ajustado, através de uma análise de resíduos, para saber se este é adequado para os fins em vista (previsão, por exemplo).

Caso o modelo não seja adequado, o ciclo é repetido, voltando-se à fase de identificação. Geralmente, utiliza-se a identificação de mais de um modelo, procedendo a suas respectivas estimações e especificações a fim de escolher entre estes o melhor ajustado. Na previsão, por exemplo, o que fornecer o menor erro quadrado médio. A identificação é o que levanta maiores divergências. Pesquisadores podem identificar modelos diferentes para a mesma série temporal. Portanto esta é a fase crítica do procedimento mencionado. A notação de operadores será extensivamente utilizada neste trabalho. Esta notação será útil para entendimento e facilitará a manipulação do modelo estudado. Estes operadores são: (a) operador translação para o passado, denotado por B e definido por

BZt  Zt 1 , B m Zt  Zt m ;

(3.1)

(b) operador translação para o futuro, será chamado de F e definido por

FZt  Zt 1 , F m Zt  Zt m ;

(3.2)

Zt  Zt  Zt 1  1  B  Zt .

(3.3)

(c) operador diferença,

(d) operador soma, denotado por S e definido por 

SZt   Zt  j  Zt  Zt 1  j 0

1  B  B

2



Z , t

(3.4)

portanto (3.5) SZt  1  B  Zt   1Zt

(3.6)

S   1 .

(3.7)

1

ou seja,

3.1

MODELOS LINEARES ESTACIONÁRIOS

Todos os modelos encontrados neste estudo serão casos particulares de um modelo de filtro linear. Estes modelos supõe que a série temporal seja gerada através de um filtro linear (ou sistema linear), cuja entrada, é um ruído branco, Figura 4.

15

Figura 4 – Filtro linear, com entrada at , saída Z t e função de transferência  ( B) . As saídas serão, formalmente, dadas por

Zt    at  1at 1  2 at 2 

    B  at ,

(3.8)

Em que

  B   1  1 B  2 B 2 

(3.9)

é denominada função de transferência do filtro e μ é um parâmetro determinando o nível da série. Sendo Z t fornecido por (3.8) é um processo linear (discreto). Vale lembrar que E  at   0, t , Var  at    a2 , t , E  at as   0, s  t.

Chamando Zt  Zt   , tem-se que

Zt    B  at .

(3.10)

Se a sequência de pesos  j , j  1 for finita ou infinita e convergente, o filtro é estável e Z t é estacionária. Neste caso, µ é a média do processo. Caso contrário, Z t é nãoestacionária e µ não tem significado específico, a não ser como um ponto de referência para o nível da série. De (3.8) tem-se que    E  Zt     E  at   j at 1  j 1  

e como E  at   0 , para todo t, tem-se que E  Zt    se a série





 j convergir.

j 1

Z t pode ser escrito em uma forma alternativa, como uma soma ponderada de valores passados Zt 1 , Zt 2 ,

mais um ruído at :

Zt  1Zt 1   2 Zt 2 



 at    j Zt  j  at j 1

(3.11)

16

3.2

MODELOS AUTO-REGRESSIVOS (AR) Se em (3.11)  j  0 , j  p , será obtido um modelo auto-regressivo de ordem p,

denotado por AR  p 

Zt  1Zt 1  2 Zt 2 

  p Zt  p  at ,

(3.12)

renomeando os pesos de  j para  j , de acordo com a notação usual. Se for definido o operador auto-regressivo estacionário de ordem p

 ( B)  1  1B  2 B2 

 p B p ,

(3.13)

então, pode-se escrever

  B  Zt  at .

(3.14)

As principais características de um processo representado pelo modelo AR  p  , serão vistas a seguir. A) Estacionaridade e Invertibilidade Como  ( B)   ( B)  1  1B  ...   p B p é finito, não há restrições sobre os parâmetros para assegurar invertibilidade de Z t . Sejam Gi1 , i  1,

, p , as raízes da equação característica   B   0 ; então pode-se

escrever

 ( B)  (1  G1B)(1  G2 B)

(1  Gp B)

e, expandindo-se em frações parciais, 

Ai . i 1 1  Gi B

  B    1  B   

Se   B  deve convergir para B  1 deve-se ter que Gi  1 , i  1,

(3.15) , p . Esta

condição é equivalente a que a equação característica   B   0 tenha raízes fora do círculo unitário. Esta é a condição de estacionaridade. B) Função de autocorrelação Multiplicando-se ambos os membros de (3.12) por Z t  j e tomando-se a esperança obtém-se

E (Zt Zt  j )  1E (Zt 1Zt  j )  2 E (Zt 2 Zt  j ) 

  p E (Zt  p Zt  j )  E (at Zt  j ) .

17

Devido a (3.12), Z t  j só envolve ruídos até at  j , não-correlacionados com at ,





E at Zt  j  0 , j  0 , do que resulta

 j  1 j 1  2 j 2 

  p j  p , j  0 .

(3.16)

 p  j p , j  0 ,

(3.17)

Dividindo-se por  0  Var  Zt  , obtém-se

 j  1 j 1  2  j 2  que também pode ser escrita

 ( B)  j  0 ,

(3.18)

onde o operador B agora age em j : B j   j 1 , etc. Se p

 ( B)   1  Gi B  , i 1

então pode-se demonstrar que a solução geral de (3.17) é, no caso de raízes distintas, j j  j  AG 1 1  A2G2 

 ApGpj

(3.19)

Como Gi  1 , duas situações pode ocorrer: (a) se Gi for real, o termo Ai Gi j decai geometricamente para zero (amortecimento exponencial); (b) um par de raízes complexas conjugadas contribui com um termo da forma

Ad j sen(2 fj   ) (senóide amortecida). Genericamente, a função de autocorrelação de um processo auto-regressivo é constituída de uma mistura de polinômios, exponenciais e senóides amortecidas. 3.3

MODELOS DE MÉDIAS MÓVEIS (MA) A partir de (3.8) e supondo  j  0 , j  q ; tem-se um processo de médias móveis de

ordem q, que será denotado por MA(q) (abreviação do inglês “moving average”). Posteriormente será utilizado da forma Zt    at  1at 1 

 q at q

(3.20)

e sendo Zt  Zt   , obtém-se ~

Z t  (1  1B 

 q B q )at   ( B)at

(3.21)

onde

 ( B)  1  1B  2 B2 

 p B p

(3.22)

18

é o operador de médias móveis de ordem q. Serão vistas, a diante, as principais características do modelo MA (q). A) Estacionaridade e Invertibilidade Dado que   B   1  1B 

 q Bq , não há restrições sobre os parâmetros  j para

que o processo seja estacionário. Utilizando-se de um argumento completamente similar ao que foi feito para um modelo AR (p), no caso de estacionaridade, é verificado que a condição de invertibilidade para um modelo MA (q) é que as raízes da equação característica   B   0 estejam fora do círculo unitário. Nestas condições, um modelo MA (q) é equivalente a um modelo AR de ordem infinita. B) Função de autocorrelação De (3.20) tem-se que a facv é

 j  E Z t Z t  j 

q q     E   at   k at k   at  j  l at  j l   k 1 l 1   

 E  at at  j    k E  at  j at k   l E  at at  j l    kl E at k at  j l  q

q

q

k 1

l 1

q

k 1 l 1

Lembrando que

 2 , j  0,  a  j   E (at at  j )   a

 0, j  0,

logo,

 0  Var(Zt )   Z2  (1  12 

 q2 ) a2 .

(3.23)

em termos de  a  j  tem-se que a facv de Z t fica q

q

k 1

l 1

q

q

 j   a ( j )   k  a (k  j ) l a (l  j )   kl a ( j  l  k ) k 1 l 1

do que resulta 

q j



 j    j  l j l   a2 , j  0, l 1    ( j  1 j 1   2 j  2 

  q q  j ) a2 , j  1,

 0, j  q.

De (3.23) e (3.24) obtém-se a fac do processo,

, q,

(3.24)

19

  j  1 j 1   2 j  2    q  j q , j  1,  1  12   22    q2 j    0, j  q 

, q,

(3.25)

Nota-se, desta maneira, que a fac de um processo MA(q) é igual a zero para “lags” maiores do que q, de modo inverso ao que acontece com um processo AR. 3.4

MODELOS AUTO-REGRESSIVOS E DE MÉDIAS MÓVEIS (ARMA) Em algumas áreas de estudo, onde se deseja estimar o valor de alguma variável no

instante t como função de valores defasados da mesma variável, são bastante utilizados os modelos auto-regressivos. Para muitas séries observadas no dia-a-dia, se, se deseja um modelo com um número não muito grande de parâmetros, a inclusão de termos autoregressivos e de médias móveis é a solução adequada. Aparecem, então, os modelos ARMA (p,q), da forma

Zt  1Zt 1  Se  ( B)

e  ( B)

  p Zt  p  at  1at 1 

 q at q .

(3.26)

são os operadores auto-regressivos e de médias móveis,

respectivamente, introduzidos anteriormente, pode ser escrita (3.26) na forma compacta

 ( B)Zt   ( B)at .

(3.27)

A) Estacionaridade e Invertibilidade O processo é estacionário se as raízes de   B   0 caírem todas fora do círculo unitário e o processo é invertível se todas as raízes de   B   0 caírem fora do círculo unitário. B) Função de autocorrelação ~

Multiplicando ambos os membros de (3.26) por Z t  j e tomando esperanças, obtém-se

 j  E(Zt Zt  j )  E{(1Zt 1 

  p Zt  p  at  1at 1 

 q at q )Zt  j } ,

ou seja,

 j  1 j 1  2 j 2 

  p j  p   za ( j )  1 za ( j  1) 

 q za ( j  q) ,

(3.28)

onde  za  j  é a variância cruzada entre Z t e at , definida por

 za  j   E  at Z t  j  ~



~



Como Z t  j só depende de choques at ocorridos até o instante t  j , obtém-se

(3.29)

20

 0, j  0,  0, j  0.

 za  j  

(3.30)

De modo que (3.28) fica

 j  1 j 1  2 j 2 

  p j  p , j  p

(3.31)

 p  j p , j  q ,

(3.32)

A fac é obtida de (3.31):

 j  1 j 1  2  j 2 

e, portanto, as autocorrelações de “lags” 1, 2,..., q serão afetadas diretamente pelos parâmetros de médias móveis, mas para j  q as mesmas se comportam como nos modelos autoregressivos. Pode ser verificado que se q < p a fac consiste numa mistura de exponenciais e/ou de senóides amortecida; entretanto, se q  p , os primeiros q – p – 1 valores 0 , 1 ,

, q  p

não seguirão este padrão. 3.5

FUNÇÃO DE AUTOCORRELAÇÃO PARCIAL Os processos AR(p), MA(q) e ARMA (p,q) apresentam fac com características

especiais. Assim (i)

um processo AR(p) tem fac que decai de acordo com exponenciais e/ou senóides amortecidas, infinita em extensão;

(ii)

um processo MA(q) tem fac finita, no sentido que ela apresenta um corte após o “lag” q;

(iii)

um processo ARMA(p,q) tem fac infinita em extensão, a qual decai de acordo com exponenciais e/ou senóides amortecidas após o “lag” q  p .

Estas observações serão úteis no procedimento de identificação do modelo aos dados observados; calculando-se as estimativas das fac , que acredita-se reproduzir adequadamente as verdadeiras fac desconhecidas e comparando seu comportamento com o descrito acima, para cada modelo, tentará se escolher um (ou mais) modelo (modelos, respectivamente) que descreva(m) o processo estocástico. [12] propõem a utilização de um outro instrumento para facilitar este procedimento de identificação: a função de autocorrelação parcial (facp). Será denotado por kj o j-ésimo coeficiente de um modelo AR(k), de tal modo de kk seja o último coeficiente. Sabe-se que

 j  k1 j 1  k 2  j 2 

 kk  j k , j  1,

, k,

(3.33)

21

de [10], operando os coeficiente e colocando-os na forma matricial , obtém-se:

kk 

Pk* Pk

,

(3.34)

onde Pk é a matriz de autocorrelações e Pk* é a matriz Pk* com a última coluna substituída pelo vetor de autocorrelações. A quantidade kk , encarada como função de K, é chamada função de autocorrelação parcial (facp). Pode-se demonstrar [12] que, para os processos estudados, tem-se: (i)

um processo AR(p) tem facp kk  0 , para k  p e kk  0 , para k  p ;

(ii)

um processo MA(q) tem facp que se comporta de maneira similar à fac de um processo AR(p): é denominada por exponenciais e/ou amortecidas;

(iii)

um processo ARMA(p,q) tem facp que se comporta como a facp de um processo MA puro.

Será necessário estimar a facp de um processo AR, MA ou ARMA. Uma maneira consiste em estimar, sucessivamente, modelos auto-regressivo de ordens p  1, 2,3,

por

mínimos quadrados e tomar as estimativas do último coeficiente de cada ordem. 3.6

MODELOS NÃO-ESTACIONÁRIOS Os modelos vistos até agora (lineares estacionários) são indicados para descrever

séries estacionárias, isto é, séries que se desenvolvem no tempo ao redor de uma média constante. Muitas séries encontradas na prática não são estacionárias. Várias séries são nãoestacionárias, mas quando diferençadas tornam-se estacionárias. Os modelos tratados neste trabalho são apropriados para representar séries cujo comportamento seja não explosivo, em particular séries que apresentam alguma homogeneidade em seu “comportamento não-estacionário”. Séries Z t tais que, tomando-se um número finito de diferenças, d, tornam-se estacionárias, são chamadas não-estacionárias homogêneas, ou ainda, são portadoras de raízes unitárias. Outras séries não-estacionárias, não-explosivas, são aquelas apresentando uma tendência determinística, como em

X t  0  1t   t , onde  t

RB  0,  2  , que é um processo “trend-stationary”. Então tem-se:

22

(i)

E  X t   t  0  1t ;

(ii)

Tomando-se uma diferença,

X t  X t 1  1   t   t 1 , que é um modelo ARMA (1,1), com     1 , portanto tem-se um modelo não-estacionário e não-invertível. (iii)

Se Wt  X t  X t 1  1  B  X t  X t ,

Wt  X t  1   t , que é um modelo MA(1), estacionário, mas não-invertível (iv)

Extraindo-se a tendência de X t obtém-se

Yt  X t  1t  0   t , que é estacionário. 3.7

MODELOS ARIMA Se Wt   d Zt for estacionária, pode-se representar Wt por um modelo ARMA (p,q), ou

seja,

 ( B)Wt   ( B)at .

(3.35)

Se Wt for uma diferença de Z t , então Z t é uma integral de Wt , pode se dizer que Z t segue um modelo auto-regressivo, integrado, de médias móveis, ou modelo ARIMA,

 ( B)d Zt   ( B)at

(3.36)

de ordem (p,d,q) e escreve-se ARIMA(p,d,q), se p e q são as ordens de   B  e   B  , respectivamente. Assim, o modelo (3.36) supõe que a d-ésima diferença da série Z t pode ser representada por um modelo ARMA, estacionário e invertível. Na maioria dos casos usuais, d=1 ou d=2, que correspondem a dois casos interessantes e comuns de não-estacionaridade homogênea: (a) séries não-estacionárias quanto ao nível; oscilam ao redor de um nível médio durante algum tempo e depois saltam para outro nível temporário. Para tornálas estacionárias é suficiente tomar uma diferença; (b) séries não-estacionárias quanto à inclinação: oscilam numa direção por algum tempo e depois mudam para outra direção temporária. Para torná-las estacionárias é necessário tomar a segunda diferença.

23

O modelo ARIMA (p,d,q) é um caso especial de um processo integrado. Em geral, diz-se que X t é integrado de ordem d se  d X t for estacionário, e escreve-se X t

I  d  . Se

d  0 , X t é estacionário. Nota-se que os modelos ARIMA englobam todos os outros vistos anteriormente, sendo este o mais genérico. Pode-se muito facilmente chegar aos modelos anteriores, bastando para isso zerar uma ou mais ordens do ARIMA:

ARIMA  p,0,0   AR  p  ; ARIMA  0,0, q   MA  q  ;

ARIMA  p,0, q   ARMA  p, q  . 3.8

IDENTIFICAÇÃO DE MODELOS ARIMA Especificada a classe geral dos modelos ARIMA, os estágios do ciclo iterativo do

método de Box e Jenkins são a identificação, a estimação e a verificação. Nesta parte do trabalho será tratada a fase mais crítica do método, que é a identificação do particular modelo ARIMA a ser ajustado aos dados. O objetivo da identificação é determinar os valores p, d e q do modelo ARIMA(p,d,q), e ainda estimativas preliminares dos parâmetros a serem usados no estágio de estimação, de acordo com [7]. O procedimento de identificação consiste de três partes: (a) verificar se a série original tem necessidade de sofrer uma transformação, com o objetivo de estabilizar sua variância. Tal identificação pode ser realizada com o auxílio de gráficos, como visto na seção 2.2; (b) tomar diferenças da série, obtida no item anterior, quantas vezes se fizer necessário para se obter uma série estacionária, de modo que o processo  d Zt seja reduzido a um ARMA(p,q). O número de diferenças, d, necessárias para que o processo se torne estacionário, é alcançado quando a fac amostral de Wt   d Zt decresce rapidamente para zero.

(c) identificar o processo ARMA (p,q) resultante, através da análise das autocorrelações e autocorrelações parciais estimadas, cujos comportamentos devem imitar os comportamentos das respectivas quantidades teóricas.

24

3.9

MÉTODOS BASEADOS EM UMA FUNÇÃO PENALIZADORA A identificação do modelo a ser usado, na construção de modelos ARIMA, é um dos

maiores empecilhos à utilização da metodologia disposta até aqui. Vários pesquisadores, utilizando a mesma série, podem identificar modelos diferentes. Outras propostas de identificação têm sido apresentadas na literatura. Além do método gráfico, outro procedimento para identificação do modelo ARMA é baseado na escolha das ordens k e l que minimizem a quantidade P(k , l )  ln ˆ k2,l  (k  l )

C(N ) N

(3.37)

em que ˆ k2,l é uma estimativa da variância residual obtida ajustando um modelo ARMA(k,l) às N observações da série e C(N) é uma função do tamanho da série. Deste modo o processo de identificação do modelo se faz de forma autônoma, ou seja, sem a necessidade de intervenção humana, o que o torna vantajoso em relação ao método gráfico. A quantidade (k  l )

C(N ) , denominada termo penalizador, aumenta quando o número N

de parâmetros aumenta, enquanto que a variância residual ˆ k2,l diminui. Assim, minimizar (3.37) corresponde a identificar as ordens k e l que equilibrem esse comportamento. Intuitivamente pode ser notado que as ordens selecionadas aumentem quando N cresce. Para tal [13] sugere limites superiores dados por  ln N  , 0     . 

Um resultado bastante importante de [7] será visto a seguir quando será utilizado um critério tipo função penalizadora para identificar os parâmetros p e q de um modelo ARMA(p,q). Teorema: Considere um modelo ARMA(p,q) dado pela expressão (3.26). Sejam pˆ e qˆ os valores que minimizam a expressão (3.37), onde C(N) satisfaz as seguintes condições: lim C ( N )   , lim

N 

N 

C(N )  0. N

Então, pˆ e qˆ são fracamente consistentes. Dois procedimentos de identificação que minimizam funções penalizadoras particulares serão introduzidos adiante. No estudo de caso apresentado neste trabalho eles servirão como alternativa à estimação dos parâmetros p e q dos modelos ARMA.

25

3.9.1

Critério de informação de Akaike

[14] e [15], critério de informação Akaike (AIC), sugere escolher o melhor modelo em que as ordens k e l minimizam o critério AIC(k , d , l )  N ln ˆ a2 

N 2(k  l  1   d 0 )  N ln 2  N , N d

(3.38)

em que 1, d  0 0, d  0

d 0  

e ˆ a2 é o estimador de máxima verossimilhança de  a2 . Para a comparação de vários modelos, com N fixado, os dois últimos termos de (3.38) podem ser abandonados. Levando-se em conta que, geralmente, é identificada a série apropriadamente diferençada, obtém-se AIC (k , l )  N ln ˆ a2  2(k  l  2)

(3.39)

como critério para a determinação das ordens p e q. O que se faz, então, é estipular valores limites superiores K e L para k e l e calcular (3.39) para todas as possíveis combinações (k,l) com 0  k  K e 0  l  L . Em geral, K e L são funções de N, por exemplo K  L  ln N . Pode-se reescrever o AIC da seguinte forma AIC(k , l )  ln ˆ k2,l 

2(k  l ) N

(3.40)

pois os valores de k e l que minimizam esta última expressão são os mesmos que minimizam (3.39). Dependendo dos valores de K e L, muitos modelos têm que ser ajustados, a fim de se obter o mínimo de AIC. Para o caso de modelos AR(p), o critério AIC reduz-se a

AIC (k , l )  N ln ˆ k2,l  2k , k  K .

(3.41)

Existem várias correções para melhorar o comportamento do AIC, no sentido de diminuir a probabilidade de selecionar uma ordem maior do que a verdadeira. [16] propõem uma correção para o AIC dada por AICc (k )  AIC(k ) 

2(k  1)(k  2) ,k  K N k 2

e utilizando simulações, mostram que esta correção é útil quando N é pequeno ou quando K é uma fração “moderadamente grande” de N.

26

3.9.2

Critério de informação Bayesiano

[17], [18] e [19] sugerem minimizar o Critério de Informação Bayesiano (BIC), que no caso de um modelo ARMA é dado por BIC(k , l )  ln ˆ k2,l  (k  l )

ln N , N

(3.42)

em que ˆ k2,l é a estimativa de máxima verossimilhança da variância residual do modelo ARMA (k,l). [20] e [13] mostra que, sob determinadas condições, as estimativas p e q que minimizam a expressão (3.42) são fortemente consistentes. [21] considera o comportamento do BIC quando um limite superior (K,I) para as ordens do modelo ARMA não são especificados. Ele mostra que se for ajustado somente modelos ARMA(r,r), r  0,1, 2,

, utilizando máxima verossimilhança, então o valor de r

que minimiza o BIC é um estimador fortemente consistente de r0  max  p, q  . 3.10 ESTIMAÇÃO DE MODELOS ARIMA Após a identificação de um modelo provisório para a série temporal, o procedimento a seguir é a estimação de seus parâmetros. Em determinado momento far-se-á necessário utilizar um procedimento iterativo de estimação não linear de mínimos quadrados e as estimativas preliminares encontradas na fase de identificação serão usadas como valores iniciais neste procedimento. No modelo SARIMA (p,d,q) considerado será colocado seus p  q  1 parâmetros no vetor   ( , ,  a2 ) , onde   (1 ,...,  p ) ,   (1 ,..., p ) . Assim, se d > 0, supõem-se   0. Senão,  é incluído como mais um parâmetro a ser estimado e irá se ter p  q  2 parâmetros. Seja   ( , ) . Para a estimação de  um método empregado será o da máxima verossimilhança. Este consiste em, dadas as N observações Z1 , L  Z1 ,

Z N , considerar a função de verossimilhança

, Z N  encarada como função de  . Os estimadores de máxima verossimilhança

(EMV) de  serão os valores que maximizam L ou

 log L .

Supondo que o processo at é normal, ou seja, para cada t , at serão aproximadamente estimadores de mínimos quadrados (EMQ).

N  0,  a2  , os EMV

27

A série deverá alcançar a estacionaridade, para tal toma-se d diferenças ficando com

n  N  d observações W1 ,

,Wn , onde Wt   d Zt . Como o modelo ARMA(p,q) resultante é

estacionário e invertível, o processo pode ser escrito como

at  Wt  1Wt 1 

  pWt  p  1at 1 

 q at q

(3.43)

onde W  Wt   . Este trabalho utilizará o método da máxima verossimilhança. Neste, para calcular at , através de (3.43), é necessário obter valores iniciais para os W ' s e para os a ' s . Estes valores iniciais podem ser obtidos através de dois procedimentos: um, condicional, em que os valores iniciais desconhecidos são substituídos por valores que supõem-se serem razoáveis; outro, incondicional, em que os valores iniciais são estimados utilizando um procedimento denominado “backforecasting”.

4

REDES NEURAIS ARTIFICIAIS O comportamento da carga no tempo, a chamada curva de carga, não segue um padrão

claro e é, geralmente, não-linear. Podem variar de acordo com a hora do dia, a temperatura, o dia da semana e muitos outros fatores. Métodos lineares, portanto, não obtêm êxito na modelagem do problema. Modelos neurais tem tido sucesso nessa área, como a literatura tem mostrado. A capacidade de aproximação desse tipo de modelo e a sua característica autônoma o tornam um grande candidato à resolução do problema de previsão. O elevado esforço computacional, porém, há algum tempo atrás, era um importante fator negativo à sua utilização. Com o rápido avanço tecnológico, e consequente alto ritmo de crescimento da capacidade de processamento, aliado a algoritmos de treinamento mais velozes possibilitou a utilização de redes neurais como alternativa viável e promissora. As redes neurais artificiais foram inspiradas na maneira como o cérebro humano processa as informações, diferenciando-se de como um computador digital convencional o faz. Ela pode ser definida, segundo [22], como um processador maciçamente paralelamente distribuído constituído de unidades de processamento simples, que têm a propensão natural para armazenar conhecimento experimental e torná-lo disponível para o uso. Estas unidades de processamento simples são chamadas de neurônios. Um neurônio artificial é esquematicamente representado na Figura 5 como em [23]. O neurônio recebe informação (numérica) através de um número de nós de entrada (quatro, Figura 5), processa estes internamente e devolve uma resposta de saída. O processo é executado normalmente em dois estágios: primeiro, os valores de entrada são linearmente combinados, então o resultado é usado como o argumento da função de ativação não-linear.

29

Figura 5 – Um Neurônio Artificial

Esta estrutura tem a habilidade de desenvolver suas próprias regras através do que usualmente é denominada “experiência”. A disposição desses neurônios define a arquitetura da rede e as conexões entre eles têm um “peso”, chamado de sinapses ou pesos sinápticos. A combinação usa os pesos wi atribuídos a cada conexão e um termo constante “bias”. O procedimento utilizado para realizar o processo de aprendizagem é chamado de algoritmo de aprendizagem, cuja função é modificar os pesos sinápticos da rede de uma forma ordenada para alcançar um objetivo de projeto desejado. A saída Y é dada pela equação, como pode ser visto em [1]:

 n  Y     i xi  b   i 1 

(4.1)

Nesta equação, Y representa a saída do neurônio,   conjunto de pesos sinápticos que ligam as entradas x 

n

n

,   1 ... n  , o t

, x   x1 ... xn  , deste neurônio, t

que podem vir da saída de outro neurônio ou da própria camada de entrada, b  associado e   :



o bias

a sua função de ativação. Aquele que será usado neste trabalho é do

tipo perceptron de múltipla camada ou, em inglês, multilayer perceptron (MLP), no qual os neurônios são organizados em camadas. Os neurônios em cada camada podem compartilhar as mesmas entradas, mas não são conectados uns aos outros. Basicamente, existem duas

30

estruturas, as redes alimentadas adiante, feedforward e as recorrentes. Neste trabalho trabalhará apenas com o modelo feedforward. Neste tipo as saídas de uma camada são utilizadas nas entradas da camada seguinte. As camadas entre os nós de entrada e a camada de saída são chamadas de camadas escondidas.

Figura 6 – Uma Rede Neural Alimentada Adiante com dupla camada A Figura 6 mostra o exemplo de uma rede com quatro nós de entrada, duas camadas, sendo uma escondida, e dois neurônios de saída. Os parâmetros dessa rede são a matriz de pesos W34 (contendo o peso wi , j que conecta o neurônio i à entrada j ), a matriz de pesos

U 23 , e o vetor bias 51 (as conexões do bias não foram representadas na figura). O uso de redes neurais oferece, de acordo com [1], as seguintes propriedades úteis e capacidades: 

Não-linearidade. Um neurônio artificial pode ser linear ou não-linear. Uma rede neural constituída por conexões de neurônios não-lineares é ela mesma não-linear. A não-linearidade é muito importante nos casos em que o sinal de entrada ser inerentemente não-linear. Entretanto, esta vantajosa característica pode ser prejudicial na presença de dados ruidosos.



Mapeamento de Entrada-Saída. Apresenta-se para a rede um conjunto de pares entrada-saída e os pesos sinápticos (parâmetros livres) da rede são modificados

31

para minimizar a diferença entra a resposta desejada e a resposta real da rede. O treinamento é repetido até que não haja modificações significantes nos pesos. 

Adaptabilidade. Possuem elevada capacidade de adaptar seus pesos sinápticos a modificações das condições do ambiente em que foi originalmente treinada, bastando para isso “retreinar” a rede para as novas condições operativas do ambiente.



Implementação Simples. Já que os modelos feedforward podem ser vistos como grafos orientados, a implementação destas estruturas é extremamente simples quando comparada com o grau de complexidade dos modelos que podem ser gerados.

4.1

MODELO DE UM NEURÔNIO O diagrama de blocos da Figura 7 mostra o modelo de um neurônio, que forma a base

para o projeto de redes neurais (artificiais).

Figura 7 – Modelo não-linear de um neurônio. (retirado de [22]) Em três elementos básicos pode ser dividido o modelo neuronal: 1. Um conjunto de sinapses, cada uma caracterizada por um peso independente. Especificamente, um sinal x j na entrada da sinapse j conectada ao neurônio

k é multiplicado pelo peso sináptico wkj . Os índices do pesos sinápticos se referem: o primeiro ao neurônio e o segundo à entrada da sinapse;

32

2. Um somador para somar os sinais de entrada, ponderados pelas respectiva sinapse, no modelo em questão um combinador linear; 3. Uma função de ativação para restringir a amplitude da saída do neurônio. Restrição que manterá a amplitude do sinal de saída em um valor finito. Tipicamente, o intervalo normalizado da amplitude da saída de um neurônio é escrito como o intervalo unitário fechado [0,1] ou, alternativamente [-1,1]. O bias aplicado externamente, denotado por bk , na representação do modelo neural da Figura 7, tem o efeito de aumentar ou diminuir a entrada líquida da função de ativação, dependendo dele ser positivo ou negativo, respectivamente. Pode ser descrito matematicamente um neurônio k com o seguinte par de equações: m

uk   wkj x j

(4.2)

yk    uk  bk 

(4.3)

j 1

e

onde:

x1 , x2 ,

, xm são os sinais de entrada;

wk1 , wk 2 ,

, wkm são os pesos sinápticos do neurônio k ;

uk é a saída do combinador linear devido aos sinais de entrada; bk é o bias;



 é a função de ativação;

yk é o sinal de saída do neurônio. O uso do bias bk tem o efeito de aplicar uma transformação afim à saída uk do combinador linear, assim

vk  uk  bk

(4.4)

A dependência do bias bk ser positivo ou negativo, em particular, faz com que a relação entre o potencial de ativação vk ou campo local induzido do neurônio k e a saída do combinador linear uk seja deslocada da forma ilustrada na Figura 8. Pode ser notado que, como resultado dessa transformação afim, o gráfico vk  uk não passa mais pela origem.

33

Figura 8 – Transformação afim produzida pela presença de um bias [24] 4.2

TIPOS DE FUNÇÃO DE ATIVAÇÃO Representada por   v  , a função de ativação define a saída de um neurônio em

termos do potencial de ativação v . Três tipos básicos de função de ativação são identificados: a função de limiar (Threshold), a função linear por partes (Piecewise-linear) e a função sigmóide. 4.2.1

Função de Limiar

Para este tipo de função de ativação, mostrado na Figura 9, tem-se

1 se v  0 0 se v  0

 v  

(4.5)

Na literatura de engenharia, esta forma de função de limiar é normalmente referida como função de Heaviside. Correspondentemente, a saída do neurônio k que emprega esta função de limiar é expressa como 1 se vk  0 yk  v    0 se vk  0

(4.6)

onde vk é o potencial de ativação do neurônio; isto é m

vk   wkj  bk j 1

(4.7)

34

1,2 1

φ(v)

0,8 0,6 0,4 0,2 0 -2,0

-1,5

-1,0

-0,5

0,0

0,5

1,0

1,5

2,0

v Figura 9 – Função de Limiar A literatura se refere a este neurônio como o modelo de McCulloch-Pitts [22]. Neste modelo a saída de neurônio assume o valor 1, se o potencial de ativação do neurônio é nãonegativo, e 0 em caso contrário. Esta definição descreve a propriedade tudo-ou-nada do modelo de McCulloch-Pitts. 4.2.2

Função Linear por Partes (Piecewise-linear)

Para este tipo de função de ativação, visto na Figura 10, tem-se

1,    v    v, 0, 

v   12  12  v  - 12

(4.8)`

v  - 12

onde o fator de ampliação dentro da região linear de operação é assumido unitário. Esta forma de função de ativação pode ser vista como uma aproximação de um amplificador não-linear. As duas situações a seguir podem ser vistas como formas especiais da função linear por partes: 

Se a saturação não é atingida mantendo-se a região linear de operação, aparece um combinador linear.



Se o fator de amplificação da região linear é feito infinitamente grande, a função de limiar é reduzida à função linear por partes.

35

1,2 1

φ(v)

0,8 0,6 0,4 0,2 0 -2,0

-1,5

-1,0

-0,5

0,0

0,5

1,0

1,5

2,0

v Figura 10 – Função linear por partes (Piecewise-linear) 4.2.3

Função Sigmóide

A função sigmóide, que possui gráfico em forma de S, é sem dúvida a forma mais comum de função de ativação utilizada na construção de redes neurais artificiais. De acordo com [22], ela é definida como uma função estritamente crescente que exibe um balanceamento adequado entre comportamento linear e não-linear. Um exemplo de função sigmóide é a função logística, definida por

 v 

1 1  exp  av 

(4.9)

onde a é o parâmetro de inclinação da função sigmoidal. Com a variação do parâmetro a obtém-se funções sigmóides com diferentes inclinações, como pode ser visto na Figura 11. Na origem a curva obedece à inclinação de a . Em um caso estremo, quando o parâmetro de 4 inclinação se aproxima do infinito, a função sigmóide se torna simplesmente uma função de limiar. Enquanto a função de limiar assume o valor de 0 ou 1, uma função sigmóide assume um intervalo contínuo de valores entre 0 e 1. Deve ser notado também que a função sigmóide é diferenciável, enquanto que a função de limiar não o é. As funções de ativação de limiar, linear por partes e sigmóide, definidas nas equações (4.5), (4.8) e (4.9) se estendem de 0 a +1. Algumas vezes é desejável que a função de ativação cubra o intervalo de -1 a +1, assumindo, nesta situação, uma forma anti-simétrica em relação a origem, ou seja, a função de ativação é uma função ímpar do potencial de ativação.

36

1,2 1

φ(v)

0,8 0,6 a = 0,5 a = 1,0 a = 1,5

0,4 0,2 0 -10,00 -8,00

-6,00

-4,00

-2,00

0,00

2,00

4,00

6,00

8,00

10,00

v Figura 11 – Função Sigmóide para parâmetro de inclinação a variável Especificamente, a função de limiar da equação (4.5) é definida agora como

 1 se v  0    v   0 se v  0 1 se v  0 

(4.10)

esta é denominada função sinal. É desejável que a função sigmóide assuma valores negativos. O fato disso não ser possível com a equação original é uma motivação para se procurar uma forma correspondente, assim, a função tangente hiperbólica satisfaz essa necessidade, e é definida por

  v   tanh  v 

(4.11)

esta nova função de ativação traz benefícios analíticos importantes. 4.3

ARQUITETURAS DE REDE Os algoritmos de aprendizagem utilizados no projeto de redes neurais são classificados

como estruturados devido ao fato dos neurônios que a compõem estarem também estruturados e estes estarão intimamente ligados ao algoritmo. 4.3.1

Redes Alimentadas Adiante (feed forward) com Camada Única

Em uma rede neural em camadas, os neurônios estão dispostos de modo a formar camadas. A forma mais simples de uma rede em camadas consiste em uma camada de entrada de nós de fonte que se projeta sobre uma camada de saída de neurônios (nós computacionais) , mas nunca ao contrário, ou seja, não há realimentação, o percurso será sempre da entrada para

37

a saída. Nesta organização a rede é estritamente do tipo feed forward ou acíclica. A Figura 12 traz um caso com quatro nós tanto na camada de entrada como na de saída. Esta rede é conhecida como rede de camada única, e é designada assim pois, apesar de existir uma camada de entrada de nós de fonte, nesta não é realizada qualquer operação portanto não é contada, sendo a designação “camada única” referente apenas à camada de saída de nós computacionais (neurônios).

Figura 12 – Rede Feed Forward ou acíclica com única camada de neurônios 4.3.2

Redes Feed Forward com Múltiplas Camadas

A segunda classe de rede neural feed forward tem como diferencial a presença de uma ou mais camadas oculta, cujos nós computacionais são chamados correspondentemente de neurônios ocultos ou unidades ocultas. A função dos neurônios ocultos, de acordo com [22], é intervir entre a entrada externa e a saída da rede de forma útil. Adicionar uma ou mais camadas ocultas, torna a rede capaz de extrair estatísticas de ordem elevada. Em um sentido bastante livre, a rede adquire uma perspectiva global apesar de sua conectividade local, devido ao conjunto de conexões sinápticas e da dimensão extra de interações neurais [25]. A habilidade de os neurônios ocultos extraírem estatísticas de ordem elevada é particularmente valiosa quando o tamanho da camada de entrada é grande.

38

Os nós de fonte da camada de entrada da rede provêem os respectivos elementos do padrão de ativação (vetor de entrada), que formam os sinais de entrada aplicados aos neurônios na segunda camada (ou primeira camada oculta). Os sinais da segunda camada servem de entrada para a terceira camada, e esse processo se repete para o resto da rede até que se chega à camada de saída. A resposta global da rede, para o padrão de ativação fornecido pelos nós de fonte da camada de entrada, é constituída pelo conjunto de sinais de saída dos neurônios da camada final. O grafo arquitetural na Figura 13 ilustra a planta de uma rede neural de múltiplas camadas feed forward com uma única camada oculta. Por concisão, a rede na Figura 13 é chamada de rede 9-4-2 pois ela apresenta 9 neurônios de fonte, 4 neurônios ocultos e 2 neurônios de saída. Nos casos em que a rede apresentar mais neurônios na camada oculta sua representação será de maneira análoga, bastando para isso acrescentar o número de neurônios desta nova camada em seu respectivo lugar.

Figura 13 – Rede Feed Forward ou acíclica com uma única camada de neurônios Nos casos em que cada um dos nós de uma camada da rede está conectado a todos os nós da camada adjacente seguinte, a rede neural é dita totalmente conectada. Entretanto, se alguma conexão sináptica estiver faltando na rede, diz-se que a rede é parcialmente conectada.

39

4.3.3

Redes Recorrentes

De [22], uma rede neural recorrente se distingue de uma rede neural feed forward por ter, pelo menos, um laço de realimentação (feedback loop). Uma rede recorrente pode consistir, por exemplo, de uma única camada de neurônios, portanto sem neurônios ocultos, com cada neurônio alimentando seu sinal de saída de volta para as entradas de todos os outros neurônios, como mostrado no grafo arquitetural da Figura 14. Nesta estrutura representada não há laços de auto-realimentação (entende-se por auto-realimentação a situação em que a saída do neurônio é realimentada para sua própria entrada).

Figura 14 – Rede recorrente sem laços de auto-realimentação e sem neurônios ocultos Uma outra classe de redes recorrentes, com neurônios ocultos, é ilustrada na Figura 15. As conexões de realimentação mostradas se originam dos neurônios ocultos assim como dos neurônios de saída. A presença de laços de realimentação, quer seja na estrutura recorrente da Figura 14 ou naquela da Figura 15, tem uma influência significativa na capacidade de aprendizagem da rede e no seu desempenho. Além disso, os laços de realimentação envolvem o uso de ramos particulares compostos de elementos de atraso unitário (representado por z 1 ), o que resulta em um comportamento dinâmico não linear, admitindo-se que a rede neural contenha unidades não-lineares.

40

Figura 15 – Rede recorrente com neurônios ocultos 4.4

PROCESSOS DE APRENDIZAGEM A propriedade fundamental para uma rede neural é, sem sombra de dúvida, sua

inerente capacidade de aprender a partir do seu ambiente e aprimorar seu desempenho através da aprendizagem. A melhoria do desempenho acontece pela seguida utilização da rede e de acordo com algum parâmetro preestabelecido. Uma rede neural aprende sobre seu ambiente através de um processo interativo de ajustes aplicados a seus pesos sinápticos e níveis de bias. Idealmente, a rede se torna mais instruída a cada iteração do processo de aprendizagem. De acordo com [22], o interesse que a rede neural sugere à “aprendizagem” pode ser obtida adaptando a definição dada por [26]: Aprendizagem é um processo pelo qual os parâmetros livres de uma rede neural são adaptados através de um processo de estimulação pelo ambiente no qual a rede está inserida. O tipo de aprendizagem é determinado pela qual a modificação dos parâmetros ocorre. Esta definição implica na seguinte sequência de eventos: 1. A rede neural é estimulada por um ambiente.

41

2. A rede neural sofre modificações nos seus parâmetros livres como resultado desta estimulação. 3. A rede neural responde de uma maneira nova ao ambiente, devido às modificações ocorridas na sua estrutura interna. Um conjunto de regras previamente definidas para a solução de um problema de aprendizagem é denominado um algoritmo de aprendizagem. O projeto de redes neurais não se restringe a um único algoritmo, ele possui um conjunto de algoritmos de aprendizagem que desempenham o papel de diversificar as “ferramentas”, cada uma trazendo, na sua especificidade, vantagem na solução do problema. Os algoritmos de aprendizagem se diferem uns dos outros pela maneira como são ajustados os pesos sinápticos. Basicamente existem cinco regras de aprendizagem: aprendizagem por correção de erro, aprendizagem baseada em memória, aprendizagem hebbiana, aprendizagem competitiva e aprendizagem de Boltzmann. Os paradigmas de aprendizagem cuidam do problema de atribuição de crédito. Este é um procedimento básico para o processo de aprendizagem. Dois destes paradigmas serão tratados a seguir: a aprendizagem com um professor e a aprendizagem sem um professor. 4.4.1

Aprendizagem com um Professor

A aprendizagem com um professor, ou aprendizagem supervisionada, é aquela em que o ambiente é desconhecido pela rede neural, apenas o professor detém o conhecimento sobre o ambiente, e este é inserido através de um conjunto de exemplos de entrada-saída. A Figura 16 mostra um diagrama de blocos que ilustra esta forma de aprendizagem. No processo de aprendizagem o professor e rede neural são expostos a um vetor de treinamento retirado do ambiente. Devido ao fato do professor ter um conhecimento prévio à resposta desejada, ele é capaz de fornecer à rede uma resposta esperada para aquele vetor de treinamento. O sinal do erro é definido como a diferença entre a resposta desejada e a resposta real da rede. Os parâmetros da rede são ajustados iterativamente através da combinação do vetor de treinamento e do sinal de erro, com o objetivo de a rede neural simular a intervenção do professor. Assim, o conhecimento que o professor possui do ambiente, é transferido à rede neural através de treinamento, da forma mais conveniente possível. Quando esta condição é alcançada, pode ser então dispensada a intervenção externa e deixar a rede neural lidar com o ambiente inteiramente autônoma.

42

Figura 16 – Diagrama em blocos da aprendizagem com um professor A forma de aprendizagem supervisionada descrita é por correção de erro e é um sistema realimentado de laço fechado, porém o ambiente desconhecido não está no laço. Pode ser tomado o erro médio quadrado ou a soma dos erros quadrados sobre a amostra de treinamento, como uma medida de desempenho para o sistema, e defini-lo como uma função dos parâmetros livres do sistema. Esta função pode ser mostrada como uma superfície multidimensional de desempenho de erro com os parâmetros livres como coordenadas. Esta superfície de erro é, na realidade, obtida pela média sobre todos os exemplos possíveis de entrada-saída como descrito em [22]. Qualquer operação de sistema sob supervisão do professor é representada como um ponto sobre a superfície de erro. Para que o sistema melhore o seu desempenho ao longo do tempo e, portanto, aprenda com o professor, o ponto de operação deve ser movido para baixo sucessivamente em direção a um ponto mínimo da superfície de erro; o ponto mínimo pode ser um mínimo local ou um mínimo global . Um sistema de aprendizagem supervisionada é capaz de fazer isto com a informação útil que ele tem sobre o gradiente da superfície de erro, correspondente ao comportamento corrente do sistema. O gradiente de uma superfície de erro em qualquer ponto é um vetor que aponta na direção da descida mais íngreme. Na verdade, no caso da aprendizagem supervisionada por exemplos, o sistema pode usar a estimativa instantânea do vetor gradiente, supondo que os índices dos exemplos sejam os mesmos dos instantes de tempo. O uso de tal estimativa resulta

43

em um movimento do ponto de operação sobre a superfície de erro que se dá tipicamente na forma de uma “caminhada aleatória”. Apesar disso, dados um algoritmo projetado para minimizar a função de custo, um conjunto adequado de exemplos de entrada-saída e tempo suficiente para realizar o treinamento, um sistema de aprendizagem supervisionada é normalmente capaz de realizar tarefas como classificação de padrões e aproximação de funções. 4.4.2

Aprendizagem sem um Professor

Como já era esperado, devido a implicação do nome, no paradigma conhecido como aprendizagem sem um professor, não há um professor para supervisionar o processo de aprendizagem. Neste paradigma não há exemplos rotulados da função a ser aprendida pela rede e pode ser subdividido em dois, a aprendizagem por reforço/Programação neurodinâmica e a aprendizagem não supervisionada. Este trabalho tratará apenas do segundo. 4.4.2.1

Aprendizagem não-supervisionada

Na aprendizagem não-supervisionada ou auto-organizada, não há um professor externo ou um crítico para supervisionar o processo de aprendizado, como indicado na Figura 17. Em vez disso são dadas condições para realizar uma medida independente da tarefa da qualidade da representação que a rede deve aprender, e os parâmetros livres da rede são otimizados em relação a esta medida. Uma vez que a rede tenha se ajustado às regularidades estatísticas dos dados de entrada, ela desenvolve a habilidade de formar representações internas para codificar as características da entrada e, desse modo, de criar automaticamente novas classes [27].

Figura 17 – Diagrama de blocos da aprendizagem não supervisionada Para realização da aprendizagem não-supervisionada pode ser utilizada a regra da aprendizagem competitiva. Em uma rede neural de duas camadas por exemplo, com uma camada de entrada e uma competitiva, na camada de entrada ocorre a recepção dos dados disponíveis e na camada competitiva os neurônios competem, de acordo com as regras de

44

aprendizagem, pela “oportunidade” de responder às informações oriundas dos dados de entrada. Em outras palavras, a rede funciona de acordo com uma estratégia do tipo “o vencedor leva tudo”, o neurônio de maior entrada total “ganha” a disputa e se torna ligado, todos os outros neurônios se tornam desligados. 4.5

ETAPAS GERAIS DE PROJETOS DE REDES NEURAIS De acordo com [28], um projeto de redes neurais pode ser dividido em etapas. A

Figura 18 traz uma rotina a ser seguida visando a implantação de um projeto.

Figura 18 – Diagrama em Blocos das Etapas Gerais de um Projeto de Redes Neurais Cada uma das caixas acima exerce papel fundamental na elaboração de um projeto. A definição do espaço de entrada escolhe as informações a serem utilizadas e representa as variáveis em um formato conveniente (no caso das variáveis de entradas estarem em um formato indesejado). No pré-processamento ocorre a normalização das entradas e saídas e a extração e seleção das características. Na escolha da estrutura se define, pelo teorema da aproximação universal, o número de neurônios por camada, as máquinas de vetor suporte e os critérios para escolha. No treinamento do modelo escolhe-se o algoritmo a ser utilizado e neste os critérios de parada. A avaliação do modelo faz uma validação cruzada e são tomadas medidas analíticas de avaliação do modelo. Caso o modelo escolhido cumpra com os critérios de avaliação, ele estará apto a ser utilizado, caso contrário escolhe-se novamente uma estrutura e o resto do processo se repete. Serão detalhadas a seguir cada um dessas etapas.

45

4.5.1

Definição do espaço de entrada

A definição do espaço de entrada deve ser de acordo com as informações disponíveis. As variáveis numéricas estão de acordo com o esperado da entrada e as variáveis nãonuméricas recebem uma representação numérica, como, por exemplo, onde se deseja definir a intensidade da tonalidade de um pixel, pode-se usar a seguinte correspondência:

Pixel branco  0 Pixel cinza  0,5 Pixel preto  1 São designados índices para a mensuração do relacionamento entre os dados de entrada. Quando é encontrada uma correlação linear, este padrão pode facilmente ser percebido em processos computacionais com a aplicação dos métodos clássicos de identificação linear de sistema, não será necessária a utilização de informação redundante. Na literatura de modelos neurais, os métodos de seleção de variáveis de entrada são divididos em dois grupo, a saber: as técnicas de filtragem e as encapsuladas. Para problemas de previsão a última é mais indicada, como pode ser visto em [1], porém com a desvantagem de necessidade de intervenção. Assim, para uma tomada de decisão automática, técnicas empíricas foram necessárias. Os modelos para previsão encontrados adiante nesse trabalho utilizaram técnicas vindas da teoria do caos. 4.5.2

Pré-processamento

Nesta fase acontecem transformações aplicadas às variáveis de entrada. As transformações ocorrem respeitando as características das funções de ativação definidas para a rede neural buscando uniformizar características estatísticas dos dados. Identificadas as características das variáveis de entrada, são retiradas as que apresentam informação redundante, fato louvável, pois é reduzida a dimensão do vetor de entrada para cada padrão entrada-saída. As transformações são úteis pelo fato de trazerem as informações, tanto das variáveis de entrada quanto das de saída, para a região de operação dos neurônios de modo que estes não operem na sua região de saturação. Normalmente é feita uma transformação linear mudando a faixa de operação dos sinais para ficarem compreendidos entre 0 e 1 [0,1], e este fato não altera comportamento da variável. Uma transformação logarítmica é também utilizada, porém com menor frequência. O processo de extração das características com consequente redução da dimensão pode ser feita através da análise de componentes principais (Principal Component Analysis – PCA).

46

Um problema comum em reconhecimento estatístico de padrões, de acordo com [22], é a seleção das características ou extração de características. A seleção de características se refere a um processo no qual um espaço de dados é transformado em um espaço de características que, em teoria, tem exatamente a mesma dimensão do espaço original de dados. Entretanto, a transformação é projetada de tal forma que o conjunto de dados pode ser representado por um número reduzido de características “efetivas” e ainda reter a maioria do conteúdo de informação intrínseco dos dados. Em outras palavras, o conjunto de dados sofre uma redução de dimensionalidade. De forma mais prática, seja um vetor x de dimensão m e deseja-se transmiti-lo usando l números, onde l  m . Se simplesmente o vetor x for trucado, um erro médio quadrado, igual à soma das variâncias dos elementos eliminados de x ,será causado. Uma transformação linear Q , no sentido de otimizar o erro médio quadrado, deve ter a propriedade que alguns de seus componentes tenham baixa variância. A análise de componentes principais, conhecida na teoria da comunicação como transformação de Karhunen-Loève, maximiza a taxa de redução da variância e acelera a convergência do algoritmo de retropropagação de erro (MLPs). Os componentes principais z vêem da equação (4.12) e são ortogonais entre si. z  Qx

(4.12)

A Figura 19, traz a análise de componentes principais. Quando m  n são retiradas as redundâncias e quando m  n ocorre a redução da dimensionalidade. Na filtragem das saídas devem ser ainda identificados os dados faltantes e os dados discrepantes (outliers). Este último podendo utilizar a hipótese de distribuição gaussiana.

47

Figura 19 – Análise de Componentes Principais 4.5.3

Escolha da Estrutura

Nesta etapa é escolhida a arquitetura da rede baseando-se no conhecimento e na capacidade de elaboração do “projetista” sobre a estrutura, de acordo com o problema a ser enfrentado. Ele define o número de camadas ocultas, com o número de neurônio por camada, e os parâmetros iniciais. São exemplos de estrutura: perceptrons de múltiplas camadas (MLP); redes recorrentes; modelos baseados em kernel. Na definição dos parâmetros iniciais, os MLP e as redes recorrentes apresentam espaço de busca para número de neurônios, já nos modelos baseados em kernel o espaço de busca é para hiperparâmetros (SVM). O teorema da aproximação universal [29] afirma que os modelos feedforward podem aproximar com precisão arbitrária qualquer função contínua F  x  :

n



. Para tal, a

estrutura deve apresentar ao menos uma camada oculta contendo neurônios com função de ativação contínua, não-constante, limitada, e uma saída linear, representando a aproximação de F  x  gerada pelo modelo. Portanto, modelos feedforward com uma única camada escondida contendo um número suficiente de neurônios com função de ativação com as características anteriormente citadas podem aproximar qualquer função contínua. Esta característica constitui a principal motivação para utilização ao longo deste trabalho de modelos com uma única camada escondida.

48

4.5.4

Treinamento do Modelo

Nos MLPs e nas redes recorrentes a definição do algoritmo de treinamento se dá através da equação (4.13) w  k   w  k  1    k    k     k  w  k  2

(4.13)

onde:  é a taxa de aprendizagem,  o ajuste e  a taxa de momento. A taxa de aprendizagem pode ser interpretada como o passo de atualização do algoritmo, quanto maior esta for, maior será a velocidade da resposta, no entanto um aumento traz consigo uma instabilidade maior. A taxa de momento evita mudanças bruscas ao longo do processo de busca, contribuindo dessa forma para o aumento da estabilidade (decréscimo de instabilidade). Devem ser definidos critérios de parada já que não existe garantia analítica de convergência do modelo. Em [28] são três estes critérios: magnitude do gradiente, variação do erro entre interações e parada antecipada do treinamento. 4.5.5

Avaliação do Modelo

A avaliação do modelo pode ser feita via validação cruzada e através de medidas analíticas de avaliação de modelos. 4.5.5.1

Inferência Bayesiana

A inferência bayesiana pode ser utilizada para seleção da melhor estrutura em uma série de hipóteses   H1 , H 2 ,

, H K  . Pela regra de Bayes, a distribuição de probabilidade

a posteriori p  H h Y  da hipótese H h é dada por: p  H h Y  

p Y H h  p  H h  p Y 

. Visto que

p Y  é um fator de normalização e admitindo que todas as hipóteses são equiprováveis a

priori, a evidência p Y H h  pode ser utilizada para a avaliação de modelos, sendo selecionado aquele com maior probabilidade a posteriori p  H h Y  , ou seja, maior evidência [29]. 4.5.5.2

Validação Cruzada

De acordo com [22] A essência da aprendizagem por retropropagação é codificar um mapeamento de entrada-saída (representado por um conjunto de exemplos rotulados) nos

49

pesos sinápticos e bias de uma MLP. É esperado que a rede se torne bem-treinada de modo que aprenda suficiente o passado para generalizar no futuro. Desta perspectiva, o processo de aprendizagem se transforma em uma escolha da parametrização da rede para este conjunto de dados. Mais especificamente, o problema da seleção da rede pode ser visto como a escolha, dentre um conjunto de estruturas de modelos candidatas (parametrizações), a “melhor” de acordo com um certo critério. Neste contexto, uma ferramenta padrão da estatística conhecida como validação cruzada fornece um princípio orientador atraente [30][31]. Primeiramente, o conjunto de dados disponível é dividido aleatoriamente em um conjunto de treinamento e um conjunto de teste. O conjunto de treinamento é dividido em dois subconjuntos disjuntos: 

Subconjunto de estimação, usado para selecionar o modelo.



Subconjunto de validação, usado para testar ou validar o modelo.

A motivação aqui é validar o modelo com um conjunto de dados diferente daquele usado para estimar os parâmetros. Desta forma, pode-se usar o conjunto de treinamento para avaliar o desempenho de vários modelos candidatos e, assim, escolher o “melhor”. Há, entretanto, uma possibilidade considerável de que o modelo assim selecionado, com os valores de parâmetros com melhor desempenho, possa acabar ajustando excessivamente o subconjunto de validação. Para desviar-se dessa possibilidade, o desempenho da generalização do modelo selecionado é medido sobre o conjunto de teste, que é diferente do subconjunto de validação. O uso de validação cruzada é atrativo particularmente quando se tem que projetar uma rede neural grande cujo objetivo seja uma boa generalização. Pode-se, por exemplo, utilizar a validação cruzada para determinar a MLP com o melhor número de neurônios ocultos e quando é melhor para o treinamento evitar o ajuste excessivo (ou overfitting), que constitui uma das principais tarefas no desenvolvimento de modelos neurais.

5

RESULTADOS Este capítulo apresenta as séries de carga utilizadas, algumas de suas características e

os processamentos efetuados nelas. Será descrita a metodologia para a aplicação de Análise Clássica de séries temporais e de Rede Neural. Antes da descrição das séries, vale ressaltar que neste trabalho serão consideradas previsões de carga para o curto prazo, previsões realizadas para horizontes variando de um dia a até um mês à frente, discretizadas em base diária. Visto que este trabalho busca a comparação entre metodologias e não o desenvolvimento de um sistema completo de previsão de carga para curto prazo, as estruturas a serem utilizadas foram definidas buscando este objetivo. 5.1

DADOS DO HISTÓRICO DA CARGA A base de dados de carga utilizada foi obtida de uma concessionária de energia

européia, a East-Slovakia Power Distribution Company, disponibilizados pelo European Network on Intelligent Technologies for Smart Adaptive Systems, popularmente conhecido pela sigla EUNITE, para realização de uma competição entre modelos de previsão em 2001, e encontrados em [32]. Para esta base de dados, dos anos de 1997 e 1998 divididos em intervalos de 30 minutos. Os modelos desenvolvidos objetivam a previsão da carga máxima diária um passo à frente, com estas previsões sendo utilizadas para previsão das 31 cargas máximas verificadas em janeiro de 1999. As demandas máximas que serviram de base de dados foram colocadas em um gráfico para uma visualização de suas características. A fim de comparar os anos em questão (1997 e 1998) os dois foram sobrepostos, como mostra a Figura 20. A dinâmica de curto prazo da carga é fortemente influenciada pelo período do ano [33], em função das condições climáticas relacionadas com as estações do ano e pelo dia da semana, em virtude da diferença do padrão de consumo entre dias úteis e finais de semana. A influência da época do ano no comportamento da demanda a curto prazo pode ser verificada facilmente. Na Figura

51

20, a diferença do comportamento da carga no inverno e no verão é nítida, com a estação do ano influenciando fortemente o pico de carga diário. Nota-se um “vale” para os dias no centro da curva. 900

Carga Máxima Diária

850 800 750

[MW]

700 650 600 550

1997

500

1998

450

Dias do Ano Figura 20 – Demandas Máximas Diárias de 1997 e 1998 (Dias do Ano x Demanda [MW]) A influência do dia da semana na dinâmica de curto prazo da carga pode ser verificada na Figura 21, observando-se claramente os altos e baixos periódicos (com 7 dias de intervalo). 900

Carga Máxima Diária

850 800

[MW]

750 700

Mês de Janeiro Mês de Julho

650 600 550 500 450

1

8

15

Dias do Mês

22

Figura 21 – Mês de Janeiro e Julho de 1997

29

52

5.2

DEFINIÇÃO DOS MODELOS 5.2.1

Modelos SARIMA

Para a aplicação de análise clássica de séries temporais para a previsão de potência ativa em barramentos, foram usados algoritmos programados direcionados ao Matlab, fazendo uso de toolbox – GARCH (Generalized Autoregressive Conditional heteroskedasticity) que é especifica para aplicação de Análise de Séries Temporais. Adiante será apresentada a rotina para a construção do modelo, Figura 22.

Figura 22 – Rotina de Construção de um Modelo Para a modelagem da série será proposta a seguinte configuração para a série, vista na equação (5.1).

Zt  Tt  St  at

(5.1)

Esta série apresentada sofrerá um processamento para a retirada da tendência ( Tt ) e da sazonalidade ( St ), com a retirada espera-se que a série se constitua apenas da componente aleatória at , que é estacionária, ou espera-se que seja, já que esta é uma característica desejável à analise da série. 5.2.1.1

Retirada da Tendência

A primeira transformação aplicada à série de carga original Z  t  é dada pela equação (5.2). O intuito é a retirada da tendência local de crescimento desta série, fazendo com que agora as estações do ano não tenham influência na previsão da demanda.

53

Z1  t   Z  t   Z  t  1

(5.2)

Foi utilizado o método das diferenças sucessivas. Apenas a primeira diferença foi necessária. 5.2.1.2

Retirada da Sazonalidade

Como pôde ser notado nos dados, a série apresenta uma forte influência do dia da semana na demanda máxima. Sendo o horizonte de previsão de um mês, com discretização diária, deve ser retirada essa sazonalidade semanal. A série processada não será mais a inicial, e sim a que foi gerada após a primeira diferença, que será: Z2  t   Z1  t   Z1  t  7 

(5.3)

Agora a nova série processada ficará sem a tendência Tt e sem a sazonalidade St , restando apenas o sinal do processo estacionário at . 5.2.1.3

Definição do Modelo

Neste trabalham não foram desenvolvidos os modelos de previsão, portanto não serão detalhadas passo a passo as tentativas de identificação do mesmo. Na ferramenta toolbox – GARCH (Generalized Autoregressive Conditional heteroskedasticity), com o qual o Matlab processa seus cálculos, é necessário especificar a variância, uma vez que o modelo GARCH não é linear e supõe que a variância é variável. Para o caso em que a variância é especificada como constante, o modelo GARCH torna-se um modelo ARMAX, onde o X depende de fatores exógenos, que não estão disponíveis para este estudo. Deve-se, portanto, ser utilizada outra ferramenta, o garchset, que tratará o problema da variância. Para o método gráfico, serão utilizadas as funções de autocorrelação (fac) e autocorrelação parcial (facp), para estimação dos parâmetros p e q do modelo ARMA(p,q) pelo método da máxima verossimilhança. Com a utilização da função garchfit são estimados os parâmetros da série. Serão obtidos portanto os coeficientes do modelo ARMA, 1 , 2 ,

, n

, que correspondem a parcela auto-regressiva (AR), onde n é define a ordem p do modelo,

1 , 2 , , m da parcela da média móvel (MA), com m sendo a ordem q do modelo, a variância (constante) (K) e a média (C). Ainda há uma série denominada resíduo (innovations), que será obtida a partir da série processada Z 2  t  , essa nova série verificará o ajuste do modelo.

54

Outra forma de analisar a ordem do modelo ARMA é pelo critério AIC, ou critério de informação de Akaike (Akaike’s Information Criterion), e BIC, ou critério de informação bayesiano (Bayesian Information Criterion), esse critério verifica a equação penalizadora, de forma que haja um equilíbrio entre o modelo ajustado e o número de coeficientes do modelo. A função aicbic do Matlab, fornece essas ordens, bastando informar a ordem máxima a ser testada. Vale ser ressaltado que quanto maior a ordem maior o esforço computacional para testar todas as combinações p e q possíveis. Em contrapartida, neste método não é necessária uma intervenção externa. Após a definição dos parâmetros, os resíduos gerados, com a aplicação das funções fac e facp, serão analisados e se estiverem fora dos limites esperados serão estimados novamente os parâmetros. Esse processo se repete até que se atinja a condição desejada. No caso em estudo os parâmetros serão, para o método gráfico p  7 e q  7 e pelo método automático p  2 e q  7 . Portanto para o método gráfico e para o método automático foram escolhidos o

ARMA(7,7) e o ARMA(2,7), respectivamente. 5.2.2

Modelos Neurais

Os modelo neurais utilizados para o desenvolvimento deste trabalho trarão características distintas quanto ao tratamento da sazonalidade, à seleção do conjunto inicial de entradas via teoria do caos e com e sem clusterização. Não será retirada a tendência linear da série, nem tampouco serão utilizadas variáveis exógenas. O modelo de previsão será o Perceptron de múltipla camada Bayesiano (ou Bayesian MLP) e o algoritmo de clusterização será o RPEM. A função utilizada para a previsão foi a FACCBNF (Fully Automatic Chaotic Clustering Bayesian Neural Forecaster que seria, na tradução ao pé da letra, o seguinte: Previsão Neural Bayesiana com Clusterização Caótica Totalmente Automática). Serão três os tipos de tratamento da sazonalidade. O primeiro será sem tratamento. O segundo fará o processamento mensal por meio de séries de Fourier. No terceiro, e último, a sazonalidade será processada por um modelo neural utilizando somente binárias como entrada. O preparo dos dados de entrada será tratado pela teoria do caos, e apresentarão lags (atrasos) ou definidos pelo usuário ou automaticamente. A clusterização poderá ser sem agrupamento, ou seja, com um único cluster, ou ser com até três agrupamentos para modelagem do atrator.

55

5.2.2.1

Sem Tratamento para a Sazonalidade

Na série sem tratamento de sazonalidade a rotina seguida para a previsão será mostrada na Figura 23. O modelo recebe os dados originais, monta os padrões de entradasaída sem levar em consideração a sazonalidade, treina a rede e finalmente gera a previsão.

Figura 23 – Rotina da Previsão da Rede Neural sem Sazonalidade 5.2.2.2

Tratamento da Sazonalidade com Análise Espectral via Séries de Fourier

No tratamento via séries de Fourier é utilizada uma função para modelagem determinística de sazonalidades presentes em séries temporais. São utilizados testes de periodicidade para identificação das frequências [7]. A Figura 24 ilustra a rotina. A primeira caixa do modelo – Dados Originais – lê os dados de entrada. A segunda – Análise Espectral – faz o processamento da sazonalidade utilizando a análise espectral via séries de Fourier. Para tal é usada uma função para modelagem determinística de sazonalidades presentes em séries temporais. São feitos testes de periodicidade para identificação das frequências, no caso, será utilizado o teste de Whittle [7]. Na terceira caixa – Cálculo do Resíduo – é retirada a componente sazonal, sobrando o resíduo. Na quarta é montado o padrão de entrada-saída, são dados de entrada a série sem tendência (ou deveria ser, porém nesta simulação não houve retirada da tendência) e o nível de significância. A quinta treina a rede neural. A sexta faz a modelagem da componente sazonal via séries de Fourier. Na sétima caixa tem-se a previsão para o resíduo. Finalmente a última gera a previsão somando o resíduo à previsão sazonal.

56

Figura 24 – Rotina da Previsão com Análise Espectral para modelagem da Sazonalidade 5.2.2.3

Sazonalidade Processada por Rede Neural com Entradas Binárias

O último caso de tratamento da sazonalidade, o processamento desta é por meio de variáveis binárias retira todas as frequências relevantes dos sinais, ao contrário da análise automática via séries de Fourier. Para séries diárias, a inclusão da binária para o número da semana no mês faz com que o sinal resultante não tenha nenhuma frequência. A rotina segue o esquema da Figura 25. A primeira tarefa da rotina será o recebimento dos dados – Dados Originais. Montagem do Padrão Entrada-Saída – a entrada nesse caso será binária e a saída será como dados da carga (formato original). Treinamento da Rede – ocorre o treinamento da rede considerando somente os dados normais. Cálculo do Resíduo – são retirados os dados referentes a outliers (fora do padrão usual) previamente detectados. Como os dados relacionados com outliers serão filtrados (substituídos) pelas estimativas do modelo sazonal, para estes pontos, o valor da série processada e igual a zero. Portanto, o resíduo fica

57

Resíduo = Dados - Filtro(RN) . Montagem do Padrão Entrada-Saída – Na entrada serão só os

atrasos e na saída os resíduos. Como as binárias foram utilizadas para o processamento anterior da sazonalidades, elas devem ser retiradas do modelo neural. Treinamento da Rede – a rede é treinada agora normalmente como nos modelos anteriores. Previsão para o Resíduo – é efetuada a previsão para o resíduo. Previsão para o Filtro – faz-se a previsão para o filtro via rede neural. Previsão (Resíduo + Filtro) – retorna a componente do filtro ao resíduo e gera-se a previsão final.

Figura 25 – Rotina da Previsão por meio de Variáveis Binárias 5.2.2.4

Seleção do Conjunto Inicial de Entrada pela Teoria do Caos

Como dito anteriormente a seleção do conjunto inicial de entrada pela teoria do caos poderá ter seus atrasos (lags) definidos pelo usuário ou então de forma automática pelo primeiro mínimo da função de informação mútua, na função takens_embedding. Nos atrasos definidos pelo usuário foram utilizados lags de um a sete [1 2 3 4 5 6 7].

58

Já na função para determinação dos atrasos de uma dada série temporal a serem utilizados como entradas dos modelos de previsão, a escolha e baseada no teorema de Takens, que afirma que o espaço de estados de um sistema caótico determinístico pode ser reconstruído através do vetor d-dimensional de atrasos da serie temporal, dado por:

x  t    y  t  y  t  T 

y t   d  1 T 

T – atraso da imersão;

Sendo:

d – dimensão da imersão. Para determinação do atraso da imersão, será utilizado o primeiro mínimo da função de auto-informação mútua [34]. Para determinação da dimensão da imersão, será utilizado o método dos falsos vizinhos mais próximos proposto em [35]. 5.3

COMPARAÇÃO DOS RESULTADOS OBTIDOS Para os modelos utilizados para a previsão do pico de carga, utilizando os

procedimentos discutidos até aqui para a seleção dos parâmetros que definem os modelos SARIMA e os modelos neurais, foram obtidas as previsões apresentadas na Tabela 1. Na primeira colunas traz os dias do mês a ser previsto, janeiro de 1999. Na segunda coluna são mostrados os valores medidos em mega-watts [MW], que seriam os resultados ideais a serem atingidos. Na terceira tem-se a primeira previsão realizada, utilizando o método gráfico do modelo SARIMA, neste caso foi um SARIMA(7,7). Na coluna seguinte, quarta, são mostrados os resultados obtidos pelo SARIMA automático, baseado na função penalizadora, AICBIC, resultando o SARIMA(2,7). As colunas seguintes trarão as previsões dos modelos de redes neurais (NN). Os índices mostrados entre parênteses definem os parâmetros usados nas funções de previsão, a seguir será explicado como fazer a sua leitura. O primeiro, (*,*,*), é referente a sazonalidade. 

0 (zero) – série sem tratamento de sazonalidade;



1 (um) – tratamento da sazonalidade via séries de Fourier;



3 (três) – processamento da sazonalidade é por meio de variáveis binárias.

O segundo, (*,*,*), é referente à seleção da entrada pela teoria do caos. 

0 (zero) – lags definidos pelo usuário;



1 (um) – lags escolhidos pela função.

O Terceiro, (*,*,*), é referente a clusterização dos dados. 

0 (zero) – sem clusterização;



1 (um) – com clusterização.

59

Dia do Mês

Valores Medidos [MW]

SARIMA Gráfico

SARIMA AICBIC

NN (0,0,0)

NN (0,0,1)

NN (0,1,0)

NN (0,1,1)

NN (1,0,0)

NN (1,0,1)

NN (1,1,0)

NN (1,1,1)

NN (3,0,0)

NN (3,0,1)

NN (3,1,0)

NN (3,1,1)

Tabela 1 – Previsão dos Modelos em Estudo

1

751

728

700

732

730

735

736

729

730

736

735

730

733

731

737

2

703

708

691

707

706

712

713

704

706

713

712

709

703

716

719

3

677

673

700

666

665

669

670

665

665

670

669

673

675

686

687

4

718

738

736

739

739

743

744

738

739

743

742

741

743

761

761

5

738

740

740

746

746

749

752

744

745

750

749

748

748

769

766

6

709

751

749

749

750

753

755

747

749

752

752

751

756

778

769

7

745

735

730

743

745

748

749

742

743

745

748

746

751

783

766

8

749

729

697

739

739

739

740

735

737

737

740

740

749

779

764

9

734

709

689

712

712

714

715

708

711

711

714

717

723

760

737

10

679

671

699

669

669

668

669

666

668

665

669

679

690

724

698

11

748

736

734

743

743

739

741

740

742

737

741

748

758

794

765

12

739

739

738

749

750

743

746

746

749

742

746

754

764

799

765

13

756

749

748

752

753

747

749

749

752

743

749

757

769

803

765

14

763

733

729

746

749

742

743

744

746

738

744

753

765

801

760

15

752

728

696

742

743

732

734

737

740

729

736

747

762

793

757

16

738

707

688

715

717

705

707

711

714

702

709

724

735

770

733

17

699

669

697

673

674

658

661

669

672

656

664

686

701

733

696

18

782

735

733

747

748

730

733

742

746

728

737

754

769

799

766

19

782

737

737

753

754

735

739

748

752

733

742

761

775

802

770

20

792

748

746

756

757

738

743

751

755

736

746

764

779

802

772

21

801

731

727

750

752

734

737

746

750

731

741

759

774

801

769

22

781

726

695

746

747

725

729

739

744

723

733

753

770

793

763

23

731

705

686

719

720

698

703

713

718

696

707

730

743

770

738

24

708

668

696

676

678

652

658

671

676

651

662

692

708

731

700

25

789

733

731

750

751

725

730

744

750

723

735

760

775

794

768

26

798

735

735

756

758

730

737

750

755

729

741

766

780

794

769

27

791

746

745

759

760

734

740

753

759

731

745

769

784

796

771

28

776

730

726

752

756

729

735

748

753

726

741

764

778

794

767

29

792

724

693

748

750

720

727

741

747

718

733

758

774

784

762

30

763

704

685

722

724

694

701

715

721

691

707

735

746

759

735

31

743

666

694

679

681

650

657

673

679

647

663

697

712

721

696

Métodos de Previsão Estudados (resultados em MW)

Com os valores obtidos podem ser calculados o erro percentual absoluto médio (MAPE – Mean Absolute Percentage Error) e o erro absoluto máximo (MAE – Maximum Absolute Error) dos modelos. Para avaliação dos modelos serão utilizados os mesmos critérios que o EUNITE usou como medida de desempenho, o MAPE e o MAE, anteriormente mencionados. Além destes ainda poderá ser avaliado o esforço computacional,

60

que, para o estudo, foi utilizada a contagem do tempo em segundos como critério para comparação, quanto menor o tempo mais preferível é o modelo, já que para o tipo de previsão tratada neste trabalho, curto prazo, o tempo de resposta é um fator preponderante. O programa foi processado em um computador com um processador Intel Celeron M CPU 530 @ 1.73GHz e 1,99GB de RAM. Os resultados obtidos estão na Tabela 2. Em destaque, os melhores resultados para cada critério.

Tabela 2 – MAPE, MAE e Tempo para Previsão

NN (3,1,1)

NN (3,1,0)

NN (3,0,1)

NN (3,0,0)

NN (1,1,1)

NN (1,1,0)

NN (1,0,1)

NN (1,0,0)

NN (0,1,1)

NN (0,1,0)

NN (0,0,1)

NN (0,0,0)

SARIMA AICBIC

SARIMA Gráfico

Métodos de Previsão Estudados

MAPE 4,39 5,27 3,24 3,11 4,80 4,45 3,68 3,29 5,01 4,12 2,38 1,66 3,54 2,56 [%] MAE 77,2 99,0 63,8 61,8 93,3 86,0 70,4 63,7 95,6 80,0 46,3 47,4 69,5 59,7 [MW] Tempo para 1,66 4,00 7,72 323 10,6 343 9,58 316 11,4 348 13,4 242 13,8 264 Previsão [s]

Entre os modelos SARIMA, o método gráfico possui menores erros e menor esforço computacional, portanto, aparentemente é melhor que o método autônomo, aparente pois, o modelo gráfico necessita de intervenção do usuário. Entre os modelos neurais, o NN(3,0,0) obteve o menor MAE e o NN(3,0,1) o menor (MAPE). Pode ser destacado que, quando processada a sazonalidade através da rede neural com entradas binárias, os resultados são nitidamente melhores que os demais. Quanto à seleção das entradas pela teoria do caos pode ser dito que o que utiliza os atrasos definidos pelo usuário fica a frente da forma automática, onde também não deve esquecido que o usuário defina o número de atrasos. A clusterização, em todos os modelos, melhora o resultado da previsão, mas onera muito o tempo de processamento, devendo ser levado em conta na escolha do modelo. Como não se sabe em que situação será aplicado o modelo de previsão, para comparação dos resultados não será levado em conta o esforço computacional.

61

Em competição promovida pelo EUNITE em 2001, os resultados obtidos pelos dez primeiros colocados, disponibilizados no site da competição [32], estão apresentados na Tabela 3.

Tabela 3 – Resultado dos 10 Primeiros Colocados na competição da EUNITE.

1º 2º 3º 4º 5º 6º 7º 8º 9º 10º

MAPE [%]

MAE [MW]

1,9824 2,1489 2,4979 2,8733 2,9846 3,2234 3,2644 3,3797 3,3880 3,3887

51,4 40,0 60,5 65,8 68,6 55,2 77,0 74,0 109,0 74,5

Entre os modelos estudados neste trabalho, se os melhores, hipoteticamente, estivessem competindo teriam conseguido uma posição de destaque. O NN(3,0,0) ficaria em segundo tanto para o MAPE quanto para o MAE. O NN(3,0,1) obteria a primeira colocação referente ao MAPE e na segunda no quesito MAE (apresenta um tempo de resposta 18 (dezoito) vezes maior que o NN(3,0,0), mas isto não é considerado no concurso). A Figura 26 mostra os gráficos com as previsões destes modelos.

62

810 790

Carga Máxima Diária

770 750

[MW]

730 710 690 670 650

1

8

Valores Medidos [MW] 2º NN (3,0,0)

15

22

29

Dias do Mês

1º SARIMA Gráfico NN (3,0,1)

Figura 26 – Previsões de Janeiro A competição contou com participantes de vários países, ao todo foram 56 participantes de 21 países, a Figura 27 ilustra essa distribuição.

56 Competidores registrados de 21 Países EUA 5 Ucrânia 1 Taiwan 1

Austrália 3

Bélgica 1 Canadá 1 Rep. Tcheca1 Inglaterra 5

Espanha 6

Finlândia 1 França 2

Eslováquia 4

Alemanha 6

Rússia 2 Polônia5 Noruega 1

Holanda 4

Coréia 1

Grécia 1 Israel 1 Itália 4

Figura 27 – Distribuição dos Competidores de Cada País

6

CONCLUSÃO Este trabalho teve por objetivo a comparação de modelos de previsão de carga máxima

diária para o horizonte de curto prazo (um mês), utilizando a análise clássica de séries temporais e a aplicação da Rede Neural Artificial (RNA) feedforward, mais especificamente do algoritmo de treinamento bayesiano de MLP’s. Conforme evidenciado pelos resultados apresentados no capítulo 5, os métodos utilizados apresentam resultados bastante satisfatórios, pois, quando comparados a modelos bem qualificados internacionalmente podem superá-los. O NN(3,0,1) apresentou um MAPE de 1,66 %, deixando boa margem em relação ao primeiro colocado na competição. A análise clássica de séries temporais apresentou desempenho inferior aos oriundos de RNA’s e não se classificaria entre os dez melhores na competição em questão. Fazendo uma análise local nos dois modelos utilizados, o que apresenta melhor rendimento, menor MAPE, menor MAE e menor tempo de processamento, é o modelo gráfico, mas deixando destacado que este necessita de intervenção especializada, desvantagem não presente no modelo automático (AICBIC). Foi visualizado que quando processada a sazonalidade da série por redes neurais utilizando entradas binárias, há uma melhora significativa nos resultados obtidos em relação à que não recebe tratamento sazonal ou, ainda, tratada com a análise espectral via séries de Fourier, em todos os casos houve melhora no desempenho. A seleção do conjunto inicial de entrada pela teoria do caos apresentou melhor resultado quando houve intervenção do usuário para a seleção dos lags (atrasos) em todos as situações. A clusterização, em até três agrupamentos para modelagem do atrator, mostrou trazer um aumento de dezoito a até quarenta vezes o esforço computacional, mostrando-se pouco robusto para previsões que necessitem de agilidade, mas leva consigo uma melhora nos MAPE’s. Em linhas gerais, todos os métodos exibiram um avanço no auxílio da tomada de decisão baseada na carga futura de curto prazo. Além das aplicações práticas mencionadas

64

acima, as competições promovidas recentemente por diversas entidades ao redor do mundo objetivando o estudo e a avaliação de modelos automáticos de previsão evidenciam a relevância do assunto abordado.

7

REFERÊNCIAS BIBLIOGRÁFICAS

[1]. FERREIRA, Vitor Hugo. Desenvolvimento de modelos neurais autônomos para previsão de carga elétrica. 2008. Tese de Doutorado. COPPE/UFRJ. [2]. ALVES, Amando. Análise Econômica de Sistemas Elétricos. Departamento de Engenharia de Produção, Universidade Federal Fluminense. Niterói : s.n., 2002. [3]. HOBBS, B.F., et al. Analysis of the Value for Unit Commitment of Improved Load Forecasts. 1999. pp. 1342-1348. Vol. 14. [4]. RANAWEERA, D.K., KARADY, G.G. e FARMER, R.G. Economic Impact Analysis of Load Forecasting. 1997. pp. 1388-1392. Vol. 13. [5]. DOUGLAS, A.P., et al. Risk Due to Load Forecast Uncertainty in Short Term Power Systems Planning. 1998. pp. 1493-1499. Vol. 13. [6]. VALENZUELA, J., MAZUMDAR, M. e KAPOOR, A. Influence of Temperature and Load Forecast Uncertainty on Estimates of Power Generation Costs. 2000. pp. 668674. Vol. 15. [7]. MORETTIN, Pedro A. e TOLOI, Clélia M. C. Análise de Séries Temporais. São Paulo : Edgard Blücher, 2004. [8]. ONS, OPERADOR NACIONAL DO SISTEMA ELÉTRICO. Procedimentos de Rede. ONS - Operador Nacional do Sistema. [Online] 18 de Junho de 2010. [Citado em: 24 de Junho de 2011.] Módulo 5 - Consolidação da previsão de carga. http://www.ons.org.br. [9]. PIERCE, D. A. Seazonal adjustment when both deterministic and stocastic seasonality are present. Washington, D.C. : Arnold Zellner, 1979. pp. 242-269.

66

[10]. PEREIRA, Dorivan da Cruz. Previsão de Carga Utilizando Análise Clássica de Séries Temporais. Niterói : s.n., 2010. TCC, UFF. [11]. BOX, G. E. P. and JENKINS. Time Series Analysis: Forecasting and Control. Revised edition, 1976. São Francisco : Holden-Day, 1970. [12]. BOX, G. E. P., JENKINS, G.M. and REINSEL, G. Time Series Analysis: Forecasting and Control. Third Edition. Englewood Cliffs : Prentice Hall, 1994. [13]. HANNAN, E. J. Testing for autocorrelation and Akaike's criterion. Sheffield : s.n., 1982. pp. 403-412. Vol. special volume 19A. [14]. AKAIKE, H. Maximum likehood identification of Gaussian autoregressive moving average models. s.l. : Biometrika, 60, 1973. pp. 255-265. [15]. —. A new look at the statistical model indentification. s.l. : IEEE Trans-actions on Automatic Control, AC-19, 1974. pp. 716-723. [16]. HURVICH, C. M. e TSAI, C. L. Regression and time series model selection in small samples. 1989. pp. 297-307. [17]. AKAIKE, H. On Entropy maximization principle. P.R.Krishnaiah. Amsterdan : NorthHolland, 1977. pp. 27-41. [18]. RISSANEM, J. Modeling by shortest data description. 9178. pp. 465-471. [19]. SCHWARZ, G. Estimating the dimension of a model. 1978. pp. 461-464. [20]. HANNAN, E. J. The estimation of the order of an ARMA process. 1980. pp. 1071-1081. [21]. PÖTSCHER, B. M. Estimation of autoregressive moving-average order given an infitite number of models and approximation of spectral densities. 1990. pp. 165-179. [22]. HAYKIN, S. Redes neurais: princípios e prática. 2ª ed. Porto Alegre : Bookman, 2001. [23]. HIPPERT, Henrique Steinherz, PEDREIRA, Carlos Eduardo e SOUZA, Reinaldo Castro. Neural Network for Short-Term Load Forecasting: A Review and Evaluation. 2001. Vol. 16.

67

[24]. BRITO, Pedro Marcondes de. Estimativa da Margem de Estabilidade de Tensão Através de Redes Neurais Artificiais e Simulação Rápida no Tempo. Niterói : s.n., 2010. TCC. UFF. [25]. CHURCHLAND, P. S. e SEJNOWSKI, T. J. The Computational Brain. Cambridge : MIT Press, 1992. [26]. MENDEL, J. M. and MCLAREN, R. W. Reinforcement-learning control and pattern recognition systems. New York : Academic Press. pp. 287-318. Vol. 66. [27]. BECKER, S. Unsupervised learning procedures for neural networks. s.l. : Iternational Journal of Neural Systems, 1991. pp. 267-277. Vol. 2. [28]. FERREIRA, Vitor Hugo. Notas de Aula. 2010. Aulas 16-22. [29]. BISHOP, C. M. Neural Networks for Pattern Recognition. Oxford : Oxford University Press, 1995. [30]. STONE, M. Cross-Validation: A review. 1978. pp. 127-139. Vol. 9. [31]. STONE, M. Cross-Validatory choice and assessment of statistical predictions. 1974. pp. 111-133. Vol. B36. [32]. World-wide competition within the EUNITE network. Eletricity Load Forecast using Intelligent Adaptative Technology. [Online] [Citado em: 25 de Junho de 2011.] http://neuron.tuke.sk/competition/. [33]. FERREIRA, Vitor Hugo. Técnicas de Regularização de Modelos Neurais Aplicadas à Previsão de Carga a Curto Prazo. Rio de Janeiro : s.n., 2005. COPPE/UFRJ. [34]. FRASER, Andrew M. e SWINNEY, Harry L. Independent Coordinates for Strange Attractors from Mutual Information. Austin : The American Physical Society, 1986. Vol. 33. 2. [35]. CAO, Liangyue. Practical method for determining the minimum embedding dimension of a scalar time series. Maio de 1997.