Representação da Incerteza em Modelos de …...UNIVERSIDADE ESTADUAL DE CAMPINAS aculFdade de Engenharia Mecânica ernandaF Nakano Kazama Representação da Incerteza em Modelos

UNIVERSIDADE ESTADUAL DE CAMPINASFaculdade de Engenharia Mecânica

Fernanda Nakano Kazama

Representação da Incerteza em Modelos deProgramação Dinâmica Estocástica Através de Latisse

Binomial: Análise na Perspectiva do DECOMP

CAMPINAS2017

Agência(s) de fomento e nº(s) de processo(s): CNPq, 155510/2015-0

Ficha catalográficaUniversidade Estadual de Campinas

Biblioteca da Área de Engenharia e ArquiteturaLuciana Pietrosanto Milla - CRB 8/8129

Kazama, Fernanda Nakano, 1992- K189r KazRepresentação da incerteza em modelos de programação dinâmica

estocástica através de latisse binomial : análise na perspectiva do DECOMP /Fernanda Nakano Kazama. – Campinas, SP : [s.n.], 2017.

KazOrientador: Paulo de Barros Correia. KazDissertação (mestrado) – Universidade Estadual de Campinas, Faculdade

de Engenharia Mecânica.

Kaz1. Otimização. 2. Sistemas de energia elétrica hidrotérmica. 3.

Programação estocástica. 4. Programação dinâmica. 5. Conditional value atrisk (CVaR). I. Correia, Paulo de Barros,1954-. II. Universidade Estadual deCampinas. Faculdade de Engenharia Mecânica. III. Título.

Informações para Biblioteca Digital

Título em outro idioma: Representation of uncertainty in stochastic dynamic programmingmodels through binomial lattice : analysis from the perspective of DECOMPPalavras-chave em inglês:OptimizationHydrothermal electric power systemsStochastic ProgrammingDynamic programmingConditional valeu at risk (CVaR)Área de concentração: Planejamento de Sistemas EnergéticosTitulação: Mestra em Planejamento de Sistemas EnergéticosBanca examinadora:Paulo de Barros Correia [Orientador]Sergio Valdir BajayMônica de Souza ZambelliData de defesa: 17-02-2017Programa de Pós-Graduação: Planejamento de Sistemas Energéticos

Powered by TCPDF (www.tcpdf.org)

UNIVERSIDADE ESTADUAL DE CAMPINAS

FACULDADE DE ENGENHARIA MECÂNICA

COMISSÃO DE PÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICA

PLANEJAMENTO DE SISTEMAS ENERGÉTICOS

DISSERTAÇÃO DE MESTRADO ACADÊMICO

Representação da Incerteza em Modelos deProgramação Dinâmica Estocástica Através de Latisse

Binomial: Análise na Perspectiva do DECOMP

Autor: Fernanda Nakano Kazama

Orientador: Prof. Dr Paulo de Barros Correia

A Banca Examinadora composta pelos membros abaixo aprovou esta Tese:

Prof. Dr. Paulo de Barros Correia, PresidenteUniversidade Estadual de Campinas

Prof. Dr. Sérgio Valdir BajayUniversidade Estadual de Campinas

Dra. Mônica de Souza ZambelliCompanhia Paulista de Força e Luz

A Ata da defesa com as respectivas assinaturas dos membros encontra-se no processo devida acadêmica do aluno.

Campinas, 17 de Fevereiro de 2017.

Dedicatória

Dedico este trabalho à minha família que sempre acreditou em mim e me apoiou achegar até aqui.

Agradecimentos

Este trabalho não poderia ser terminado sem a ajuda de diversas pessoas às quais prestominha homenagem:

Ao Prof. Dr. Paulo de Barros Correia, orientador deste trabalho.

Aos meus colegas de sala, que sempre se predispuseram a me auxiliar nos pequenos deta-lhes.

À minha família, pelo apoio e incentivo dado.

Aos meus amigos, pelo apoio e companheirismo.

Ao CNPq, pelo apoio �nanceiro.

�Não arriscar nada é arriscar tudo.�Al Gore

Resumo

O planejamento da operação de sistemas hidrotérmicos de geração de energiaelétrica no Brasil utiliza modelos computacionais de otimização (NEWAVE e DECOMP)para determinar operação do sistema, buscando atender a demanda de forma econômicae con�ável. Como o sistema está sujeito a várias incertezas (a�uências, demanda,...),utilizam-se adicionalmente mecanismos de aversão ao risco, como o CVaR (ConditionalValue at Risk), para manter em níveis aceitáveis o risco de não atendimento da demanda.Além disso, a principal estocasticidade do problema é tratada nos modelos NEWAVE eDECOMP como árvores de a�uências, onde cada ramo apresenta uma probabilidade deevolução de um determinado cenário, com um cenário completo representado pelo cami-nho do nó origem até uma das folhas da árvore. Esta representação torna o problemacomputacionalmente muito complexo com o crescimento dos períodos de tempo analisa-dos, uma vez que os ramos da árvore de a�uências aumentam de forma exponencial como tempo, podendo cair na chamada maldição da dimensionalidade.

Assim, esta dissertação propõe uma melhoria para a representação da incer-teza no modelo DECOMP, através da substituição da árvore de a�uências por latisse, queapresenta um crescimento linear a cada período analisado, diminuindo o esforço compu-tacional e possibilitando a análise de mais períodos de tempo.

Palavras-chave: Otimização, Sistemas de Energia Elétrica Hidrotérmica, ProgramaçãoEstocástica, Programação Dinâmica, Conditional Value at Risk (CVaR), Árvore, Latisse.

Abstract

The planning of hydrothermal systems operation of electric power generation inBrazil uses computational optimization models (NEWAVE and DECOMP) to determinethe system operation, seeking to attend the demand in an economically and reliably way.As the system is subject to various uncertainties (a�uent energy, demand,...), it is usedin addition risk aversion mechanisms, such as CVaR (Conditional Value at Risk), tomaintain the risk of non-ful�llment of the demand in acceptable levels. Furthermore, thestochasticity of the problem is treated by NEWAVE and DECOMP models though in�owtrees, where each branch represents a probability of evolution of a given scenario, and acomplete scenario is represented by the path from the root node to one of the leaves of thetree. This representation makes the problem more computationally complex as greateras the number of periods of time analyzed, since the branches of the in�ow tree increaseexponentially with time, which can result in what is called curse of dimensionality.

Thus, this thesis proposes an improvement of the representation of uncertaintyin DECOMP model, by replacing the in�ow tree for lattice, which features a linear growthin the analyzed period, reducing the computational e�ort and making possible to analyzemore periods of time.

Key-words: Optimization, Hydrothermal Electric Power Systems, Stochastic Program-ming, Dynamic Programming, Conditional Value at Risk (CVaR), Tree, Lattice.

Lista de Figuras

1.1 Sistema hidrotérmico equivalente. . . . . . . . . . . . . . . . . . . . . . . . 171.2 Esquema de representação dos cenários de a�uência no modelo DECOMP. 191.3 Análise do número de ramos e nós por estágio de uma árvore binária. . . . 191.4 Estrutura de uma latisse binomial. . . . . . . . . . . . . . . . . . . . . . . 201.5 Análise do número de ramos e nós por estágio de uma latisse binomial. . . 212.1 Esquema de uma usina hidroelétrica (ANEEL, 2008). . . . . . . . . . . . . 232.2 Representação dos parâmetros de uma hidroelétrica. . . . . . . . . . . . . . 242.3 Esquema de usinas hidroelétricas situadas numa mesma bacia. . . . . . . . 252.4 Geração de energia por meio de uma termoelétrica em ciclo à vapor. . . . . 262.5 Geração de energia por meio de uma termoelétrica com turbina a gás

(ANEEL, 2008). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.6 Geração de energia através de uma termoelétrica em ciclo combinado (FUR-

NAS, 2016). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.7 Representação do custo de geração térmica (PEREIRA, 2000). . . . . . . . 292.8 Curva de carga típica de um dia (ALVES, 2007). . . . . . . . . . . . . . . . 302.9 Curva de carga do subsistema nordeste em 2015 (ONS, 2016). . . . . . . . 302.10 Variação da incerteza ao longo do tempo (BERTHO, 2013). . . . . . . . . . 312.11 Divisão dos subsistemas brasileiros (BERTHO, 2013). . . . . . . . . . . . . 322.12 Esquema da geração de séries sintéticas de energia. . . . . . . . . . . . . . 332.13 Esquema de um reservatório equivalente. . . . . . . . . . . . . . . . . . . . 342.14 Entradas e saídas de dados nos modelos NEWAVE e DECOMP (CPAMP,

2013b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.1 Séries em paralelo (pente) (CEPEL, 2010). . . . . . . . . . . . . . . . . . . 373.2 Séries em árvore (CEPEL, 2010). . . . . . . . . . . . . . . . . . . . . . . . 383.3 Representação dos cenários hidrológicos através do modelo GEVAZP. . . . 383.4 Distribuição de probabilidades dos valores históricos de energia natural

a�uente de Janeiro para a bacia do Tocantins e sua respectiva árvore. . . . 393.5 Esquema de cálculo do balanço hídrico na árvore binária. . . . . . . . . . . 403.6 Esquema de cálculo do valor esperado do custo de operação na árvore binária. 413.7 Esquema de cálculo do valor esperado do custo de operação na árvore ternária. 423.8 Esquema de cálculo do balanço hídrico na árvore multiária. . . . . . . . . . 433.9 Esquema de cálculo do valor esperado do custo de operação na árvore

multiária. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.10 Esquema de cálculo do balanço hídrico na latisse binomial. . . . . . . . . . 453.11 Esquema de cálculo do valor esperado du custo de geração na latisse binomial. 463.12 Esquema de cálculo do valor esperado do custo de operação na latisse tri-

nomial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.13 Esquema de cálculo do balanço hídrico na latisse multinomial. . . . . . . . 48

3.14 Esquema de cálculo do valor esperado do custo de operação na latisse mul-tinomial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4.1 Representação de uma árvore binomial de dois estágios. . . . . . . . . . . . 514.2 Representação de uma árvore binária multi-estágios. . . . . . . . . . . . . . 524.3 Representação dos subproblemas de uma árvore binomial multi-estágios. . 544.4 Representação de uma latisse binomial multi-estágios. . . . . . . . . . . . . 574.5 Representação dos subproblemas de uma latisse binomial multi-estágios. . . 595.1 Distribuição de retornos (CCEE, 2012). . . . . . . . . . . . . . . . . . . . . 635.2 Representação do CVaR numa curva de distribuição de retornos (CCEE,

2012). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635.3 Distribuição de perdas (CCEE, 2012). . . . . . . . . . . . . . . . . . . . . . 645.4 Distribuição de a�uências. . . . . . . . . . . . . . . . . . . . . . . . . . . . 645.5 Distribuição de custos esperados. . . . . . . . . . . . . . . . . . . . . . . . 655.6 Representação dos cenários mais críticos para uma árvore binária multi-

estágios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 665.7 Representação dos subproblemas de uma árvore binomial multi-estágios

aversa ao risco. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.8 Representação dos cenários mais críticos para uma latisse binomial multi-

estágios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.9 Representação dos subproblemas de uma latisse binomial multi-estágios

aversa ao risco. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 706.1 Usinas hidroelétricas na bacia do Tocantins e suas respectivas capacidades

máximas de geração (ONS, 2014). . . . . . . . . . . . . . . . . . . . . . . . 736.2 Esquema do problema a ser resolvido. . . . . . . . . . . . . . . . . . . . . . 756.3 Esquema da árvore de a�uências utilizada no problema. . . . . . . . . . . . 766.4 Esquema da latisse binomial utilizada no problema. . . . . . . . . . . . . . 766.5 Solução para o caso indiferente ao risco utilizando árvore. . . . . . . . . . . 776.6 Solução para o caso averso ao risco utilizando árvore. . . . . . . . . . . . . 776.7 Solução para o caso indiferente ao risco utilizando latisse. . . . . . . . . . . 786.8 Solução para o caso averso ao risco utilizando latisse. . . . . . . . . . . . . 786.9 Grá�co comparativo entre os volumes do reservatório da árvore e da latisse

para o caso neutro ao risco. . . . . . . . . . . . . . . . . . . . . . . . . . . 796.10 Grá�co comparativo entre os volumes do reservatório para o caso neutro e

averso ao risco. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 796.11 Solução para o caso indiferente ao risco utilizando árvore, com volume

inicial de 9.000MWmed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.12 Solução para o caso indiferente ao risco utilizando latisse, com volume

inicial de 9.000MWmed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.13 Grá�co comparativo entre os volumes do reservatório da árvore e da latisse

para o caso neutro ao risco, com volume inicial de 5.000MWmed. . . . . . 81

A.1 Representação das probabilidades de uma árvore binária multi-estágios. . . 86A.2 Representação das probabilidades de uma latisse binomial multi-estágios. . 87

Lista de Tabelas

1.1 Número de nós �nais e ramos de cada estágio de uma árvore. . . . . . . . . 191.2 Número de nós �nais de cada estágio de uma latisse. . . . . . . . . . . . . 211.3 Número de ramos de uma latisse a cada estágio. . . . . . . . . . . . . . . . 216.1 Custo de geração das termoelétricas. . . . . . . . . . . . . . . . . . . . . . 736.2 Dados da demanda de energia. . . . . . . . . . . . . . . . . . . . . . . . . . 746.3 Cenários de a�uência. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

B.1 Resultados obtidos pela otimização utilizando árvore para o caso neutro aorisco, com volume inicial equivalente à 10.000MWmed. . . . . . . . . . . . 89

B.2 Resultados obtidos pela otimização utilizando latisse para o caso neutro aorisco, com volume inicial equivalente à 10.000MWmed. . . . . . . . . . . . 89

B.3 Resultados obtidos pela otimização utilizando árvore para o caso averso aorisco, com volume inicial equivalente à 10.000MWmed. . . . . . . . . . . . 90

B.4 Resultados obtidos pela otimização utilizando latisse para o caso averso aorisco, com volume inicial equivalente à 10.000MWmed. . . . . . . . . . . . 90

B.5 Resultados obtidos pela otimização utilizando árvore para o caso neutro aorisco, com volume inicial equivalente à 9.000MWmed. . . . . . . . . . . . . 91

B.6 Resultados obtidos pela otimização utilizando latisse para o caso neutro aorisco, com volume inicial equivalente à 9.000MWmed. . . . . . . . . . . . . 91

Lista de Abreviaturas e Siglas

Letras Latinas

a - queda bruta da usinaA - a�uênciaAinc - a�uência incrementalc - custo de geração das usinas por unidade de energiacf - nível do canal de fugaCTO - custo total da operaçãoD - demandaEV - evaporaçãoF - in�ltraçãog - aceleração da gravidadeGk - potência máxima da usina kGH - geração hidroelétrica (vazão de água turbinada ou

energia gerada)GT - geração termoelétricaH - potência (energia instantânea)I - vertimentoK - número de usinasn - número de cenáriosN - número de estágiosr - distribuição de retornos esperadosS - valor esperado do custo de operaçãoT - temperaturaV - volume do reservatório da hidroelétricax - vetor coluna de geração de eletricidade das usinasX - elemento do vertor x (geração de uma usina)x∗ - resposta encontraday - função objetivoZ - vazão devido à agua turbinada ou vertida nos reser-

vatórios à montante

Letras Gregas

α - custo esperado da geração futuroα e λ - parâmetros do CVaRαβ - value at riskβ - nível de con�ançaγ - fator de conversão de vazão para volumeη - rendimentoη - rendimento médioηmax - rendimento máximoπ - multiplicador simplex

Subscritos

f - frioij - ramoi ou j - nók - usinam,M - posicionamento dos ramosn - número de ramosq - quentes - índice que representa o tipo de cenário

Subrescritos

t - estágio

Siglas

ANEEL - Agência Nacional de Energia ElétricaCCEE - Câmara de Comercialização de Energia ElétricaCVAR - Conditional Value at RiskCPAMP - Comissão Permanente para Análise de Metodologias e Programas Com-

putacionais do Setor ElétricoCVU - Custo Variável UnitárioENA - Energia Natural A�uenteMAR - Mecanismo de Aversão ao RiscoPDDE - Programação Dinâmica Dual EstocásticaVAR - Value at Risk

Sumário

1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.1 Contexto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.2 Objetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201.3 Organização da Dissertação . . . . . . . . . . . . . . . . . . . . . . 22

2 OPERAÇÃO DE SISTEMAS HIDROTÉRMICOS . . . . . . . . . . . . . . 232.1 Geração Hidroelétrica . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.1.1 Volume Armazenado . . . . . . . . . . . . . . . . . . . . . 242.1.2 Energia A�uente . . . . . . . . . . . . . . . . . . . . . . . 242.1.3 Energia Gerada . . . . . . . . . . . . . . . . . . . . . . . . 25

2.2 Geração Termoelétrica . . . . . . . . . . . . . . . . . . . . . . . . . 262.2.1 Capacidade Agregada de Geração . . . . . . . . . . . . . . 282.2.2 Custo Agregado de Geração . . . . . . . . . . . . . . . . . 29

2.3 Curva de Carga . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.4 Modelos de Otimização Utilizados no Brasil . . . . . . . . . . . . . 31

2.4.1 NEWAVE . . . . . . . . . . . . . . . . . . . . . . . . . . . 322.4.2 DECOMP . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3 REPRESENTAÇÃO DOS CENÁRIOS HIDROLÓGICOS . . . . . . . . . 373.1 GEVAZP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.2 Representação por Árvore . . . . . . . . . . . . . . . . . . . . . . . 39

3.2.1 Árvore Binária . . . . . . . . . . . . . . . . . . . . . . . . 403.2.2 Árvore Ternária . . . . . . . . . . . . . . . . . . . . . . . . 413.2.3 Árvore Multiária . . . . . . . . . . . . . . . . . . . . . . . 42

3.3 Representação por Latisse . . . . . . . . . . . . . . . . . . . . . . . 443.3.1 Latisse Binomial . . . . . . . . . . . . . . . . . . . . . . . 443.3.2 Latisse Trinomial . . . . . . . . . . . . . . . . . . . . . . . 463.3.3 Latisse Multinomial . . . . . . . . . . . . . . . . . . . . . 47

4 MODELO DINÂMICO DE OTIMIZAÇÃO COM INDIFERENÇA AORISCO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.1 Programação Dinâmica Dual Estocástica . . . . . . . . . . . . . . . 504.2 Modelo Integrado com Árvore Binária . . . . . . . . . . . . . . . . 524.3 Modelo Decomposto com Árvore Binária . . . . . . . . . . . . . . . 544.4 Modelo Integrado com Latisse Binomial . . . . . . . . . . . . . . . 574.5 Modelo Decomposto com Latisse Binomial . . . . . . . . . . . . . . 59

5 MODELO DINÂMICO DE OTIMIZAÇÃO COM AVERSÃO AO RISCO . 625.1 Conditional Value at Risk . . . . . . . . . . . . . . . . . . . . . . . 625.2 Modelo Integrado com Árvore Binária . . . . . . . . . . . . . . . . 655.3 Modelo Decomposto com Árvore Binária . . . . . . . . . . . . . . . 675.4 Modelo Integrado com Latisse Binomial . . . . . . . . . . . . . . . 68

5.5 Modelo Decomposto com Latisse Binomial . . . . . . . . . . . . . . 696 ESTUDO DE CASO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

6.1 Metodologia e Dados do Problema . . . . . . . . . . . . . . . . . . 726.2 Resolução e resultados . . . . . . . . . . . . . . . . . . . . . . . . . 76

7 CONCLUSÕES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 827.1 Sugestões para Trabalhos Futuros . . . . . . . . . . . . . . . . . . . 83

Referências . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

APÊNDICE 86A Cálculo das Probabilidades de Ocorrência de Cada Cenário . . . . . . . . . 86

A.1 Árvore de A�uências . . . . . . . . . . . . . . . . . . . . . . . . . . 86A.2 Latisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

B Resultados Obtidos para o Estudo de Caso . . . . . . . . . . . . . . . . . . 89

17

1 INTRODUÇÃO

1.1 Contexto

Segundo Castro et al. (2012), a essencialidade do �bem energia� no desenvol-

vimento das atividades econômicas e sociais do mundo contemporâneo é inquestionável,

sendo imprescindível o planejamento do setor energético para atender à crescente de-

manda por energia. Uma maior complexidade deste planejamento apresenta-se no âmbito

do setor elétrico, uma vez que o produto �nal (eletricidade) não é estocável e, portanto,

deve haver um equilíbrio instantâneo entre oferta e demanda.

A energia elétrica é obtida através da transformação de diferentes formas de

energia em eletricidade. Em uma usina hidroelétrica, por exemplo, a eletricidade é resul-

tante do aproveitamento da energia potencial da água de um rio, apresentando, no Brasil,

um custo de operação quase nulo. Em usinas termoelétricas, por outro lado, utiliza-se

a energia proveniente da queima de combustíveis diversos, tendo seu custo de operação

variando conforme o preço dos insumos e o rendimento da termoelétrica. Além destas,

existem ainda outras fontes de obtenção de energia elétrica, como a eólica e a fotovoltaica,

mas, de modo geral, a maior parte da eletricidade no Brasil é gerada a partir de usinas

hidroelétricas e termoelétricas, que são conectadas aos centros de carga através de um

sistema de transmissão, formando um sistema hidrotérmico, como mostrado na Figura

1.1.

Hidroeletrica Termoeletrica

Rede deTransmissao

Centrosde

Carga

Figura 1.1: Sistema hidrotérmico equivalente.

18

A �m de diminuir os custos de geração, o planejamento da operação de siste-

mas hidrotérmicos apresenta como objetivo econômico substituir, na medida do possível,

a geração de origem termoelétrica por geração de origem hidroelétrica. Para isso, é pre-

ciso levar em consideração que a otimização do sistema apresenta características de um

problema dinâmico, interconectado, estocástico, não linear e de grande porte (SOARES,

1987).

O problema é dinâmico, uma vez que os recursos hidroelétricos de geração,

representados pela água armazenada nos reservatórios, são limitados, e a sua disponibi-

lidade num dado momento depende do grau de sua utilização anterior; interconectado,

pois as usinas hidroelétricas situadas na mesma bacia hidráulica apresentam acoplamento

operativo entre si; estocástico, devido à incerteza sobre as vazões a�uentes futuras e a

demanda de energia; não linear, devido à função de produção hidráulica e ao custo de ge-

ração das usinas termoelétricas; e de grande porte, devido ao número de usinas e restrições

a serem consideradas na resolução.

Assim, de modo a otimizar os custos esperados de geração e ainda atender à

demanda de energia ao longo de um certo período de tempo, é necessário determinar,

para cada etapa, quais usinas devem ser acionadas e suas respectivas metas de geração,

levando em conta todas as restrições envolvidas no problema.

Os modelos NEWAVE e DECOMP (CEPEL, 2012a), utilizados para o pla-

nejamento da operação a médio e curto prazo, respectivamente, no Brasil, consideram o

problema de otimização com o critério de mínimo custo, onde o objetivo é obter uma po-

lítica de operação que minimize o valor esperado da geração térmica e eventuais cortes de

carga, considerando-se um dado conjunto de possíveis cenários de a�uências futuras aos

reservatórios. Entretanto, a natureza estocástica das vazões a�uentes faz com que nem

todos os cenários resultem em garantia de suprimento. Logo, para auxiliar os modelos na

mitigação de risco de não suprimento foram adicionados ferramentas e critérios de aversão

ao risco. Essa ferramenta utilizada atualmente é o Valor Condicionado a um Dado Risco,

ou do inglês, Conditional Value at Risk (CVaR) (CPAMP, 2013a).

A estocasticidade das vazões é introduzida nos modelo DECOMP ao �nal do

primeiro período de estudos através da construção de cenários de vazões a�uentes às

usinas do sistema. Estes cenários hidrológicos podem ser representados através de uma

árvore de a�uências, com probabilidades de ocorrência associadas a cada ramo. Para a

montagem da árvore, deve-se determinar o número de cenários (ramos) à ser analisado em

cada período de tempo (estágio). A Figura 1.2 apresenta uma estrutura de cenários para

um horizonte de estudo de 4 meses, em que o primeiro mês está dividido em 4 semanas

com a�uências supostas conhecidas.

19

1 2 3 4 5 6 7 8Primeiro mes Demais meses

Figura 1.2: Esquema de representação dos cenários de a�uência no modelo DECOMP.

Analisando-se a parte estocástica da representação, pode-se observar que a

cada estágio da árvore de a�uências ocorre um crescimento exponencial do número de

nós e ramos, havendo apenas um possível caminho do nó inicial para cada nó �nal, como

mostrado na Figura 1.3 e na Tabela 1.1.

2

48

1

8

24

Figura 1.3: Análise do número de ramos e nós por estágio de uma árvore binária.

Tabela 1.1: Número de nós �nais e ramos de cada estágio de uma árvore.

Quantidade de ramospor nó �nal

Estágios1 2 3 4 N

2 2 4 8 16 2N

3 3 9 27 81 3N

4 4 16 64 256 4N

n n n2 n3 n4 nN

20

Este crescimento exponencial dos cenários di�culta a análise do problema para

períodos muito longos devido ao grande esforço computacional exigido. Assim, a �m de

contornar esta situação, o modelo NEWAVE utiliza a técnica de Programação Dinâmica

Dual Estocástica (PDDE) (CEPEL, 2012b) para subdividir o problema, aumentando-se

assim, o conjunto de problemas que são viáveis de se solucionar computacionalmente.

Porém, mesmo com um maior número de problemas do despacho hidrotérmico

sendo viáveis de se solucionar graças à utilização da técnica de PDDE, o crescimento do

esforço computacional continua sendo exponencial com o aumento dos períodos analisados,

tornando-se interessante estudar outros métodos ou técnicas de resolução do problema de

otimização do despacho hidrotérmico.

1.2 Objetivo

Dado o contexto, esta dissertação irá abordar a otimização do despacho hi-

drotérmico, dando ênfase à forma como a incerteza é representada no modelo DECOMP,

propondo a substituição da árvore de a�uências pela latisse de volumes, cuja estrutura é

mostrada na Figura 1.4.

1 2 3 4 5 6 7 8Primeiro mes Demais meses

Figura 1.4: Estrutura de uma latisse binomial.

Ao contrário da estrutura de árvore que apresenta um crescimento exponencial,

a latisse apresenta uma tendência linear de aumento do número de nós e ramos com o

passar do tempo, havendo Ct,s possíveis caminhos do nó inicial até cada um dos nós �nais,

sendo C a combinação de t estágios de acordo com o posicionamento s de cada nó, como

mostrado na Figura 1.5 e nas Tabelas 1.2 e 1.3.

21

2

46

1

2

34 C3,0 =

3!3!·0! = 1

C3,1 =3!

2!·1! = 3

C3,2 =3!

1!·2! = 3

C3,3 =3!

0!·3! = 1

Figura 1.5: Análise do número de ramos e nós por estágio de uma latisse binomial.

Tabela 1.2: Número de nós �nais de cada estágio de uma latisse.

Quantidade de ramospor nó inicial

Estágio1 2 3 4 N

2 2 3 4 5 N + 13 3 5 7 9 2 ·N + 14 4 7 10 13 3 ·N + 1n n 2 · n− 1 3 · n− 2 4 · n− 3 (n− 1) ·N + 1

Tabela 1.3: Número de ramos de uma latisse a cada estágio.

Quantidade de ramospor nó inicial

Estágio1 2 3 4 N

2 2 4 6 8 2 ·N3 3 9 15 21 6 ·N − 34 4 16 28 40 12 ·N − 8n n n2 2 · n2 − n 3 · n2 − 2 · n (N − 1) · n2 − (N − 2) · n

Comparando-se as duas estruturas, a árvore binária apresenta, ao �nal do

terceiro estágio, 8 possíveis cenários de a�uências (caminhos diferentes do primeiro nó

até um nó �nal) e 8 possíveis estados �nais (nós terminais), pois há apenas um caminho

que leva a cada nó. Enquanto a latisse binomial, no mesmo estágio, apresenta também 8

possíveis cenários, porém apenas 4 estados �nais, uma vez que há mais de um caminho

diferente que pode levar à um mesmo nó.

Desse modo, se provada a viabilidade da substituição da estrutura de árvore

de a�uências pela latisse, poder-se-ia resolver o problema de otimização do despacho

hidrotérmico utilizando-se um menor esforço computacional, mantendo-se ainda o número

22

de cenários analisados e a continuidade dos valores assumidos em cada estado. Assim, o

objetivo desta dissertação é veri�car a viabilidade desta substituição.

Os testes de melhoria no gasto de tempo computacional não estão no escopo

desta dissertação, �cando como sugestão para trabalhos futuros.

1.3 Organização da Dissertação

O primeiro capítulo contém o contexto, o objetivo e a estrutura da dissertação.

O capítulo 2, por sua vez, apresenta a revisão bibliográ�ca sobre a operação de sistemas

hidrotérmicos e, em seguida, o capítulo 3 é abordado apresentando a representação dos

cenários hidrológicos através dos conceitos de árvore e latisse. Os capítulos 4 e 5 apresen-

tam os modelos dinâmicos de otimização com indiferença e com aversão ao risco, sendo

mostrado para isso uma breve explicação dos conceitos de Programação Dinâmica Dual

Estocástica e Conditional Value at Risk. Dados os fundamentos, o capítulo 6 apresenta

um estudo de caso à �m de avaliar a proposta de substituição da árvore de a�uências por

latisse. Por �m, no capítulo 7 são apresentadas as conclusões desta dissertação.

23

2 OPERAÇÃO DE SISTEMAS HIDROTÉRMICOS

Um sistema hidrotérmico é formado por hidroelétricas e termoelétricas que

operam de forma conjugada de modo a atender a demanda de energia. Assim, para a

melhor compreensão da operação de sistemas hidrotérmicos no Brasil, é necessário ter-se

conhecimento sobre a geração hidroelétrica, termoelétrica, a curva de carga e os modelos

de planejamento da operação (NEWAVE e DECOMP), que são apresentados nas seções

à seguir.

2.1 Geração Hidroelétrica

Segundo Kligerman (1992) uma usina hidroelétrica é composta basicamente

por: uma barragem formadora do reservatório, que represa um curso d'água; condutos

forçados que levam a água do reservatório até a casa de força; a casa de força, onde estão

instalados os grupos turbina-gerador; um canal de fuga, para onde a água é restituída

ao seu curso natural após turbinada; e o vertedouro, por onde há liberação da água

diretamente, sem passar pelas turbinas; como mostrado no per�l esquemático da Figura

2.1.

Figura 2.1: Esquema de uma usina hidroelétrica (ANEEL, 2008).

Os reservatórios das usinas hidroelétricas, conforme sua capacidade de regula-

rização, são classi�cados como de compensação, quando tem volume útil su�ciente apenas

para regularizar descargas de poucos dias; ou de acumulação, quando são capazes de pro-

mover a regularização por um período de um a vários meses. As usinas com reservatórios

de compensação são chamadas usinas a �o d'água. Aquelas com reservatório de acumu-

lação são chamadas usinas de reservatório (KLIGERMAN, 1992).

24

Assim, a geração hidroelétrica ocorre por meio da conversão da energia poten-

cial da água no reservatório em energia cinética e de pressão na turbina, que aciona o

gerador e, por sua vez, converte a energia cinética e de pressão em elétrica.

Desse modo, os principais parâmetros da operação de uma usina hidroelétrica

estão relacionadas ao volume de água armazenada nos reservatórios, à energia a�uente à

usina e à energia gerada, que serão apresentadas a seguir.

2.1.1 Volume Armazenado

No caso das usinas com reservatório de acumulação, a equação dinâmica que

determina a evolução do nível de água armazenada no reservatório deve satisfazer o princí-

pio da conservação de massa, ou seja, o volume de água no reservatório deve ser equivalente

ao volume de água que já está no reservatório somado ao volume de água que chega neste,

subtraído do volume de água que sai do mesmo.

Assim, para a representação esquemática da Figura 2.2, onde foram descon-

sideradas as perdas de água do reservatório devido à desvios e usos múltiplos, tem-se a

Equação (2.1), em que V t representa o volume do reservatório no instante t e V t−1 repre-

senta o volume no instante anterior, At a vazão a�uente ao reservatório até o instante t,

GH t a vazão de água turbinada para produzir energia no instante t, I t a vazão de água

vertida, EV t a vazão de água evaporada, F t a vazão de água in�ltrada e γ o fator de

conversão de vazão para volume.

At

V t

EV t

F t

It

GHt

Figura 2.2: Representação dos parâmetros de uma hidroelétrica.

V t − V t−1 = (At −GH t − I t − EV t − F t) · γ (2.1)

2.1.2 Energia A�uente

No modelo NEWAVE, a vazão de água Atk que chega à um reservatório de uma

usina k num instante t é convertida em unidade de energia, sendo chamada de energia

a�uente. Essa energia é proveniente da vazão Ztk, correspondente à agua turbinada ou

25

vertida nos reservatórios das usinas k à montante, sendo portanto esta parcela controlável;

e da vazão incremental Atinc,k, caso exista algum rio a�uente ao reservatório, essa vazão

varia conforme a ocorrência de chuvas, sendo portanto estocástica e não-controlável.

Assim, para as usinas 1 e 3 apresentadas na Figura 2.3, as vazões At1 e At3apresentam apenas a parcela estocástica e não-controlável devido à ocorrência de chuvas,

e para as usinas 2 e 4, as vazões a�uentes At2 e At4 apresentam tanto a parcela controlável

como a não controlável, sendo representadas pelas Equações (2.2) e (2.3).

Usina 1

Usina 2

Usina 3

Usina 4

At1Zt

1

Atinc,2

At2Zt

2

Atinc,4

At4

Zt4

At3Zt

3

Figura 2.3: Esquema de usinas hidroelétricas situadas numa mesma bacia.

At2 = Zt1 + Atinc,2 (2.2)

At4 = Zt2 + Zt

3 + Atinc,4 (2.3)

E a energia que se obtém quando essas vazões são turbinadas nas usinas si-

tuadas rio-abaixo, a partir de um ponto de observação, é chamada de Energia Natural

A�uente (ENA).

2.1.3 Energia Gerada

Por �m, a potência H t gerada por uma usina hidroelétrica é proporcional à

energia potencial da água no reservatório e é dada pela Equação (2.4), em que ηt é o

rendimento da usina (conjunto turbina gerador), GH t é a massa de água turbinada para

produzir energia no instante t, g é a aceleração da gravidade (considerada constante) e

até a queda líquida da usina, que é função do volume de água armazenado V t no instante

t e do nível do canal de fuga cf t.

26

H t = ηt(GH t, at) ·GH t · g · at(V t, cf t) (2.4)

Como simpli�cações, pode-se considerar que para uma hidroelétrica com re-

servatório de grande porte, a altura de queda at é constante no curto prazo e, além disso,

pode-se de�nir um rendimento médio η para o período analisado, obtendo-se assim a

Equação (2.5).

H t = η ·GH t · g · a (2.5)

2.2 Geração Termoelétrica

Quanto à geração termoelétrica, há diversos tipos de usinas, variando de acordo

com o combustível, a turbina utilizada e o ciclo termodinâmico aplicado.

Nas usinas termoelétricas em ciclo a vapor, a primeira etapa consiste na queima

de um combustível fóssil, como carvão, óleo ou gás, transformando a água em vapor com

o calor gerado na caldeira. A segunda etapa consiste na utilização deste vapor, em alta

pressão, para girar a turbina que, por sua vez, aciona o gerador elétrico. Na terceira

etapa, o vapor é condensado, transferindo o resíduo de sua energia térmica para um

circuito independente de refrigeração, retornando a água à caldeira, completando o ciclo,

como mostra a Figura 2.4 (FURNAS, 2016).

Figura 2.4: Geração de energia por meio de uma termoelétrica em ciclo à vapor.

Já nas usinas termoelétricas com turbina a gás, o ar é comprimido no compres-

sor e direcionado à câmara de combustão, onde é misturado ao combustível, possibilitando

sua queima. Ao sair da câmara de combustão, os gases, à alta pressão e temperatura, se

27

expandem conforme passam pela turbina, acionando o gerador elétrico, como mostrado

no esquema de uma usina termoelétrica a gás natural da Figura 2.5.

Figura 2.5: Geração de energia por meio de uma termoelétrica com turbina a gás (ANEEL,2008).

Tem-se que a e�ciência máxima teórica ηmax de uma usina termoelétrica é

dada pela equação de Carnot representada na Equação (2.6), onde Tf e Tq representam a

temperatura na saída do condensador e na entrada da turbina, respectivamente.

ηmax = 1− TfTq

(2.6)

Desse modo, como as usinas termoelétricas em ciclo a vapor apresentam a

temperatura Tf variando entre 320K e 350K e Tq variando entre 640K e 700K, a e�ciência

máxima teórica destas usinas varia entre 37 e 54%. Já as usinas termoelétricas com

turbina a gás apresentam a temperatura Tf variando entre 500K e 550K e Tq variando

entre 950K e 1000K, assim, a e�ciência máxima teórica destas usinas estaria entre 35 e

44% (TOLMASQUIM, 2016).

Um modo de se obter uma maior e�ciência seria por meio de uma usina termo-

elétrica em ciclo combinado, que utiliza uma turbina a vapor e outra a gás, apresentando

valores mais baixos de Tf e valores mais altos de Tq, e assim, consequentemente, obtendo-

se uma maior e�ciência máxima teórica entre entre 50 e 61% (TOLMASQUIM, 2016).

O esquema da Figura 2.6 mostra o funcionamento de uma usina termoelétrica em ciclo

combinado.

28

Figura 2.6: Geração de energia através de uma termoelétrica em ciclo combinado (FUR-NAS, 2016).

Assim, as principais características da operação de um parque de usinas ter-

moelétricas são a capacidade agregada de geração e o custo agregado, apresentados a

seguir.

2.2.1 Capacidade Agregada de Geração

A capacidade agregada de geração do parque de usinas termoelétricas tem

limitação superior da sua produção dada pela soma da capacidade máxima de todas as

usinas que compõe o parque e inferior dada pela soma das energias mínimas que as usinas

29

são obrigadas a produzir.

O limite inferior nem sempre é zero devido às restrições operativas das usinas,

por exemplo, caso uma termoelétrica esteja desligada, para que ela volte a funcionar é

necessário custos iniciais com combustíveis e tempo até que se alcance a temperatura de

operação, desse modo, torna-se mais vantajoso para algumas usinas que estas estejam

sempre operando.

2.2.2 Custo Agregado de Geração

As usinas de um parque termoelétrica são ordenadas em termo de custos mar-

ginais crescentes e a entrada em operação de cada termoelétrica é dada por ordem de

custo unitário, em que o último gerador alocado é chamado de gerador marginal, uma vez

que fazendo-se o ajuste marginal de carga, este de�ne o custo marginal de operação.

Todos os geradores com custo inferior ao marginal operam em sua capacidade

máxima e os geradores com custo de operação superior ao do marginal não entram em

operação (CEPEL, 2012a), ou seja, uma unidade térmica só entra em operação quando a

de custo imediatamente inferior estiver utilizando sua capacidade máxima. Desta forma,

o custo de operação do sistema termoelétrico no tempo t pode ser obtido por meio da

função de geração das termoelétricas que, em geral, é representada por uma função linear

por partes, como ilustrado na Figura 2.7 (PEREIRA, 2000).

Figura 2.7: Representação do custo de geração térmica (PEREIRA, 2000).

2.3 Curva de Carga

A oferta de eletricidade deve coincidir com a demanda, havendo um equílibrio

instantâneo, dado a di�culdade do armazenamento de energia elétrica, que se encontra no

fato da demanda de potência de um sistema sofrer variações de intensidade e característica

ao longo de um período, visto que a utilização das cargas dentro desse sistema é dinâmica,

havendo períodos de maior ou menor demanda.

30

A carga de um sistema é de�nida como a potência média recebida durante um

determinado intervalo de tempo denominado de intervalo de carga. Os valores usuais de

intervalo de carga são 15 minutos, 30 minutos, 1 hora, ou até maiores, sendo o primeiro

mais comum. A relação carga versus tempo é denominada curva de carga e pode ser

diária, semanal, mensal ou anual, conforme seja o período a que se re�ram. A área sob a

curva de carga corresponde a energia consumida no período representado (ALVES, 2007).

Existem grandes variações de carga ao longo do dia. Normalmente, o horário

de menor consumo é durante a madrugada e o de maior consumo começa às 18 horas,

como mostrado na Figura 2.8. Porém, com o aumento da temperatura, o horário de pico

tem se deslocado para o meio da tarde.

Figura 2.8: Curva de carga típica de um dia (ALVES, 2007).

Há ainda uma variação da curva de carga conforme os dias da semana. E com

relação a variação da carga durante os meses, a Figura 2.9 mostra a curva de carga média

mensal para o subsistema nordeste para o ano de 2015.

Figura 2.9: Curva de carga do subsistema nordeste em 2015 (ONS, 2016).

31

Assim, para que se possa ter o atendimento da demanda é necessário ter-se

conhecimento sobre a curva de carga.

2.4 Modelos de Otimização Utilizados no Brasil

Tendo conhecimento do funcionamento dos sistemas hidroelétricos e termo-

elétrico e a carga a ser atendida, pode-se agora apresentar os modelos de otimização

utilizados no Brasil para planejamento da operação de sistemas hidrotérmicos.

Em cada etapa do planejamento da operação de sistemas hidrotérmicos de

geração de energia elétrica são utilizados modelos com diferentes graus de detalhamento

para representação do sistema, abrangendo períodos de estudos com horizontes distintos

(médio prazo, curto prazo e programação diária), uma vez que com a diminuição do

horizonte de planejamento, a incerteza das vazões a�uentes são reduzidas, como mostra

a Figura 2.10.

26

3 Modelos Utilizados no Planejamento da Operação

Energética no Brasil

Esta seção apresenta uma breve discussão acerca dos modelos

computacionais utilizados atualmente no planejamento da operação energética

brasileira.

A designação planejamento energético é utilizada quando o aspecto energético

da operação é o foco do planejamento. Com isso, sua principal meta é o atendimento

da demanda de energia envolvida na operação do sistema, para horizontes de

planejamento de médio e curto prazo, ou seja, horizontes plurianuais. O planejamento

elétrico da operação diz respeito ao planejamento no horizonte de curtíssimo prazo, e

possui um maior nível de detalhamento na modelagem do sistema, considerando

inclusive as restrições advindas da operação elétrica do sistema de forma detalhada. A

divisão do planejamento em médio, curto e curtíssimo prazo deve-se à dificuldade de

manuseio simultâneo da estocasticidade das vazões afluentes com a representação

individualizada das usinas hidrelétricas.

Baseando-se na incerteza das vazões afluentes (Figura 5), foi proposta a atual

decomposição do problema do planejamento da operação em modelos diferenciados

Figura 5: Incerteza nas vazões x horizonte de planejamento. Figura 2.10: Variação da incerteza ao longo do tempo (BERTHO, 2013).

O objetivo básico do planejamento da operação é obter, para cada etapa, as

metas de geração de cada usina (hidro e termoelétrica) do sistema de forma a atender a

carga e minimizar o valor esperado do custo de operação ao longo do período de planeja-

mento. Este custo é composto pelo custo variável de combustível das usinas termoelétricas

e pelo custo atribuído às interrupções de fornecimento de energia, representado por uma

função de penalização dos dé�cits de energia (custo do dé�cit).

O modelo NEWAVE foi desenvolvido para o planejamento de médio prazo

(até 5 anos). É um modelo de planejamento da operação de sistemas hidrotérmicos

com representação individualizada do parque termoelétrico e representação agregada, por

meio de reservatórios equivalentes de energia, do parque hidroelétrico. Um dos principais

resultados obtidos na estratégia de solução do modelo NEWAVE é a função de custo

32

futuro. É por meio dessa função que o encadeamento com o modelo de curto prazo

(DECOMP) é realizado.

O modelo DECOMP realiza a otimização do despacho energético no horizonte

de curto prazo (até 12 meses), para as usinas do sistema de forma individualizada. Com

base nas curvas de custo futuro obtidas na etapa de médio prazo, determinam-se as

metas individuais de geração das usinas hidráulicas e térmicas do sistema, bem como os

intercâmbios de energia entre subsistemas.

Finalmente na programação diária, de�ne-se uma programação horária de gera-

ção que atenda as metas estabelecidas na etapa anterior, sujeita às condições operacionais

da rede elétrica (CEPEL, 2012a).

A seguir serão apresentados de forma sucinta os modelos NEWAVE e DE-

COMP.

2.4.1 NEWAVE

O modelo NEWAVE foi desenvolvido para tratar o planejamento da operação

energética do SIN (Sistema Integrado Nacional) quando são considerados horizontes de

planejamento de médio prazo, ou seja, horizontes plurianuais. Nesse tipo de horizonte, a

estocasticidade das vazões a�uentes passa a ser um fator de grande in�uência, pois imerge

o problema em um ambiente de grande incerteza (BERTHO, 2013).

O SIN é representado por quatro subsistemas, mostrados na Figura 2.11, onde

a usina de Itaipu é considerada à parte por ser dividida com o Paraguai e a região Norte

Isolado não participa do planejamento da operação do SIN (BERTHO, 2013). A existência

de interligações com os sistemas vizinhos permite uma redução dos custos de operação,

por meio do intercâmbio de energia e um aumento da con�abilidade de fornecimento, via

repartição das reservas de energia (CEPEL, 2012a).

Figura 2.11: Divisão dos subsistemas brasileiros (BERTHO, 2013).

33

A maneira encontrada pelo setor elétrico para introduzir a estocasticidade

intrínseca ao problema do planejamento de médio prazo na modelagem computacional foi

gerar séries sintéticas de energias, como mostrado na Figura 2.12, para depois empregar

técnicas de programação dinâmica dual estocástica (PDDE).

Mod

elod

eO

timizacao/S

imu

lacao

Dados Historicos

ModeloEstocastico Serie

Sintetica 1

SerieSintetica 2

SerieSintetica 3

Vh

V1

V2

V3

Figura 2.12: Esquema da geração de séries sintéticas de energia.

A técnica da PDDE é utilizada para determinar a energia que será gerada

pelas usinas hidroelétricas e pelas usinas térmicas e os intercâmbios entre os subsistemas,

atendendo à carga de energia no país num horizonte de cinco anos. A função objetivo

do modelo é minimizar o valor esperado do custo operacional do sistema, composto por

duas parcelas, uma que re�ete o custo presente do combustível, operação e manutenção

das usinas térmicas e outra que re�ete o custo futuro do não atendimento da demanda de

energia em função da decisão de gerar diferentes parcelas de energia hidráulica e térmica

num determinado estágio (chamada de Função de Custo Futuro) (LOPES, 2007).

Objetivando a viabilidade do modelo computacional, para que se evite a cha-

mada �maldição da dimensionalidade�, uma importante simpli�cação do SIN é realizada

no NEWAVE, onde os reservatórios das usinas hidroelétricas são agregados em reservató-

rios equivalentes de energia, esquematizado pela Figura 2.13, reduzindo-se a dimensão do

sistema interligado brasileiro.

34

Energia Afluente

Afluencia Controlavel

Reservatorio Equivalentede Energia

Perdas Evaporativas

Energia Hidroeletrica

AfluenciaEstocastica

Figura 2.13: Esquema de um reservatório equivalente.

Com o processo de agregação dos reservatórios, surge a necessidade de realizar

uma adequação dos valores da vazão a�uente nas usinas hidroelétricas para a compatibi-

lização dessa abordagem com as grandezas veri�cadas no sistema. Dessa forma, a vazão

a�uente deve deixar de ser abordada como quantidade de água que chega ao sistema,

mas sim como quantidade de energia, conhecida como energia natural a�uente (ENA). O

processo de obtenção da ENA é realizado calculando-se o aproveitamento energético da

vazão a�uente na usina hidroelétrica e em todas as usinas à jusante (BERTHO, 2013).

As usinas térmicas são representadas no NEWAVE por grupos de usinas com

custos semelhantes, chamados de classes térmicas. As classes térmicas são de�nidas ba-

sicamente por meio dos valores de geração máxima e mínima das usinas térmicas, assim

como seus custos incrementais de operação. O dé�cit de energia é uma ocorrência de

grande relevância, e foi abrangido no modelo como uma unidade térmica de capacidade

igual à demanda de energia.

Entre os resultados do NEWAVE estão: os despachos de geração térmica e

hidráulica nos subsistemas equivalentes, as energias fornecidas e recebidas entre os subsis-

temas e a Função Custo Futuro, esta última é usada como entrada de um segundo modelo,

o chamado DECOMP (LOPES, 2007). Além disso, o modelo NEWAVE é rodado uma

vez por mês, enquanto o DECOMP é rodado toda semana.

2.4.2 DECOMP

Na cadeia de modelos utilizada no Brasil para a realização do planejamento

da operação energética do sistema de potência, a alocação ótima da geração de energia

elétrica no horizonte de curto prazo é desempenhada pelo DECOMP. O modelo DECOMP

35

formula o problema de despacho hidrotérmico como um problema de programação linear

no qual a função objetivo consiste em minimizar o custo de operação ao longo do período

de planejamento, dado o conjunto de informações disponíveis (carga, vazões, disponibi-

lidades, limites de transmissão entre subsistemas, função de custo futuro do NEWAVE).

Assim, para o planejamento de curto prazo condizer com o planejamento de médio prazo, é

necessário um acoplamento entre os modelos NEWAVE e DECOMP, que é esquematizado

na Figura 2.14.

Centro de Pesquisas de Energia Elétrica - CEPEL 24 de Julho de 2013

Preço do Mercado de Curto Prazo ou

Preço de Liquidação de Diferenças (PLD)

Baseia-se no cálculo do CMO e é determinado semanalmente para cada submercado (Norte, Nordeste, Sudeste/Centro-Oeste e Sul) e para cada patamar de carga (pesada, média e leve), sem considerar as usinas em teste e as restrições de transmissão internas a cada submercado, incorporadas pelo ONS no planejamento da operação

O PLD é limitado ainda por um preço máximo e um mínimo, estabelecidos pela ANEEL

O modelo NEWAVE é rodado uma vez ao mês, enquanto o modelo DECOMP é rodado toda semana

Figura 2.14: Entradas e saídas de dados nos modelos NEWAVE e DECOMP (CPAMP,2013b).

Diferentemente do modelo NEWAVE, o DECOMP apresenta um maior deta-

lhamento do sistema em sua modelagem matemática, passando a representar as usinas

hidroelétricas de forma individualizada. Isso possibilita a inserção de uma série de novas

características no modelo de otimização, como enchimento de volume morto, volume de

espera, evaporação, desvios de água, tempo de viagem da vazão de�uente, entre outras

características inerentes a operação das usinas hidroelétricas.

A estocasticidade das vazões é introduzida no modelo DECOMP a partir do

�nal do primeiro mês do planejamento, por meio da construção de cenários de vazões men-

sais a�uentes às usinas do sistema. Estes cenários hidrológicos podem ser representados

por uma árvore de a�uências, cuja teoria é abordada no capítulo 3.

36

Assim, o resultado da otimização realizada pelo DECOMP de curto prazo,

indicará para cada usina, semanalmente no primeiro mês do horizonte de planejamento e

mensalmente nos demais, a alocação ótima da geração de energia respeitando, na medida

do possível, as restrições operativas de cada usina.

37

3 REPRESENTAÇÃO DOS CENÁRIOS HIDROLÓ-

GICOS

Dado o objetivo da dissertação de veri�cação da viabilidade de substituição

da árvore de a�uências por latisse, este capítulo irá apresentar o modelo de geração de

séries sintéticas e vazões GEVAZP (CEPEL, 2010), responsável pela geração dos cenários

hidrológicos necessários ao modelo NEWAVE e DECOMP, o conceito e a formulação da

teoria de árvores e o conceito e a formulação adaptada da teoria de latisse para a possível

substituição.

3.1 GEVAZP

O modelo GEVAZP é o responsável pela geração dos cenários hidrológicos

necessários ao modelo NEWAVE e DECOMP, realizando a representação dos possíveis

cenários de vazões de forma diferenciada para cada etapa do processo do plajenamento

da operação.

No modelo NEWAVE, os cenários de a�uência podem ser representados em

uma estrutura paralela (pente), como mostrado na Figura 3.1, ou em uma estrutura de

árvore.

Figura 3.1: Séries em paralelo (pente) (CEPEL, 2010).

Já no modelo DECOMP, os cenários hidrológicos gerados a partir do modelo

GEVAZP são representados apenas pela estrutura de árvore de a�uências, como mostrada

na Figura 3.2, cuja teoria será abordada na seção a seguir.

38

Figura 3.2: Séries em árvore (CEPEL, 2010).

A estrutura de pente exige um menor esforço computacional para sua resolu-

ção, porém também agrega menos informações, enquanto a estrutura de árvore apresenta

informações mais densas, porém de mais complexa resolução. Dessa forma, propõe-se

aplicar também nos modelos NEWAVE e DECOMP a estrutura de latisse, que apresenta

um porte intermediário entre as estruturas já presentes no modelo GEVAZP. Para que a

estrutura de latisse seja utilizada, não seria necessário modi�car os modelos NEWAVE e

DECOMP, e sim o modelo GEVAZP, de modo que fosse possível optar para que o modelo

pudesse gerar além das estruturas já pré-existentes, também uma série com a estrutura

de latisse, como mostrado na Figura 3.3.

Modelo

GEVAZP

Modelo

NEWAVE

Modelo

DECOMP

Seriesem pente

Seriesem arvore

Seriesem latisse

Dadoshistoricos

Figura 3.3: Representação dos cenários hidrológicos através do modelo GEVAZP.

Para que essa alteração possa ser aplicada, será apresentado nas seções à seguir

a teoria de árvore e latisse.

39

3.2 Representação por Árvore

Utilizando os valores históricos de energia natural a�uente, a incerteza da

a�uência futura pode ser representada por uma distribuição de probabilidades, que, por

sua vez, pode ser representada por meio de uma árvores de a�uências, em que cada ramo

possui uma probabilidade de ocorrência relacionada à seu respectivo cenário representado,

de modo que a soma das probabilidades de cada ramo em cada estágio deve ser igual à 1,

como mostrado na Figura 3.4.

Figura 3.4: Distribuição de probabilidades dos valores históricos de energia natural a�u-ente de Janeiro para a bacia do Tocantins e sua respectiva árvore.

De forma geral, para árvores binomiais, calibra-se os ramos de forma a obter-se

dois cenários equiprováveis, sendo que um cenário representa a ocorrência de um valor

médio menor que o esperado, e o outro cenário representa a ocorrência de um valor médio

maior que o esperado. Tem-se assim que quanto maior a diferença entre os valores repre-

sentados por cada ramo, maior é a dispersão entre os pontos sob a curva de distribuição.

Para resolver o problema completo do planejamento de operação deve-se resol-

ver um problema para cada nó da árvore de a�uências. Cada nó da árvore representa uma

decisão que deve ser tomada a cada estágio, e cada decisão conduz através de seus ramos

a diferentes resultados. Se a decisão for guardar a água no presente e vierem a ocorrer

vazões elevadas no futuro, possivelmente será preciso verter o excedente, o que representa

desperdício da energia que foi guardada. Por outro lado, se a decisão for usar a água no

presente e vierem a ocorrer vazões reduzidas no futuro, possivelmente será preciso elevar

a complementação termoelétrica utilizando unidades mais dispendiosas, ou mesmo cortar

a demanda por meio de políticas de racionamento.

A árvore de a�uências é montada de forma que cada nó represente o volume

40

de água mantido no reservatório, e cada ramo represente um possível cenário de a�uência.

Desse modo, para cada ramo há um valor de geração associado e, para cada nó, há um

valor de custo de geração esperado.

O volume da água representado em cada nó é calculado por etapas forward,

através das equações de balanço hídrico, que consideram que o volume ao �nal de um

estágio é equivalente ao volume inicial, somado à a�uência do ramo que conduz à este

nó, descontado a quantidade de água usada na geração hidroelétrica neste estágio e o

vertimento ocorrido.

O valor esperado de cada nó é calculado por etapas backward, e equivale à

esperança (ou média ponderada, de acordo com a probabilidade de ocorrência) do custo

de geração dos estágios seguintes.

3.2.1 Árvore Binária

Exempli�cando para o caso mais simples, como o da árvore binária apre-

sentada na Figura 3.5, aplica-se a Equação (2.1) para os cálculos de balanço hídrico,

desconsiderando-se a quantidade de água in�ltrada e evaporada, por se tratarem de uma

parcela pequena comparadas ao volume total, obtendo-se as Equações (3.1) à (3.6), sendo

V ti o volume do reservatório no �nal do estágio t e ao inicio do estágio (t + 1) no nó

i, Atij, GHtij e I

tij a a�uência, a geração hidroelétrica e o vertimento do cenário ij do

estágio (t) que representa a transição do nó i do estágio (t− 1) para o nó j do estágio t,

respectivamente.

V 00

V 11

V 21

V 12

V 22

V 23

V 24

A102

A101

A211

A212

A224

A223

GH101

GH102

GH211

GH212

GH223

GH224

Figura 3.5: Esquema de cálculo do balanço hídrico na árvore binária.

V 11 = V 0

0 + A101 −GH1

01 − I101 (3.1)

V 12 = V 0

0 + A102 −GH1

02 − I102 (3.2)

41

V 21 = V 1

1 + A211 −GH2

11 − I211 (3.3)

V 22 = V 1

1 + A212 −GH2

12 − I212 (3.4)

V 23 = V 1

2 + A223 −GH2

23 − I223 (3.5)

V 24 = V 1

2 + A224 −GH2

24 − I224 (3.6)

Tendo calculado a geração hidroelétrica de cada estágio, tem-se um custo de

geração associado à cada ramo devido a geração por meio de termoelétricas para a com-

plementação da energia para o atendimento da demanda. Assim, considerando que os

ramos da árvore de a�uências apresentam probabilidade P e (1− P ) de ocorrência dado

a ocorrência do nó anterior, como mostrado na Figura 3.6, tem-se que o cálculo do valor

esperado do custo de operação Sti do estágio t e nó i, pode ser representado pelas Equações

(3.7), (3.8) e (3.9).

S00

S11

S21

S12

S22

S23

S24

(1− P )

P

P

(1− P )

(1− P )

P

Figura 3.6: Esquema de cálculo do valor esperado do custo de operação na árvore binária.

S11 = P · S2

1 + (1− P ) · S22 (3.7)

S12 = P · S2

3 + (1− P ) · S24 (3.8)

S00 = P · S1

1 + (1− P ) · S12 (3.9)

3.2.2 Árvore Ternária

Caso a árvore fosse ternária ao invés de binária, como mostrado na Figura

3.7, e cada um dos seus ramos tivesse a probabilidade P1, P2 e P3 de ocorrência dado a

ocorrência do nó anterior, teriam-se as equações Eqs. (3.10) à (3.13) para o cálculo do

valor esperado do custo de operação de cada nó.

42

S00

P1

P2

P3

P1

P1

P1

P2

P2

P2

P3

P3

P3

S11

S12

S13

S21

S22

S23

S24

S25

S26

S27

S28

S29

Figura 3.7: Esquema de cálculo do valor esperado do custo de operação na árvore ternária.

S11 = P1 · S2

1 + P2 · S22 + P3 · S2

3 (3.10)

S12 = P1 · S2

4 + P2 · S25 + P3 · S2

6 (3.11)

S13 = P1 · S2

7 + P2 · S28 + P3 · S2

9 (3.12)

S00 = P1 · S1

1 + P2 · S12 + P3 · S1

3 (3.13)

De forma análoga, aos cálculos feitos para a árvore binária, poderia-se realizar

os cálculos de balanço hídrico para a árvore ternária, obtendo-se doze equações para a

árvore ternária de dois estágios (uma para cada ramo da árvore), onde o volume �nal

do reservatório em cada estágio seria equivalente ao volume inicial somado à a�uência e

subtraído o gasto de água para a geração e o vertimento.

3.2.3 Árvore Multiária

Por �m, seguindo o mesmo raciocínio, generalizando para uma árvore multiá-

ria, tem-se a Figura 3.8, cujo o cálculo do balanço hídrico de cada ramo é mostrado pela

Equação (3.14), onde Atij, GHtij, e I

tij são a a�uência, a água utilizada para a geração

hidroelétrica e o vertimento de água do estágio t que levam do volume V (t−1)i do nó i do

estágio (t− 1) até o volume V tj do nó j do estágio t. Havendo assim, para os dois estágios

da árvore, um total de (n + n2) equações, onde n é o número de cenários, ou seja, uma

equação para cada ramo.

43

V 00

GH101

GH102

GH211

V 11

V 12

V 21

GH10n

GH21n

GH22(2n)

GH2n(n2)

GH22(n+1)

GH2n(2n−2)

V 22

V 2n

V 2(n+1)

V 2(n+2)

V 2(2n)

V 1n

V 2(n2)

V 2(n2−1)

V 2(n2−2)

Figura 3.8: Esquema de cálculo do balanço hídrico na árvore multiária.

V tj = V

(t−1)i + Atij −GH t

ij − I tij (3.14)

O cálculo do valor esperado do custo de operação de cada nó está representado

pela Figura 3.9 e Equação (3.15), havendo uma equação para cada nó, ou seja, (n + 1)

equações para os dois estágios da árvore.

S00

P1

P2

P1S11

S12

S21

Pn

Pn

Pn

Pn

P1

P1

S22

S2n

S2n+1

S2n+2

S22n

S1n

S2n2

S2n2−n

S2n2−n+1

Figura 3.9: Esquema de cálculo do valor esperado do custo de operação na árvore mul-tiária.

44

Sti = P1 · S(t+1)((i−1)·n+1) + P2 · S(t+1)

((i−1)·n+2) + ...+ Pn · S(t+1)(i·n) (3.15)

3.3 Representação por Latisse

A representação dos cenários de a�uência por latisse segue os mesmos princí-

pios da representação por árvore, porém é preciso �forçar� para que os nós intermediários

coincidam, de modo que os valores assumidos por estes ainda sejam contínuos.

Isto é possível, pois o problema de otimização do despacho hidrotérmico mo-

delado no DECOMP é considerado do tipo �acaso-decisão� (CEPEL, 2012a), ou seja,

sabendo-se primeiro a a�uência ocorrida, toma-se a decisão do despacho hidrotérmico de

forma que os nós intermediários coincidam.

Além disso, de forma geral, as usinas hidroelétricas brasileiras apresentam

um reservatório de grande capacidade de armazenamento, sendo possível ter um maior

controle de decisão do quanto estas usinas irão gerar em cada estágio. Sendo assim, as

decisões de cada nó podem caminhar através de diferentes ramos de forma a �amarrar� os

nós intermediários da árvore de a�uências, fazendo com que esta passe a ter a estrutura de

uma latisse. E portanto, o número de possíveis cenários cresceria de forma mais moderada.

Os cálculos para o balanço hídrico e valor esperado do custo de operação, por

sua vez, também seriam menores em termos quantitativo, e são mostrados nas subseções

a seguir.

3.3.1 Latisse Binomial

Exempli�cando para um caso mais simples, como o da latisse binomial apre-

sentada na Figura 3.10, continua-se aplicando a Equação (2.1) para os cálculos de balanço

hídrico, desconsiderando-se também a quantidade de água in�ltrada e evaporada, por se

tratarem de uma parcela pequena comparadas ao volume total, e obtendo-se de forma

simpli�cada as Equações (3.16) à (3.21), sendo V ti o volume do reservatório no �nal do

estágio t e ao inicio do estágio (t + 1) no nó i, Atij, GHtij e I

tij a a�uência, a geração

hidroelétrica e o vertimento do cenário ij do estágio t, que representa a transição do nó i

do estágio (t− 1) para o nó j do estágio t, respectivamente.

45

A211

A102

V 00

V 11

V 12

V 12

V 22

V 23

A101

A212

A222

A223

GH101

GH102

GH211

GH212

GH222

GH223

Figura 3.10: Esquema de cálculo do balanço hídrico na latisse binomial.

V 11 = V 0

0 + A101 −GH1

01 − I101 (3.16)

V 12 = V 0

0 + A102 −GH1

02 − I102 (3.17)

V 21 = V 1

1 + A211 −GH2

11 − I211 (3.18)

V22 = V 1

1 + A212 −GH2

12 − I212 (3.19)

V22 = V 1

2 + A222 −GH2

22 − I222 (3.20)

V 23 = V 1

2 + A223 −GH2

23 − I223 (3.21)

É importante chamar a atenção para os cálculos que envolvem os nós interme-

diários, como as Equações (3.19) e (3.20), que devem coincidir, uma vez que representam

os cálculos do volume de um mesmo nó, a razão para se obter duas equações é devido ao

fato de ser possível chegar ao mesmo nó através de caminhos diferentes.

Tendo calculado a geração hidroelétrica de cada estágio, tem-se um custo de

operação associado à cada ramo devido a geração por meio de termoelétricas para a

complementação da energia para o atendimento da demanda. Assim, considerando que

os ramos da árvore de a�uências apresentam probabilidade P e (1 − P ) de ocorrência

dado a ocorrência do nó anterior, como mostrado na Figura 3.11, tem-se que o cálculo do

valor esperado do custo de operação Sti do estágio t e nó i, pode ser representado pelas

Equações (3.22), (3.23) e (3.24).

46

(1− P )

S0

S11

S12

S21

S22

S23

P

P

(1− P )

(1− P )

P

Figura 3.11: Esquema de cálculo do valor esperado du custo de geração na latisse binomial.

S11 = P · S2

1 + (1− P ) · S22 (3.22)

S12 = P · S2

2 + (1− P ) · S23 (3.23)

S00 = P · S1

1 + (1− P ) · S12 (3.24)

Neste caso, o que diferencia os cálculos da latisse para a árvore seria o fato

do valor S22 ser considerado duas vezes, dado que é possível chegar neste nó através de

dois caminhos diferentes e, dessa forma, o nó intermediário apresenta uma probabilidade

maior de ocorrência que os demais nós �nais.

Pode-se observar ainda que as equações, tanto de balanço hídrico como de

cálculo do valor esperado do custo de operação, se aparentam muito com as equações

utilizadas para a árvore binária, porém, há um menor número de nós �nais, o que acar-

reta em um menor número de equações a serem calculadas caso mais um estágio fosse

adicionado ao período analisado.

3.3.2 Latisse Trinomial

Para o caso de uma latisse trinomial como a mostrada na Fig. 3.12, com

probabilidades P1, P2 e P3 de ocorrência de cado ramo dado a ocorrência do nó anterior,

utiliza-se as Eqs. (3.25) à (3.28) para o cálculo do valor esperado do custo de cada nó.

47

S00

S11

S13

S21

S22

S25

P1

P3

P2S12

S24

S23

P1

P1

P1

P2

P2

P2

P3

P3

P3

Figura 3.12: Esquema de cálculo do valor esperado do custo de operação na latisse trino-mial.

S11 = P1 · S2

1 + P2 · S22 + P3 · S2

3 (3.25)

S12 = P1 · S2

2 + P2 · S23 + P3 · S2

4 (3.26)

S13 = P1 · S2

3 + P2 · S24 + P3 · S2

5 (3.27)

S00 = P1 · S1

1 + P2 · S12 + P3 · S1

3 (3.28)

De forma análoga, aos cálculos feitos para a latisse binomial, poderia-se realizar

os cálculos de balanço hídrico para a latisse trinomial, obtendo-se três equações para o

primeiro estágio e nove equações para o segundo estágio, sendo uma equação para o

cálculo de V 21 , duas equações para o cálculo de V 2

2 , três equações para o cálculo de V 23 ,

duas equações para o cálculo de V 24 e uma para o cálculo de V 2

5 , devido ao número de

possíveis caminhos através dos ramos que levariam à estes nós.

Até o segundo estágio, o número de equações da latisse e da árvore são os

mesmos, mas, a partir do terceiro estágio, a latisse começa a apresentar menos equações

que a árvore devido ao menor número de nós e ramos.

3.3.3 Latisse Multinomial

Por �m, generalizando para o caso de uma latisse multinomial como mostrada

na Figura 3.13, tem-se que a equação para o cálculo do balanço hídrico de cada ramo

segue o mesmo padrão dos demais casos, cujo volume inicial (V tj ) do cenário ij do estágio

t é equivalente ao volume inicial V (t−1)i do nó i do estágio (t− 1) somado à a�uência Atij

e subtraído à geração hidroelétrica GH tij e o vertimento I

tij do ramo que levou do volume

de um estágio ao do outro, como mostrado na Equação (3.14) para árvore multinomial.

48

V0

V 11

GH101

GH102

GH103

GH10n

V 12

V 13

V 1n

GH211GH2

12

GH213

V 21

V 22

V 23

V 2n

V 2(2n−1)

V 2(n+1)

V 2(n+2)

GH2n(2n−1)

Figura 3.13: Esquema de cálculo do balanço hídrico na latisse multinomial.

O que difere o cálculo das equações de balanço hídrico da árvore para a latisse,

é que na latisse pode-se chegar a um mesmo nó por mais de um caminho através dos

ramos da latisse, portanto, haverá mais de uma equação de balanço hídrico para o cálculo

do volume dos nós intermediários. Assim, o número de equações de balanço hídrico para

cada nó varia entre 1 e o número de ramos n, conforme a posição i de cada nó.

Mesmo havendo mais de uma equação de balanço hídrico para cada nó, ainda

assim haverá apenas uma equação de balanço hídrico para cada ramo, e assim, a latisse

apresentará um menor número de equações que a árvore, uma vez que esta apresenta um

menor aumento do número de ramos a cada estágio.

Por �m, para o cálculo do valor esperado do custo de operação Sti de cada nó,

tem-se a Figura 3.14 e a Equação (3.29), onde t representa o estágio, i a posição do nó e

n o número de cenários. Havendo uma equação de cálculo do valor esperado do custo de

operação para cada nó da estrutura.

49

S0

S11

P1 P2

P3

Pn

S12

S13

S1n

P1

P1

P1

P1

P2

P2

P2

P3

P3

Pn

Pn

Pn

Pn

S21

S22

S23

S2n

S2(2n−1)

S2(n+1)

S2(n+2)

Figura 3.14: Esquema de cálculo do valor esperado do custo de operação na latisse mul-tinomial.

Sti = P1 · S(t+1)i + P2 · S(t+1)

(i+1) + ...+ Pn · S(t+1)(i+n−1) (3.29)

Para um maior aprofundamento sobre a teoria de árvore e latisse, o livro do

Luenberger (1998) pode ser consultado.

50

4 MODELO DINÂMICO DE OTIMIZAÇÃO COM IN-

DIFERENÇA AO RISCO

Tendo apresentado os modelos NEWAVE e DECOMP e a representação da

estocasticidades das vazões a�uentes por árvore e latisse, este capítulo objetiva apresentar

os modelos dinâmicos de otimização com indiferença ao risco. Para que haja um maior

entendimento dos modelos, antes será apresentado a teoria de Programação Dinâmica

Dual Estocástica.

4.1 Programação Dinâmica Dual Estocástica

No modelo NEWAVE, o problema de planejamento da operação energética de

médio prazo é representado por um problema de programação estocástica linear multi-

estágio, baseado na aproximação do valor esperado da função de custo futuro. A técnica

utilizada para encontrar a solução ótima é a Programação Dinâmica Dual Estocástica

(PDDE), baseada na decomposição de Benders. A formulação da PDDE atualmente

empregada no modelo leva em conta a correlação temporal das a�uências aos reservatórios.

Assim, consideram-se como variáveis de estados do problema o armazenamento no início

do período e as a�uências passadas.

O cálculo da política de operação do NEWAVE é feito iterativamente, sendo

que cada iteração é constituída por duas etapas. Uma delas é a etapa forward, em que é

estimado o custo médio imediato, dado pela geração térmica mais custo de não suprimento

de energia. A etapa backward, que é uma etapa recursiva, calcula o custo de operação do

último estágio e atualiza o valor presente do estágio anterior, até que se chegue ao estágio

inicial. Nesta etapa é calculada a função de custo futuro por meio dos cortes de Benders

(DEUS, 2010).

A teoria de Programação Dinâmica Dual Estocástica pode ser vista com mais

detalhes no Manual de Referência do Modelo NEWAVE (CEPEL, 2012b) e a Teoria de

Corte de Benders associada pode ser estudada no livro Optimization Theory for Large

Systems (LASDON, 1970).

De forma sucinta, o algoritmo de PDDE pode ser ilustrado por um problema

de dois estágios (Figura 4.1) apresentado pelo problema de otimização restrito (4.1), no

qual o primeiro estágio é considerado conhecido e o segundo estágio depende dos valores

que uma ou mais variáveis aleatórias podem assumir.

51

10 estagio 20 estagio

P1

P2

Figura 4.1: Representação de uma árvore binomial de dois estágios.

y = min(c1x101 + P1c2xt11 + P2c

2xt12)

s.a.

Q1x101≥ b11R1x101 +Q2x211≥ b21R1x101 +Q2x212≥ b22

(4.1)

Este problema pode ser dividido em dois subproblemas, um do primeiro estágio

e outro do segundo. O subproblema do primeiro estágio está representado no conjunto de

Equações (4.2), onde c1x01 representa o custo imediato, e a função α1(x101) representa o

valor esperado do custo de operação futuro, que é consequência da decisão x1∗01.

y = min(c1x101 + α1(x

101)

s.a.

Q1x101≥ b11

(4.2)

O subproblema do segundo estágio está representado no conjunto de Equações

(4.3), que depende da solução x1∗01 encontrada na resolução do primeiro estágio para ser

resolvido.

α1(x

101) = min(P1c

2x211 + P2c2x212)

s.a.

Q2x211≥ b21 −R1x1∗01

Q2x212≥ b22 −R1x1∗01

(4.3)

A solução do primeiro estágio é utilizada no segundo estágio através de uma

etapa forward, caso o valor da função objetivo encontrado no segundo estágio não coin-

cida com o valor de α1(x101) encontrado no primeiro estágio, deve-se realizar uma etapa

52

backward, adicionando ao subproblema do primeiro estágio a equação de corte de Ben-

ders associada, mostrada na Equação (4.4), onde π1 e π2 são os multiplicadores simplex

associados às restrições do problema.

P1π1(b21 −R1x101) + P2π2(b

22 −R1x101)≤ α (4.4)

Dessa forma, deve-se realizar iterações entre as etapas forward e backward até

que os valores convirjam e a solução do problema seja encontrada.

4.2 Modelo Integrado com Árvore Binária

Para a resolução do problema de otimização do despacho hidrotérmico de forma

integrada para uma árvore binária com N estágios, sendo o primeiro estágio considerado

conhecido, como mostrado na Figura 4.2, é necessário de�nir a função objetivo e as

equações de restrições.

10 estagio 20 estagio 30 estagio N-agesimo estagio

V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V 34

V N2(N−1)

V N1

V N2

V N2(N−1)−1

Figura 4.2: Representação de uma árvore binária multi-estágios.

A função objetivo é de�nida como o mínimo da soma do custo esperado de

operação de cada estágio, como mostrado na Equação (4.5), onde c representa o vetor

linha de custos de geração por unidade de energia de cada usina, que, nesse caso é con-

siderado constante no tempo, xtij e ptij representam o vetor coluna da geração de energia

de cada usina k e a probabilidade de ocorrência de cada cenário ij de um estágio t,

respectivamente. Os cálculos dos valores de ptij podem ser vistos no apêndice A.

53

y = min(cx101 + (p211cx211 + p212cx

212) + (p311cx

311 + p312cx

312 + p323cx

323 + p324cx

324) + ...

+(pN11cxN11 + pN12cx

N12 + ...pN(2(N−2))(2(N−1))cx

N(2(N−2))(2(N−1)))) (4.5)

As restrições às quais a função objetivo está sujeita são: as restrições de potên-

cia, atendimento da demanda e balanço hídrico, de acordo com os parâmetros utilizados

para a de�nição da árvore de a�uências.

As restrições de potência determinam que a geração de cada usina não pode

ultrapassar sua respectiva potência máxima, havendo assim uma equação para cada usina

k em cada cenário ij de cada estágio t, que é representada de forma genérica pela Equação

(4.6), onde X tijk representa os elementos do vetor xtij, que, por sua vez, representa a

quantidade de energia gerada por cada usina k e Gk representa a potência máxima de

geração de cada usina k.

X tijk≤ Gk (4.6)

As restrições de atendimento da demanda determinam que o total de energia

gerado X tijk por todas as usinas k em um cenário ij deve atender à demanda Dt do período

t, havendo uma equação por ramo da árvore, como representado pela Equação (4.7), onde

K representa o número total de usinas.

K∑k=1

X tijk = Dt (4.7)

Por �m, as restrições de balanço hídrico determinam, através de um balanço de

massa, o volume do reservatório em cada cenário de cada estágio, sendo representado pelo

conjunto de Equações (4.8), havendo também uma equação para cada ramo da árvore,

onde V tj representa o volume, Atij a a�uência, I tij o vertimento do cenário ij do estágio t

e GH tij representa o elemento do vetor xtij referente a quantidade de energia gerada pela

usina hidroelétrica.

54

V 11 = V 0

0 + A101 −GH1

01 − I101V 21 = V 1

1 + A211 −GH2

11 − I211V 22 = V 1

1 + A212 −GH2

12 − I212V 31 = V 2

1 + A311 −GH3

11 − I311V 32 = V 2

1 + A312 −GH3

12 − I312V 33 = V 2

2 + A323 −GH3

23 − I323V 34 = V 2

2 + A324 −GH3

24 − I324...

V N2(N−1)−1 = V

(N−1)2(N−2) + AN

(2(N−2))(2(N−1)−1) +GHN(2(N−2))(2(N−1)−1) + IN

(2(N−2))(2(N−1)−1)

V N2(N−1) = V

(N−1)2(N−2) + AN

(2(N−2))(2(N−1))+GHN

(2(N−2))(2(N−1))+ IN

(2(N−2))(2(N−1))

(4.8)

4.3 Modelo Decomposto com Árvore Binária


decomposta é necessário de�nir a estrutura dos subproblemas a serem resolvidos, de forma

que haja o encadeamento entre os problemas para que seja encontrada a solução do

problema global. A árvore binária mostrada na Figura 4.3 apresenta as divisões dos

subproblemas.


V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V 34

V N2(N−1)

V N1

V N2

V N2(N−1)−1

Figura 4.3: Representação dos subproblemas de uma árvore binomial multi-estágios.

Assim, para o primeiro subproblema deve-se de�nir a função objetivo y como

o mínimo da soma do custo esperado de geração do primeiro estágio cx101 e do custo

esperado de geração dos estágio futuros α2, como mostrado na Equação (4.9).

55

y = min(cx101 + α2) (4.9)

Quanto às restrições de potência, atendimento da demanda e balanço hídrico

que a função objetivo está submetida, estas devem ser referentes apenas ao primeiro

estágio, sendo mostradas no conjunto de Equações (4.10).

X1

01k≤ Gk∑Kk=1X

101k = D1

V 11 = V 0

0 + A101 − A1

01 − A101

(4.10)

Para o subproblema do segundo estágio, a função objetivo α2 é de�nida como

o mínimo da soma do custo esperado de geração do segundo estágio e do custo esperado

de geração dos estágio futuros α31 e α32, como mostrado na Equação (4.11). Neste caso

há dois termos referentes à geração dos estágios futuros devido ao fato de haver dois

subproblemas no terceiro estágio.

α2 = min((p211cx211 + p212cx

212) + α31 + α32) (4.11)

Quanto às restrições de potência, atendimento da demanda e balanço hídrico

que a função objetivo está submetida, estas devem ser referentes apenas ao segundo

estágio, sendo mostradas no conjunto de Equações (4.12), onde V 1∗1 representa o valor

do volume obtido por meio da solução do subproblema do primeiro estágio e j pode

assumir os valores 1 ou 2.

X21jk≤ Gk∑Kk=1X

21jk = D2

V 11 = V 1∗

1

V 21 = V 1

1 + A211 −GH2

11 − I211V 22 = V 1

1 + A212 −GH2

12 − I212

(4.12)

Por raciocínio análogo ao anterior, pode-se fazer esse mesmo processo sucessi-

vamente até os subproblemas do N-agésimo estágio, que são divididos em 2(N−2) subpro-

blemas, sendo o primeiro representado pela função objetivo da Equação (4.13) sujeito às

restrições do conjunto de Equações (4.14), com j assumindo os valores 1 ou 2.

56

αN1 = min(pN11cxN11 + pN12cx

N12) (4.13)

XN1jk≤ Gk∑Kk=1X

N1jk = DN

V(N−1)1 = V

(N−1)∗1

V N1 = V

(N−1)1 + AN11 −GHN

11 − IN11V N2 = V

(N−1)1 + AN12 −GHN

12 − IN12

(4.14)

O último subproblema representado pela função objetivo da Equação (4.15)

sujeito às restrições do conjunto de Equações (4.16), com j representando (2(N−1)− 1) ou

2(N−1). Note que a função objetivo destes subproblemas é apenas a minimização do valor

esperado do custo de geração do estágio N , uma vez que não há estágios futuros.

αN(2(N−2)) = min(pN(2(N−2))(2(N−1)−1)cxN(2(N−2))(2(N−1)−1) + pN(2(N−2))(2(N−1))cx

N(2(N−2))(2(N−1)))

(4.15)

XN(2(n−2))jk

≤ Gk∑Kk=1X

N(2(n−2))jk

= DN

V(N−1)2(N−2) = V

(N−1)∗2(N−2)

V N2(N−1)−1 = V

(N−1)2(N−2) + AN

(2(N−2))(2(N−1)−1) +GHN(2(N−2))(2(N−1)−1) + IN

(2(N−2))(2(N−1)−1)

V N2(N−1) = V

(N−1)2(N−2) + AN

(2(N−2))(2(N−1))+GHN

(2(N−2))(2(N−1))+ IN

(2(N−2))(2(N−1))

(4.16)

Desse modo, resolvendo-se o subproblema do primeiro estágio, um valor de

V 1∗1 será obtido. Este valor deve ser passado para o segundo estágio, que, ao ser resolvido,

encontrará outros valores de volumes �nais que deverão ser passados para os subproble-

mas dos estágio seguintes e assim por diante, até se chegar no último subproblema (etapa

forward). Feita esta primeira etapa, inicia-se a etapa backward, onde a partir da reso-

lução do último subproblema acrescenta-se uma equação de corte, Equação (4.17), no

subproblema anterior.

α(t+1)j − α∗(t+1)j≥ π(t+1)j · (V (t−1)j − V (t−1)∗

j ) (4.17)

57

Sendo π(t+1)j a média dos valores de πtij, que se refere ao valor que se econo-

mizaria caso houvesse uma unidade a mais de água no cenário ij do estágio t.

Estas etapas são repetidas até a solução convergir.

4.4 Modelo Integrado com Latisse Binomial


integrada para uma latisse binomial com N estágios, sendo o primeiro estágio considerado

conhecido, como mostrado na Fig. 4.4, é necessário de�nir a função objetivo e as equações

de restrições.


V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V NN

V N1

V N2

Figura 4.4: Representação de uma latisse binomial multi-estágios.

A função objetivo é de�nida como o mínimo custo esperado de operação, como

mostrado na Equação (4.18), onde c representa o vetor linha de custos de geração por

unidade de energia de cada usina, que nesse caso é considerado constante no tempo, xtije ptij representam o vetor coluna da geração de energia de cada usina k e a probabilidade

de ocorrência de cada cenário ij de cada estágio t, respectivamente. Sendo os cálculos

dos valores de ptij apresentados no apêndice A.

y = min(cx101 + (p211cx211 + p212cx

212) + (p311cx

311 + p312cx

312 + p322cx

322 + p323cx

323) + ...


N12 + ...+ pN(N−1)Ncx

N(N−1)N)) (4.18)

As restrições às quais a função objetivo está sujeita são: as restrições de potên-

cia, atendimento da demanda e balanço hídrico, de acordo com os parâmetros utilizados

58

para a de�nição da latisse binomial.

Assim como para a árvore, as restrições de potência determinam que a geração

de cada usina não pode ultrapassar sua respectiva potência máxima, havendo assim uma

equação para cada usina k em cada cenário ij de cada estágio t, que é representada de

forma genérica pela Equação (4.6), diferindo apenas nos subíndices devido à mudança na

estrutura.

As restrições de atendimento da demanda determinam que o total de energia

gerado X tijk por todas as usinas k em um cenário ij deve atender à demanda Dt do

período t, havendo uma equação por ramo da latisse, sendo representado de forma genérica

também pela Equação (4.7), diferindo mais uma vez apenas nos subíndices ij.

Por �m, as restrições de balanço hídrico determinam, através de um balanço de

massa, o volume do reservatório em cada cenário de cada estágio, sendo representado pelo

conjunto de Equações (4.19), havendo também uma equação para cada ramo da latisse,

onde V tj representa o volume, Atij a a�uência, I tij o vertimento do cenário ij do estágio t

e GH tij representa o elemento do vetor x

tij, referente a quantidade de energia gerada pela

usina hidroelétrica.

V 11 = V 0

0 + A101 −GH1

01 − I101V 21 = V 1

1 + A211 −GH2

11 − I211V 22 = V 1

1 + A212 −GH2

12 − I212V 31 = V 2

11 + A311 −GH3

11 − I311V 32 = V 2

12 + A212 −GH2

12 − I212V 32 = V 2

22 + A322 −GH3

22 − I322V 33 = V 2

23 + A323 −GH3

23 − I323...

V N(N−1) = V

(N−1)(N−1) + AN(N−1)(N−1) +GHN

(N−1)(N−1) + IN(N−1)(N−1)V NN = V

(N−1)(N−1) + AN(N−1)N) +GHN

(N−1)N + IN(N−1)N

(4.19)

É importante ressaltar que, para os volumes intermediários como o V 32 , haverá

duas equações para o cálculo do mesmo volume, uma vez que é possível chegar no mesmo

estado (valor) de volume por dois caminhos diferentes.

Então, comparando-se o modelo integrado com árvore binária com o modelo

integrado com latisse binomial, pode-se notar semelhança entre suas equações, que se

alteram somente de acordo com a estrutura dos cenários de a�uência. Assim, para a

latisse binomial, haverá uma menor quantidade de restrições, uma vez que estas variam

de acordo com o número de cenários e a latisse apresenta uma menor quantidade de

cenários que a árvore.

59

4.5 Modelo Decomposto com Latisse Binomial


decomposta é necessário de�nir a estrutura dos subproblemas a serem resolvidos, de forma

que haja o encadeamento entre os problemas, para que �nalmente seja encontrada a

solução do problema global. Assim, a latisse binomial representada na Figura 4.5 mostra

as divisões dos subproblemas.


V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V NN

V N1

V N2

Figura 4.5: Representação dos subproblemas de uma latisse binomial multi-estágios.

Devido ao fato dos nós intermediários da latisse coincidirem, cada estágio deve

ser resolvido de forma única. Como para os dois primeiros estágios a estrutura da latisse

binomial é a mesma da árvore binária, o primeiro subproblema da latisse é representado

pela função objetivo mostrada na Equação (4.9) e as restrições do conjunto de Equações

(4.10), assim como para a árvore.

O segundo subproblema é representado pela função objetivo da Equação (4.20),

que difere da Equação (4.11) apenas no termo referente ao valor esperado do custo das

gerações futuras, devido à diferença entre árvore e latisse no número de subproblemas

no estágio seguinte, estando sujeita às restrições do mesmo conjunto de Equações (4.12)

da árvore de a�uências, uma vez que até este estágio, as duas estruturas apresentam a

mesma forma.

α2 = min(p211cx211 + p212cx

212 + α3) (4.20)

60

O subproblema do terceiro estágio é representado pela função objetivo mos-

trada na Equação (4.21) sujeita às restrições do conjuntos de Equações (4.22), cuja função

objetivo apresenta o mínimo do valor esperado do custo de geração do terceiro estágio

mais o valor esperado do custo de geração dos estágios seguintes, sujeito às restrições de

potência, atendimento da demanda e balanço hídrico do terceiro estágio, onde V 2∗1 e V 2∗

2

são obtidos por meio da otimização do subproblema do segundo estágio, e i e j podem

assumir os valores 1 ou 2 e 1,2 ou 3, respectivamente, de acordo com os ramos da latisse.

α3 = min(p311cx311 + p312cx

312 + p322cx

322 + p323cx

323 + α4) (4.21)

X3ijk≤ Gk∑Kk=1X

3ijk = D3

V 21 = V 2∗

1

V 22 = V 2∗

2

V 31 = V 2

1 + A311 −GH3

11 − I311V 32 = V 2

1 + A312 −GH3

12 − I312V 32 = V 2

2 + A322 −GH3

22 − I322V 33 = V 2

2 + A323 −GH3

23 − I323

(4.22)

Por �m, o subproblema do N-agésimo estágio apresenta a função objetivo como

o mínimo do valor esperado do custo de operação do N-agésimo estágio apenas, uma vez

que não há estágios futuros, como mostrado na Equação (4.23).

αN = min(pN11cxN11 + pN12cx

N12 + ...+ pN(N−1)Ncx

N(N−1)N) (4.23)

Esta equação está sujeita às restrições de potência (uma para cada usina em

cada cenário), atendimento da demanda (uma para cada cenário), valores de volumes

iniciais obtidos do subproblema anterior (uma para cada nó inicial) e balanço hídrico

(uma para cada ramo), apresentadas de forma genérica no conjunto de Equações (4.24),

onde i pode assumir valores entre 1 e (N − 1) e j pode assumir os valores entre 1 e N ,

sendo a combinação ij feita de acordo com os ramos da latisse.

61

XNijk≤ Gk∑Kk=1X

Nijk = DN

V(N−1)1 = V

(N−1)∗1

V(N−1)2 = V

(N−1)∗2

...

V(N−1)(N−1) = V

(N−1)∗(N−1)

V N1 = V

(N−1)1 + AN11 −GHN

11 − IN11V N2 = V

(N−1)1 + AN12 −GHN

12 − IN12...

V N(N−1) = V

(N−1)(N−1) + AN(N−1)(N−1) +GHN

(N−1)(N−1) + IN(N−1)(N−1)V NN = V

(N−1)(N−1) + AN(N−1)N) +GHN

(N−1)N + IN(N−1)N

(4.24)

As iterações são realizadas da mesma forma que no método decomposto utili-

zando árvore binária, havendo uma equação de corte por cada nó �nal de cada subpro-

blema a cada iteração.

62

5 MODELO DINÂMICO DE OTIMIZAÇÃO COMAVER-

SÃO AO RISCO

No capítulo 4 foi apresentada a teoria de Programação Dinâmica Dual Esto-

cástica e os modelos integrados e decompostos utilizando árvore binária e latisse binomial.

Este capítulo tem como objetivo adicionar aos modelos o uso de uma métrica de risco

nomeada Conditional Value at Risk. Para tanto, será apresentada a seguir a teoria desta

ferramenta e os modelos aversos ao risco.

5.1 Conditional Value at Risk

O objetivo relevante dos mecanismos de aversão a risco (MARs) é encontrar

uma solução de compromisso entre o aumento da segurança e os impactos nos custos do

sistema. Em sistemas hidrotérmicos, estes mecanismos buscam antecipar o despacho de

geração térmica com custos variáveis unitários de operação (CVUs) mais baixos, com o

intuito de evitar o atingimento, no futuro, de níveis indesejáveis de armazenamento nos

reservatórios das usinas hidroelétricas e, com isso, minimizar o risco de dé�cits de energia,

mas sem onerar em demasia os custos de operação do sistema (CPAMP, 2013b).

O CVaR é uma medida de risco apropriada à otimização, visto que possui as

propriedades de monotonicidade, invariância sobre translações, homogeneidade positiva

e subaditividade, com consequente convexidade; tendo sido implementada recentemente

nos modelos NEWAVE e DECOMP, como forma de reduzir o risco na consideração dos

cenários de a�uência.

Esta medida tem a função de representar uma parcela que agrega maiores

informações quando da ocorrência de valores extremos, fazendo com que seja possível

a representação de não somente o valor esperado ou central, mas também de efeitos

relacionados à ocorrência de valores extremos das distribuições.

A teoria geral de CVaR é apresentada a seguir, e foi feita com base no artigo

de Shapiro et al. (2012) e na apresentação da CCEE (2012).

Suponha uma distribuição de probabilidade de retornos esperados r = f(u,w)

que assume valores de 0 a 1, sendo αβ, Value at Risk (VAR), é o menor valor para um

nível de con�ança β , como mostrado na Figura 5.1.

63

17

representação de não somente o valor esperado ou central, mas também de efeitos

relacionados à ocorrência de valores extremos das distribuições.

A teoria geral de CVaR é apresentada a seguir, e foi feita com base numa

apresentação da CCEE, 2012.

Imagine uma distribuição de probabilidade de retornos esperados que

assume valores de 0 a 1, onde , Value at Risk (VAR), é o menor valor para um nível de

confiança , como mostrado nas figura 2.7.

Figura 2.7 – Distribuição de Retornos [CCEE, 2012].

Matematicamente, temos que:

Equação 2.1

Equação 2.2

O CVaR ( ) de uma distribuição de retornos esperados é representado na figura

2.8, e mede a média das piores perdas esperadas para um determinado nível de confiança

, podendo ser descrito matematicamente como:

Equação 2.3

Figura 5.1: Distribuição de retornos (CCEE, 2012).

Matematicamente, tem-se que:

β =

∫ ∞α

r(u,w)dw = Pr[r≥ αβ] (5.1)

1− β =

∫ α

−∞r(u,w)dw = Pr[r≤ αβ] (5.2)

O CVaR (φβ) de uma distribuição de retornos esperados é representado na

Figura 5.2, e mede a média das piores perdas esperadas para um determinado nível de

con�ança β, podendo ser descrito matematicamente como:

18

Figura 2.8 – Representação do CVaR numa curva de distribuição de retornos

[CCEE, 2012].

Quando se trata de uma distribuição de perdas, é o maior valor para um nível de

confiança , portanto, se localiza na cauda direita da distribuição como mostrado na figura

2.9.

Figura 2.9 – Distribuição de Perdas [CCEE,2012].

Assim, se a função considerada representar as perdas esperadas de um portfólio, o

sinal da equação para o cálculo do CVaR se inverte e a equação passa a ser escrita como:

Equação 2.4

Dessa forma, em vez de minimizar o valor esperado de operação, será minimizada a

combinação convexa entre o valor esperado e o valor assumido pela medida de risco de

Figura 5.2: Representação do CVaR numa curva de distribuição de retornos (CCEE,2012).

φβ(x) = E[r|r≤ αβ(u)] =∫∞−∞ r · Pr[w|r≤ αβ(u)]dw (5.3)

Quando se trata de uma distribuição de perdas, αβ é o maior valor para um

nível de con�ança β. Portanto, se localiza na cauda direita da distribuição, como mostrado

na Figura 5.3.

64

18

Figura 2.8 – Representação do CVaR numa curva de distribuição de retornos

[CCEE, 2012].

Quando se trata de uma distribuição de perdas, é o maior valor para um nível de

confiança , portanto, se localiza na cauda direita da distribuição como mostrado na figura

2.9.

Figura 2.9 – Distribuição de Perdas [CCEE,2012].

Assim, se a função considerada representar as perdas esperadas de um portfólio, o

sinal da equação para o cálculo do CVaR se inverte e a equação passa a ser escrita como:

Equação 2.4

Dessa forma, em vez de minimizar o valor esperado de operação, será minimizada a

combinação convexa entre o valor esperado e o valor assumido pela medida de risco de

Figura 5.3: Distribuição de perdas (CCEE, 2012).

Assim, se a função considerada representar as perdas esperadas de um portfó-

lio, o sinal da equação para o cálculo do CVaR se inverte e a equação passa a ser escrita

como:

φβ(x) = E[r|r≥ αβ(u)] =∫∞−∞ z · Pr[w|r≥ αβ(u)]dw (5.4)

Dessa forma, em vez de minimizar o valor esperado de operação, será mini-

mizada a combinação convexa entre o valor esperado e o valor assumido pela medida de

risco, de acordo com um peso a ser considerado (MARCATO et al., 2013).

Assim, de modo análogo, pode-se considerar o uso do CVaR no problema

de otimização do despacho hidrotérmico de duas formas: considerando a ocorrência de

menores a�uências que acarretariam em maiores gastos devido a um maior despacho

térmico, ou considerando diretamente a ocorrência de maiores valores esperados do custo

de geração, como mostram as Figuras. 5.4 e 5.5, respectivamente.

Figura 5.4: Distribuição de a�uências.

65

Figura 5.5: Distribuição de custos esperados.

Desse modo, com a aplicação do CVaR, a função objetivo dos modelos de

otimização do despacho hidrotérmico deixa de minimizar apenas o valor esperado do custo

total de operação (CTO), e passa a minimizar o CTO com um peso (1−λ) somado à uma

parcela referente ao custo dos cenários hidrológicos mais críticos (CVaR), identi�cados

por meio de um parâmetro α (relacionado ao nível de proteção), com um peso λ. Dessa

forma, a função objetivo é dada pela Equação (5.5).

min((1− λ)E[CTO] + λCV aRα) (5.5)

A determinação dos valores dos parâmetros λ e α estão associadas ao maior

ou menor grau de aversão ao risco que se deseja adotar. A política de operação se torna

tanto mais avessa ao risco quando λ se aproximar de 1, uma vez que o CVaR apresentará

uma ponderação maior no cálculo da média, e quanto mais o percentual α se aproximar

de zero, uma vez que será considerado cenários mais críticos para o cálculo do CVaR.

Atualmente os valores dos parâmetros λ e α adotados são 25% e 50%, respec-

tivamente (CPAMP, 2013a).

5.2 Modelo Integrado com Árvore Binária

Como visto anteriormente, para a resolução do problema de otimização do

despacho hidrotérmico de forma integrada é necessário de�nir a função objetivo como o

mínimo custo esperado e inserir as equações de restrições de potência, atendimento da

demanda e balanço hídrico de acordo com os parâmetros utilizados para a de�nição da

árvore de a�uências. A mudança agora inserida se situa na alteração da função objetivo,

que passa a minimizar não apenas o custo esperado da operação, mas também uma parcela

adicional do custo dos cenários mais críticos.

Assim, para a árvore binária apresentada no capítulo 4, considerando-se o

uso do CVaR com parâmetro α = 50%, por exemplo, tem-se os cenários mais críticos

66

apresentados dentro do retângulo na Figura 5.6.


V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V 34

V N2(N−1)

V N1

V N2

V N2(N−1)−1

Figura 5.6: Representação dos cenários mais críticos para uma árvore binária multi-estágios.

O modelo utilizado para a resolução do problema averso ao risco é apresentado

pela função objetivo mostrada na (5.6), que é de�nida pela minimização do custo do pri-

meiro estágio, que é determinístico, somado à ponderação do custo esperado de operação

dos demais estágios e à ponderação do valor esperado do custo de geração dos cenários

mais críticos, estando sujeita as mesmas equações de restrição do modelo indiferente ao

risco (Equações (4.6), (4.7) e (4.8)).

y = min(cx101 + (1− λ)((p211cx211 + p212cx

212) + (p311cx

311 + p312cx

312 + p323cx

323 + p324cx

324) + ...


N12 + ...pN(2(N−2))(2(N−1))cx

N(2(N−2))(2(N−1)))) +

+λ(cx212 + (pa323cx23 + pa324cx324) + ...

+(paN(2(N−3)+1)(2(N−2)+1)cxN(2(N−3)+1)(2(N−2)+1) + ...+ paN(2(N−2))(2(N−1))cx

N(2(N−2))(2(N−1)))))

(5.6)

Sendo c o vetor linha de custos de geração por unidade de energia de cada usina,

que nesse caso é considerado constante no tempo, xtij e ptij o vetor coluna da geração de

energia de cada usina k e a probabilidade de ocorrência de cada cenário ij de cada estágio

t, respectivamente, e patij a probabilidade auxiliar considerando a ocorrência apenas dos

cenários mais secos.

67

5.3 Modelo Decomposto com Árvore Binária


decomposta considerando a aversão ao risco utiliza-se a mesma divisão dos subproblemas

mostrada na capítulo 4, devendo-se novamente mudar apenas a função objetivo, para que

o CVaR seja considerado, como mostrado na Figura 5.7.


V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V 34

V N2(N−1)

V N1

V N2

V N2(N−1)−1

Figura 5.7: Representação dos subproblemas de uma árvore binomial multi-estágios aversaao risco.

O primeiro subproblema mantém-se igual ao caso indiferente ao risco represen-

tado pela função objetivo da Equação (4.9), sujeito às restrições do conjunto de Equações

(4.10), uma vez que o primeiro estágio é considerado conhecido.

O segundo subproblema é representado pela função objetivo de�nida como

a soma ponderada do valor esperado do custo de operação com o custo de ocorrência

dos cenários mais secos, acrescido dos custos esperados dos subproblemas dos estágios

seguintes, como representado na Equação (5.7), estando sujeito às mesmas restrições do

caso indiferente ao risco representadas no conjunto de equações (4.12).

α2 = min((1− λ)(p211cx211 + p212cx

212) + λ(cx212) + α31 + α32) (5.7)

Para o N-agésimo estágio também haverá 2(N−2) subproblemas, sendo o pri-

meiro representado pela função objetivo da Equação (5.8), que considera apenas a pon-

deração por (1 − λ), uma vez que este não apresenta nenhum cenário considerado parte

dos 50% pessimistas e nem estágios futuros, estando sujeito às restrições do conjunto de

Equações (4.14)

68

αN1 = min((1− λ)(pN11cxN11 + p12

NcxN12)) (5.8)

Por �m, o último subproblema do N-agésimo estágio apresenta a função ob-

jetivo como a minimização do valor esperado do custo de geração do N-agésimo estágio,

ponderado com um fator (1 − λ) somado ao valor esperado do custo de geração dos ce-

nários mais pessimistas ponderado com um fator λ, como mostrado na Equação (5.9),

estando sujeito às restrições apresentadas no conjunto de Equações (4.24).

αN(2(N−2)) = min((1− λ)(pN(2(N−2))(2(N−1)−1)cxN(2(N−2))(2(N−1)−1) + pN(2(N−2))(2(N−1))cx

N(2(N−2))(2(N−1))) +

λ(paN(2(N−2))(2(N−1)−1)cxN(2(N−2))(2(N−1)−1) + paN(2(N−2))(2(N−1))cx

N(2(N−2))(2(N−1))))

(5.9)

As iterações são realizadas do mesmo modo que foi mostrado no capítulo 4.

5.4 Modelo Integrado com Latisse Binomial

Assim como no modelo integrado com árvore binária, a diferença entre o mo-

delo averso ou não ao risco se encontra na função objetivo, que passa a dar uma maior

importância para a ocorrência dos cenários mais críticos. Desse modo, utilizando nova-

mente o parâmetro α = 50%, tem-se os cenários mais críticos representados dentro do

retângulo na Figura 5.8.


V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V NN

V N1

V N2

Figura 5.8: Representação dos cenários mais críticos para uma latisse binomial multi-estágios.

69

O modelo utilizado na resolução do problema é representado pela função ob-

jetivo da Equação (5.10), que é de�nida pela minimização do custo do primeiro estágio,

que é determinístico, somado à ponderação do custo esperado de operação dos demais

estágios e à ponderação do valor esperado do custo de geração dos cenários mais críticos,

estando sujeita as mesmas equações de restrição do modelo indiferente ao risco (Equações

(4.6), (4.7) e (4.19)).

y = min(cx101 + (1− λ)((p211cx211 + p212cx

212) + (p311cx

311 + p312cx

312 + p322cx

322 + p323cx

323) + ...


N12 + ...+ pN(N−1)Ncx

N(N−1)N)) +

+λ(cx212 + (pa322cx322 + pa323cx

323) + ...+ (paNmMcx

NmM + ...+ paN(N−1)Ncx

N(N−1)N)))

(5.10)

Sendo c o vetor linha de custos de geração por unidade de energia de cada

usina, que nesse caso é considerado constante no tempo, xtij e ptij representam o vetor

coluna da geração de energia de cada usina k e a probabilidade de ocorrência de cada

cenário ij de cada estágio t, respectivamente, e patij representa a probabilidade auxiliar

considerando a ocorrência apenas dos cenários mais secos.

Onde os parâmetros m e M são calculados de acordo com as Equações condi-

cionais (5.11) e (5.12).

m =

{N−12

+ 1, se N for ímparN−12

+ 12, se N for par

(5.11)

M =

{N2

+ 1, se N for parN2

+ 12

se N for impar(5.12)

5.5 Modelo Decomposto com Latisse Binomial

Novamente, para a resolução do problema de otimização do despacho hidro-

térmico de forma decomposta considerando a aversão ao risco, utiliza-se a mesma divisão

dos subproblemas mostrada no capítulo 4, devendo-se novamente mudar apenas a função

objetivo, para que o CVaR seja considerado, como mostrado na Figura 5.9.

70


V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V NN

V N1

V N2

Figura 5.9: Representação dos subproblemas de uma latisse binomial multi-estágios aversaao risco.

O primeiro subproblema mantém-se igual ao caso indiferente ao risco tanto

para árvore quanto para latisse, representado pela função objetivo da Equação (4.9) e as

restrições do conjunto de Equações (4.10), uma vez que o primeiro estágio é considerado

conhecido.

O segundo subproblema é representado pela função objetivo (5.13), que difere

da Equação (5.7) apenas no termo referente ao valor esperado do custo das gerações

futuras, devido à diferença entre árvore e latisse no número de subproblemas no estágio

seguinte, estando sujeita às restrições do mesmo conjunto de Equações (4.12) da árvore

de a�uências.

α2 = min((1− λ)(p211cx211 + p212cx

212) + λ(cx212) + α3) (5.13)

O terceiro subproblema é representado pela função objetivo da Equação (5.14),

sujeito às restrições do conjunto de Equações (4.22).

α3 = min((1− λ)(p311cx311 + p312cx

312 + p322cx

322 + p323cx

323) +

+λ(pa322cx322 + pa323cx

323) + α4) (5.14)

Por �m, o subproblema do N-agésimo estágio apresenta a função objetivo

como a minimização da soma do valor esperado do custo de geração do N-agésimo estágio

com o valor esperado do custo de geração dos cenários mais secos apenas, uma vez que

71

não há estágios futuros, como mostrado na Equação (5.15), estando sujeito às restrições

apresentadas no conjunto de Equações (4.24).

αN = min((1− λ)(pN11cxN11 + pN12cx

N12 + ...+ pN(N−1)Ncx

N(N−1)N) +

+λ(paNmMcxNmM + ...+ paN(N−1)Ncx

N(N−1)N)) (5.15)

As iterações são também realizadas do mesmo modo que foi mostrado no ca-

pítulo 4.

72

6 ESTUDO DE CASO

O estudo de caso apresentado neste capítulo tem o objetivo de veri�car de

forma prática a viabilidade da substituição da árvore de a�uências pela latisse. Para isso,

foi proposto um problema de otimização do despacho hidrotérmico utilizando os dados

da bacia do Tocantins, mostrados na seção 6.1. A partir destes dados, o problema foi

resolvido utilizando primeiramente a árvore de a�uências e depois utilizando a latisse,

tanto para o caso neutro ao risco, quanto para o caso averso ao risco, variando o dado do

volume inicial do reservatório da hidroelétrica equivalente de modo a se obter uma melhor

análise dos resultados. Assim, os resultados obtidos são apresentados na seção 6.2.

6.1 Metodologia e Dados do Problema

Os dados utilizados neste estudo de caso foram obtidos através da ONS (2014)

e da CCEE (2014), sendo todos os valores convertidos para MWmed, considerado-se a

base temporal de um mês, de forma a evitar problemas de incompatibilidade de unidades

durante a otimização.

Para a geração hidroelétrica, foram consideradas as usinas hidroelétricas da

bacia do Tocantins, ilustradas na Figura 6.1 na página seguinte, que foram agregadas

numa usina hidroelétrica equivalente (usina 1 do problema), de forma a reduzir a dimensão

do problema, obtendo como restrição de geração máxima o valor de 12.821,6MWmed.

O volume do reservatório considerado foi obtido no site da ONS (2016), que

mostra, que para o mês de abril de 2015, os reservatórios do subsistema norte apresenta-

vam um volume equivalente à 12.006MWmed, correspondente à 81,06% do volume total;

dessa forma, o reservatório equivalente apresenta como restrição o volume máximo cor-

respondendo à 14.811,3MWmed. O volume inicial considerado para o reservatório foi

inicialmente de 10.000MWmed e depois alterado para 9.000MWmed, sendo que o volume

�nal não poderia ser inferior à 4.000MWmed.

Para a geração termoelétrica, foram consideradas as usinas termoelétricas pró-

ximas à bacia do Tocantins, adicionada de uma usina �ctícia, que teve função de represen-

tar o possível intercâmbio de energia entre os subsistemas. As gerações máximas e custos

de cada usina foram obtidos através da tabela de consolidação dos leilões, disponível no

site da CCEE (2014), e são mostrados na Tabela 6.1 também na página seguinte, onde as

termoelétricas Geramar I e Geramar II foram consideradas como uma única usina.

73

Figura 6.1: Usinas hidroelétricas na bacia do Tocantins e suas respectivas capacidadesmáximas de geração (ONS, 2014).

Tabela 6.1: Custo de geração das termoelétricas.

Usina Potência(MW) Tipo de Combustível Preço (R$/MWh)Maranhão III 499.2 Gás 127,04

Termomaranhão 350.0 Carvão 198,60Geramar I e II 330.0 Óleo 211,40Intercâmbio 700.0 �� 300,00

Para o cálculo da demanda a ser atendida, foi considerada a demanda do

subsistema norte e seu intercâmbio com o subsistema nordeste entre os anos de 2000 e

74

2014 (dados disponíveis no histórico de operação do ONS). A média obtida para cada mês

é apresentada na Tabela 6.2.

Tabela 6.2: Dados da demanda de energia.

MêsSubsistema Norte

(MWh/h)Intercâmbio N-NE

(MWmed)Total

(MWmed)Janeiro 3763,9 3905,7 7669,5Fevereiro 3813,6 3961,2 7774,8Março 3846,5 3996,7 7843,2Abril 3890,4 4046,8 7937,2Maio 3892,9 4044,1 7937,0Junho 3885,9 4037,5 7923,4Julho 3885,8 4061,0 7946,8Agosto 3983,0 4162,9 8145,9Setembro 3999,0 4193,6 8192,6Outubro 3973,4 4165,3 8138,7Novembro 3970,0 4162,9 8132,8Dezembro 3989,9 4179,7 8169,7

Por �m, para a inserção da estocasticidade das vazões no problema, foram

considerados os dados do histórico de ENAs (energia natural a�uente) entre os anos de

1931 à 2013 para a bacia do Tocantins. A partir destes dados, foram feitas, para cada

mês, a distribuição de probabilidade acumulada e foram selecionados valores médios para

um cenário otimista e para um cenário pessimista de vazão, de forma que cada cenário

apresenta 50% de chance de ocorrer. Estes valores estão apresentados na Tabela 6.3.

Tabela 6.3: Cenários de a�uência.

MêsCenário Pessimista

(MWmed)A�uência Média

(MWmed)Cenário Otimista

(MWmed)Janeiro 10495,7 12765,0 16676,3Fevereiro 13416,4 16672,3 21907,0Março 15528,1 18511,7 25145,1Abril 14313,3 17254,7 21255,8Maio 8526,0 10676,1 13080,9Junho 4534,5 5265,8 6598,0Julho 2934,7 3404,5 4000,5Agosto 2118,3 2463,2 2885,7Setembro 1732,2 2089,4 2509,9Outubro 2027,2 2639,2 3335,3Novembro 3534,1 4197,5 5722,1Dezembro 6304,1 8028,6 10643,0

75

Assim, o problema a ser resolvido neste estudo de caso consiste em determinar

a geração de cada usina de modo à minimizar os custos de geração, atender à demanda

e respeitar as restrições de balanço hídrico e potência, como mostrado resumidamente no

esquema da Figura 6.2

Demanda

(Dt)

c1 = 0

G1 = 12821, 65 MWmed

G2 = 499, 2 MWmed

G3 = 350 MWmed

G4 = 330 MWmed

G5 = 350 MWmed

c2 = 127, 04 R$/MWmed

c3 = 198, 60 R$/MWmed

c4 = 211, 40 R$/MWmed

c5 = 300, 00 R$/MWmed

Afluencia

Figura 6.2: Esquema do problema a ser resolvido.

Como o número de cenários de uma árvore de a�uências cresce exponencial-

mente de acordo com o número de estágios analisados, di�cultando a resolução do pro-

blema para períodos muito longos, a análise apresentada nesta dissertação irá abranger

apenas 4 estágios (meses de maio, junho, julho e agosto, respectivamente), sendo a escolha

dos meses feita com base na transição do período úmido para o seco.

A resolução do problema foi realizada pela aplicação da árvore de a�uências e

após, por meio da aplicação de latisse binomial, tanto do modo único quanto decomposto.

Finalmente, as soluções serão comparadas e será aplicado o CVaR com os parâmetros

α = 50% e λ = 25%.

A Figura 6.3 apresenta a árvore de a�uências a ser analisada, cujo primeiro

estágio é determinativo, em que cada retângulo representa um subproblemas a ser soluci-

onado pelo método de PDDE. Para a resolução por latisse binomial, a representação dos

subproblemas a serem resolvidos é mostrado na Figura 6.4.

76

1 2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

10 estagio 30 estagio 40 estagio20 estagio

V 00 V 1

1

V 21

V 22

V 31

V 32

V 33

V 34

V 41

V 42

V 43

V 44

V 45

V 46

V 47

V 48

1

2

3A

3B

4A

4B

4C

4D

Figura 6.3: Esquema da árvore de a�uências utilizada no problema.

4o estagio3o estagio2o estagio1o estagio

1 2V 21

V 22

V 31

V 32

V 33

V 41

V 42

V 44

V 43

1

2

3

4

V 00 V 1

1

3

4

5

6

7

8

9

10

11

Figura 6.4: Esquema da latisse binomial utilizada no problema.

6.2 Resolução e resultados

A resolução deste problema foi feita com base na teoria apresentada nos capí-

tulos anteriores com o auxílio do software Lingo 11.0.

Os resultados obtidos para o primeiro caso, em que o volume inicial do reserva-

tório considerado foi de 10.000MWmed, coincidiram, tanto pela solução integrada quanto

pela solução decomposta, utilizando árvore binária ou latisse binomial, e são mostrados

nas tabelas do apêndice B, sendo o valor da função objetivo para o caso indiferente ao

risco de R$ 638.781,20 e para o caso averso ao risco de R$ 692.508,00. As árvores com

77

os valores do custo de geração de cada ramo e o valor de cada volume do reservatório

obtidos durante a otimização estão mostradas nas Figuras 6.5 e 6.6.

1 2

3

4

5

6

7

8

9

10

11

12

13

14

15

16


10000 13837.5

13011.3

11627.8

9260.2

8848.4

8848.4

8148.4

4000

4000

4000

Volume do Reservatorio (MWmed)

Custo do Ramo (R$)

LEGENDA

4000

4000

4000

4000

4000

185609.2

63418.37

202690.4

24798.21

132928.4

200090.1

308740.4

0

116682.9

52315.07

202690.4

52315.07

202690.4

188442.0

412690.4

Figura 6.5: Solução para o caso indiferente ao risco utilizando árvore.

1 2

3

4

5

6

7

8

9

10

11

12

13

14

15

16


10000 13918.3

13092.1

11708.6

9260.2

8848.4

8848.4

8148.4

4000

4000

4000


Custo do Ramo (R$)

LEGENDA

4000

4000

4000

4000

4000

202690.4

63418.37

202690.4

14533.38

116881.5

183009.0

284500.4

0

116682.9

52315.07

202690.4

52315.07

202690.4

188442.0

412690.4

Figura 6.6: Solução para o caso averso ao risco utilizando árvore.

78

As latisses, são apresentadas nas Figuras 6.7 e 6.8.


1 213011.3

10627.8

9260.2

8848.4

8148.4

4000

4000

4000

4000

1000013837.5

3

4

5

6

7

8

9

10

11


185609.2

Custo do Ramo (R$)

63418.37

202690.4

24798.21

132928.4

200090.1

308740.4

0

116682.9

52315.07

202690.4

188442.0

412690.4

LEGENDA

Figura 6.7: Solução para o caso indiferente ao risco utilizando latisse.


1 213092.1

10708.6

9260.2

8848.4

8148.4

4000

4000

4000

4000

10000 13918.3

3

4

5

6

7

8

9

10

11


202690.4

Custo do Ramo (R$)

63418.37

202690.4

14533.38

116881.5

183009.0

284500.4

0

116682.9

52315.07

202690.4

188442.0

412690.4

LEGENDA

Figura 6.8: Solução para o caso averso ao risco utilizando latisse.

79

Analisando-se os resultados, a árvore de a�uências apresentou para os nós

intermediários os mesmos valores de volume, de modo que, plotando-se num grá�co o

resultado �nal dos volumes do reservatório em cada estágio, obtém-se uma sobreposição

dos nós da árvore com a latisse, como mostrado na Figura 6.9.

Volume(MWmed)

Estagio1 2 3 40

8000

9000

10000

11000

12000

13000

14000

Arvore

Latisse

4000

Figura 6.9: Grá�co comparativo entre os volumes do reservatório da árvore e da latissepara o caso neutro ao risco.

Então, comparando-se os valores de volume encontrados para o caso neutro ao

risco e averso ao risco, tanto por árvore quanto por latisse, o caso averso ao risco apresentou

valores mais elevados no reservatório, como mostrado na Figura 6.10 e, consequentemente,

maiores custos de geração, por necessitar de uma maior geração termoelétrica, uma vez

que há uma maior retenção de água.

Volume(MWmed)

Estagio1 2 3 40

4000

9000

10000

11000

12000

13000

14000

Volume neutro ao risco

Volume averso ao risco

8000

Figura 6.10: Grá�co comparativo entre os volumes do reservatório para o caso neutro eaverso ao risco.

80

Agora, para o segundo caso, em que o volume inicial do reservatório conside-

rado foi de 9.000MWmed, os resultados são mostrados nas Figuras 6.11 e 6.12 e mais

detalhadamente no apêndice B.

1 2

3

4

5

6

7

8

9

10

11

12

13

14

15

16


9000 12918.3

12442.1

11281.3

8995

8609.2

8514.2

8148.4

4000

4000

4000


Custo do Ramo (R$)

LEGENDA

4000

4000

4000

4000

4000

202690.4

132928.4

374500.4

63418.37

202690,4

202690.4

412690.4

33691.01

171699.1

93565.85

274450.4

112432.8

302950.4

188442.0

412690.4

Figura 6.11: Solução para o caso indiferente ao risco utilizando árvore, com volume inicialde 9.000MWmed.


1 212347.1

11281.3

8900

8514.2

8148.4

4000

4000

4000

4000

900012918.3

3

4

5

6

7

8

9

10

11


202690.4

Custo do Ramo (R$)

114061.4

374500.4

63418.37

202690.4

202690.4

412690.4

45759.81

191782.1

112432.8

302950.4

188442.0

412690.4

LEGENDA

Figura 6.12: Solução para o caso indiferente ao risco utilizando latisse, com volume inicialde 9.000MWmed.

81

Um grá�co comparativo entre os volumes do reservatório em cada estágio da

árvore e da latisse é mostrado na Figura 6.13.

Volume(MWmed)

Estagio1 2 3 40

8000

9000

10000

11000

12000

13000

14000

Arvore

Latisse

4000

Figura 6.13: Grá�co comparativo entre os volumes do reservatório da árvore e da latissepara o caso neutro ao risco, com volume inicial de 5.000MWmed.

Como pode-se observar, os valores obtidos não coincidiram. Para a resolução

neutra ao risco através do uso da árvore binária, o valor obtido pela função objetivo do

custo esperado da operação foi de R$ 875.517,30. Enquanto para a resolução através de

latisse binomial, o valor obtido do custo esperado da operação foi de R$ 876.023,70.

Dessa forma, constatou-se que pode haver casos em que os resultados obtidos

poderão ser diferentes. Porém, ao que tudo indica, mesmo os resultados não coincidindo,

este se apresentam bem próximos.

Assim, pelo estudo de caso aqui apresentado, a substituição da árvore de

a�uências pela latisse mostrou-se aparentemente viável, tanto para o caso neutro ao risco

quanto para caso averso ao risco, mesmo não havendo uma completa coincidência dos

valores apresentados.

82

7 CONCLUSÕES

Os programas de otimização desenvolvidos conseguiram alcançar o objetivo

proposto, isto é, são capazes de determinar a quantidade de energia que cada usina deverá

gerar, minimizando o custo esperado, atendendo a demanda, considerando a estocastici-

dade das a�uências, e respeitando as restrições de cada usina (geração máxima e volume

máximo) e de balanço hídrico.

Comparando-se a resolução através da programação única e da PDDE, tem-se

que os valores obtidos por ambos os programas coincidiram, como o esperado. Sendo

que na resolução utilizando PDDE o número de iterações �nais de cada programa foi

bem menor em relação ao número de iterações necessárias para a resolução do problema

com programação única. Desse modo, tem-se que apesar de mais complexo a resolução

utilizando a técnica de PDDE, esta requer um menor esforço computacional, tornando

viável a resolução de problemas em períodos mais longos.

Comparando os métodos de resolução que utilizam árvore e latisse, tem-se que

por latisse os subproblemas são maiores que os utilizando árvore, porém são em menor

número, uma vez que há apenas um subproblema por estágio. Caso a árvore de a�uências

fosse dividida da mesma forma que a latisse, o tamanho dos seus subproblemas seriam

maiores, dado que a quantidade de ramos e nós na árvore cresce de forma exponencial e

não linear como na latisse. Além disso, caso ainda se quisesse diminuir o tamanho dos

subproblemas da latisse, poderia se subdividi-los de forma que não se separasse cada um

dos nós �nais de cada estágio em mais de um subproblema, ou seja, poderia se ter um

subproblema para cada nó �nal de cada estágio da latisse. Desse modo, a utilização de

latisse se mostra vantajosa em relação ao uso da árvore de a�uências para a resolução do

problema da otimização do despacho hidrotérmico devido ao menor esforço computacional

exigido.

Quanto aos resultados obtidos pelo estudo de caso, mesmo não havendo a total

coincidência dos resultados encontrados pela árvore e pela latisse para o caso de menor

volume inicial, o valor obtido das funções objetivos foram bem próximos, dando margem de

diferença de menos de 0,06%, uma diferença que para alguns casos pode ser considerada

pequena comparada à vantagem computacional ganha pela utilização da latisse. Além

disso, pelos resultados obtidos pelo primeiro estudo de caso, cujo volume inicial era maior,

acredita-se que houve a coincidência dos valores obtidos do custo esperado de operação

devido há uma maior abundância de água, que permitiu uma maior �exibilidade nos

volumes dos reservatórios.

Por �m, com relação a aplicação da métrica de risco, o CVaR, tem-se que

83

com a aplicação desta ferramenta o custo esperado foi maior do que no caso em que o

CVaR não foi utilizado, uma vez que o CVaR tem o objetivo de garantir que não ocorra

um dé�cit de geração nos meses futuros devido à incerteza das vazões e, portanto, faz

com que haja um maior despacho de energia provinda de termoelétricas, resultando em

reservatórios de água mais cheios, que poderão ser utilizados no caso de ocorrência dos

cenários mais severos. Além disso, foi possível comprovar a possibilidade da aplicação do

CVaR no modelo que utiliza a latisse.

Assim, a proposta apresentada nesta dissertação de substituição da árvore de

a�uências por latisse mostrou-se viável, inclusive junto da aplicação do CVaR, que é atual

métrica de risco utilizada nos modelos NEWAVE e DECOMP, sendo necessário ainda

realizar mais estudos para que possa ocorrer de fato a implementação desta estrutura.

Além disso, como pode ser observado no decorrer desta dissertação, para que

esta proposta seja implementada, não seria necessário alterar os modelos NEWAVE ou

DECOMP, apenas adaptá-los para o recebimento da nova estrutura de dados (latisse)

que seria gerada pelo modelo GEVAZP, através da alteração do seu arquivo com as séries

geradas e o arquivo das probabilidades associadas aos cenários hidrológicos.

7.1 Sugestões para Trabalhos Futuros

O problema foi resolvido de forma simpli�cada, podendo-se ainda acrescentar

mais restrições, ou poderia-se considerar ainda mais usinas tanto termoelétricas quanto

hidroelétricas, aumentando-se assim o número de variáveis do problema e combinações

possíveis de respostas.

Caso se quisesse ir mais a fundo, poder-se-ia tentar implementar de fato esta

alteração no modelo computacional DECOMP.

Além disso, para trabalhos futuros, pode-se sugerir um estudo sobre o tempo

computacional ganho com a substituição da árvore pela latisse.

84

Referências

ALVES, H. Sistemas de Energia Elétrica. [S.l.]: Centro Federal de Educação Tecnológicado Maranhão, 2007.

ANEEL. Atlas da Energia Elétrica do Brasil. 2008. Disponível em: <http://www2.aneel.gov.br/arquivos/PDF/atlas3ed.pdf>.

BERTHO, R. J. Programação Linear com Controle de Risco para o Planejamento daOperação do SIN. Dissertação (Mestrado) � Escola de Engenharia de São Carlos, Uni-versidade de São Paulo, 2013.

CASTRO, N. J. de et al. Plano decenal de expansão de energia - pde 2020: Análise dométodo, metas e riscos. Texto de Discussão do Setor Elétrico, n. 44, p. 26, Janeiro 2012.

CCEE. Representação de Aversão a Risco no Planejamento da Operação Hidrotérmica.2012. Disponível em: <www.ccee.org.br/ccee/documentos/CCEE_042638>.

CCEE. Câmara de Comercialização de Energia Elétrica. 2014. Disponível em: <www.ccee.com.br>.

CEPEL. Manual de Referência do Modelo GERAVAZP. [S.l.], 2010.

CEPEL. Manual de Referência do Modelo DECOMP. 16.6.26.. ed. [S.l.], 2012a.

CEPEL. Manual de Referência do Modelo NEWAVE. [S.l.], 2012b.

CPAMP. Desenvolvimento, implementação e testes de validação da metodologia para in-ternalização do mecanismo de aversão a risco CVaR no programa computacional DE-COMP para estudos energéticos e formação de preço. Brasilia, 2013a.

CPAMP. Mecanismos de Aversão a Risco no Planejamento da Operação deLongo, Médio e Curto Prazos. 2013b. Disponível em: <http://www.mme.gov.br/documents/10584/1139101/1_-_Aspectos_Metodologicos_SAR_CVaR.pdf/19030c06-cb32-4010-90da-81ecb7307f35>.

DEUS, F. L. D. Geração de Cenários de Energia Natural A�uente Através da Adoção daMetodologia do Hipercubo Latino. Dissertação (Mestrado) � UFRJ, Rio de Janeiro, 2010.

FURNAS. Sistema Furnas de Geração e Transmissão. 2016. Disponível em: <http://www.furnas.com.br/hotsites/sistemafurnas/usina_term_funciona.asp>.

http://www2.aneel.gov.br/arquivos/PDF/atlas3ed.pdf

http://www2.aneel.gov.br/arquivos/PDF/atlas3ed.pdf

www.ccee.org.br/ccee/documentos/CCEE_042638

www.ccee.com.br

www.ccee.com.br

http://www.mme.gov.br/documents/10584/1139101/1_-_Aspectos_Metodologicos_SAR_CVaR.pdf/19030c06-cb32-4010-90da-81ecb7307f35



http://www.furnas.com.br/hotsites/sistemafurnas/usina_term_funciona.asp

http://www.furnas.com.br/hotsites/sistemafurnas/usina_term_funciona.asp

85

KLIGERMAN, A. S. Operação Ótima de Subsistemas Hidrotérmicos Interligados Utili-zando Programação Dinâmica Estocástica Dual. Dissertação (Mestrado) � UniversidadeEstadual de Campinas (UNICAMP), Campinas, Fevereiro 1992.

LASDON, L. S. Optimization theory for large systems. New York, NY: [s.n.], 1970.

LOPES, J. E. G. Modelo de Planejamento da Operação de Sistemas Hidrotérmicos deProdução de Energia Elétrica. Tese (Doutorado) � Escola Politécnica da Universidadede São Paulo, 2007.

LUENBERGER, D. G. Investment Science. [S.l.]: Oxford University Press, 1998. 494 p.

MARCATO, A. L. M. et al. Comparação da incorporação da medida de risco cvar entre asmetodologias pde e pdde aplicadas ao planejamento de médio prazo. Simpósio Brasileirode Pesquisa Operacional, p. 924�936, Setembro 2013. Disponível em: <http://www.din.uem.br/sbpo/sbpo2013/pdf/arq0314.pdf>.

ONS. Operador Nacional do Sistema Elétrico. 2014. Disponível em: <www.ons.org.br>.

ONS. Operador Nacional do Sistema Elétrico. 2016. Disponível em: <http://www.ons.org.br/historico/carga_propria_de_energia.aspx>.

PEREIRA, A. O. Modelagem da Operação Ótima da Indústria de Energia Elétrica noBrasil. Dissertação (Mestrado) � Universidade Federal do Rio de Janeiro, 2000.

SHAPIRO, A. et al. Risk neutral and risk averse stochastic dual dynamic programmingmethod. p. 32, 2012.

SOARES, S. Planejamento da operação de sistemas hidrotérmicos. SBA: Controle deAutomação, v. 1, p. 122�131, 1987.

TOLMASQUIM, M. T. Energia Termelétrica: Gás Natural, Biomassa, Carvão, Nuclear.Rio de Janeiro: EPE, 2016.

http://www.din.uem.br/sbpo/sbpo2013/pdf/arq0314.pdf

http://www.din.uem.br/sbpo/sbpo2013/pdf/arq0314.pdf

www.ons.org.br

http://www.ons.org.br/historico/carga_propria_de_energia.aspx

http://www.ons.org.br/historico/carga_propria_de_energia.aspx

86

APÊNDICE

A Cálculo das Probabilidades de Ocorrência de Cada

Cenário

Este apêndice visa apresentar os cálculos das probabilidades de ocorrência de

cada um dos ramos e nós da árvore de a�uências e da latisse.

A.1 Árvore de A�uências

Para uma árvore binária de a�uências, cuja probabilidade de ocorrência de um

cenário otimista é P1 e de ocorrência de um cenário pessimista é P2, a probabilidade ptj de

ocorrência de cada nó é igual a probabilidade ptij de ocorrência do ramo que leva à este nó,

e a probabilidade de ocorrência de cada ramo é igual à probabilidade P1 ou P2 multiplicada

pela ocorrência do seu nó inicial. Assim, estas probabilidades são apresentadas na Figura

A.1 e calculadas de forma genérica pelas Equações A.1 e A.2, respectivamente.


p00 = 1 p11 = 1

p21

p22

p31

p32

p33

p34

pN2(N−1)

pN1

pN2

pN2(N−1)−1

p101 = 1

p211 = P1

p212 = P2

p311

p312

p323

p324

p(N−1)1

pN11

pN12

p(N−1)2(N−2)

pN(2N−2)(2(N−1)−1)

pN(2N−2)(2(N−1))

Figura A.1: Representação das probabilidades de uma árvore binária multi-estágios.

87

ptj = ptij (A.1)

ptij = Psp(t−1)i (A.2)

Sendo:

s =

{1, se j for ímpar

2, se j for par(A.3)

A.2 Latisse

Agora, para uma latisse binomial, cuja probabilidade de ocorrência de um

cenário otimista é P1 e de ocorrência de um cenário pessimista é P2, a probabilidade ptjde ocorrência de cada nó é igual a soma das probabilidades ptij de ocorrência dos ramos

que levam à este nó, e a probabilidade de ocorrência de cada ramo é igual à probabilidade

P1 ou P2 multiplicada pela ocorrência do seu nó inicial. Assim, estas probabilidades são

apresentadas na Figura A.2 e calculadas de forma genérica pelas Equações A.4 e A.5,

respectivamente.


p00 p11

p21

p22

p31

p32

p33

pNN

pN1

pN2

p101 = 1

p211 = P1

p212 = P2

p311

p312

p323

p324

pN11

pN12p(N−1)1

p(N−1)(N−1))

pN(N−1)N

Figura A.2: Representação das probabilidades de uma latisse binomial multi-estágios.

88

ptj =∑i

ptij (A.4)

ptij = PsP(t−1)i (A.5)

Sendo:

s =

{1, se (i+ j) for par

2, se (i+ j) for ímpar(A.6)

B Resultados Obtidos para o Estudo de Caso

Nas Tabelas B.1 à B.4, são apresentados os valores obtidos pela resolução do

estudo de caso do capítulo 6, com o volume inicial de 6.000MWmed, por meio da utilização

de árvore ou latisse, tanto para o caso neutro ao risco quanto averso ao risco.

Tabela B.1: Resultados obtidos pela otimização utilizando árvore para o caso neutro aorisco, com volume inicial equivalente à 10.000MWmed.

Estágio (t) Ramo(ij) GH tij X t

ij1 X tij2 X t

ij3 X tij4 Custotij V t

j

1 01 6838,6 499,2 350,0 249,2 0 185609,2 13837,5

211 7424,2 499,2 0 0 0 63418,4 13011,312 6744,2 499,2 350,0 330,0 0 202690,4 10627,8

3

11 7751,6 195,2 0 0 0 24798,2 9260,212 7097,6 499,2 350,0 0 0 132928,4 8848,423 6779,9 499,2 350,0 317,7 0 200090,1 8848,424 6414,1 499,2 350,0 330,0 353,5 308740,4 8148,4

4

11 8145,9 0 0 0 0 0 400012 7378,5 499,2 268,2 0 0 116682,9 400023 7734,1 411,8 0 0 0 52315,1 400024 6966,7 499,2 350,0 330,0 0 202690,4 400035 7734,1 411,8 0 0 0 52315,1 400036 6966,7 499,2 350,0 330,0 0 202690,4 400047 7034,1 499,2 350,0 262,6 0 188442,0 400048 6266,7 499,2 350,0 330,0 700,0 412690,4 4000

Tabela B.2: Resultados obtidos pela otimização utilizando latisse para o caso neutro aorisco, com volume inicial equivalente à 10.000MWmed.


ij1 X tij2 X t


j

1 01 6838,6 499,2 350,0 249,2 0 185609,2 13837,5

211 7424,2 499,2 0 0 0 63418,4 13011,312 6744,2 499,2 350,0 330,0 0 202690,4 10627,8

3

11 7751,6 195,2 0 0 0 24798,2 9260,212 7097,6 499,2 350,0 0 0 132928,4

8848,422 6779,9 499,2 350,0 317,7 0 200090,123 6414,1 499,2 350,0 330,0 353,5 308740,4 8148,4

4

11 8145,9 0 0 0 0 0 400012 7378,5 499,2 268,2 0 0 116682,9

400022 7734,1 411,8 0 0 0 52315,123 6966,7 499,2 350,0 330,0 0 202690,4

400033 7034,1 499,2 350,0 262,6 0 188442,034 6266,7 499,2 350,0 330,0 700,0 412690,4 4000

90

Tabela B.3: Resultados obtidos pela otimização utilizando árvore para o caso averso aorisco, com volume inicial equivalente à 10.000MWmed.


ij1 X tij2 X t


j

1 01 6757,8 499,2 350,0 330,0 0 202690,4 13918,3

211 7424,2 499,2 0 0 0 63418,4 13092,112 6744,2 499,2 350,0 330,0 0 202690,4 11708,6

3

11 7832,4 114,4 0 0 0 14533,4 9260,212 7178,4 499,2 269,2 0 0 116881,5 8848,423 6860,7 499,2 350,0 236,9 0 183009,0 8848,424 6464,9 499,2 350,0 330,0 272,7 284500,4 8148,4

4

11 8145,9 0 0 0 0 0 400012 7378,5 499,2 268,2 0 0 116682,9 400023 7734,1 411,8 0 0 0 52315,1 400024 6966,7 499,2 350,0 330,0 0 202690,4 400035 7734,1 411,8 0 0 0 52315,1 400036 6966,7 499,2 350,0 330,0 0 202690,4 400047 7034,1 499,2 350,0 262,6 0 188442,0 400048 6266,7 499,2 350,0 330,0 700,0 412690,4 4000

Tabela B.4: Resultados obtidos pela otimização utilizando latisse para o caso averso aorisco, com volume inicial equivalente à 10.000MWmed.


ij1 X tij2 X t


j

1 01 6757,8 499,2 350,0 330,0 0 202690,4 13918,3

211 7424,2 499,2 0 0 0 63418,4 13092,112 6744,2 499,2 350,0 330,0 0 202690,4 11708,6

3

11 7832,4 114,4 0 0 0 14533,4 9260,212 7178,4 499,2 269,2 0 0 116881,5

8848,422 6860,7 499,2 350,0 236,9 0 183009,023 6464,9 499,2 350,0 330,0 272,7 284500,4 8148,4

4

11 8145,9 0 0 0 0 0 400012 7378,5 499,2 268,2 0 0 116682,9

400022 7734,1 411,8 0 0 0 52315,123 6966,7 499,2 350,0 330,0 0 202690,4

400033 7034,1 499,2 350,0 262,6 0 188442,034 6266,7 499,2 350,0 330,0 700,0 412690,4 4000

Os valores de geração e volume estão apresentados em MWmed e os valores de

custo em reais.

Por �m, nas Tabelas B.5 e B.6, são apresentados os valores obtidos para cada

variável durante a resolução do problema do estudo de caso que apresenta o volume inicial

91

do reservatório equivalente à 9.000MWmed, por meio da utilização de árvore e latisse,

respectivamente, para o caso neutro ao risco.

Tabela B.5: Resultados obtidos pela otimização utilizando árvore para o caso neutro aorisco, com volume inicial equivalente à 9.000MWmed.


ij1 X tij2 X t


j

1 01 6757,8 499,2 350,0 330,0 0 202690,4 12918,3

211 7074,2 499,2 350,0 0 0 132928,4 12442,112 6171,5 499,2 350,0 330,0 572,7 374500,4 11281,3

3

11 7447,6 499,2 0 0 0 63418,37 8995,012 6767,6 499,2 350,0 330,0 0 202690,4 8609,223 6767,6 499,2 350,0 330,0 0 202690,4 8514,224 6067,6 499,2 350,0 330,0 700,0 412690,4 8148,4

4

11 7880,7 265,2 0 0 0 33691,01 400012 7113,3 499,2 350,0 183,4 0 171699,1 400023 7494,9 499,2 151,8 0 0 93565,85 400024 6727,5 499,2 350,0 330,0 239,2 274450,4 400035 7399,9 499,2 246,8 0 0 112432,8 400036 6632,5 499,2 350,0 330,0 334,2 302950,4 400047 7034,1 499,2 350,0 262,6 0 188442,0 400048 6266,7 499,2 350,0 330,0 700,0 412690,4 4000

Tabela B.6: Resultados obtidos pela otimização utilizando latisse para o caso neutro aorisco, com volume inicial equivalente à 9.000MWmed.


ij1 X tij2 X t


j

1 01 6757,8 499,2 350,0 330,0 0 202690,4 12918,3

211 7169,2 499,2 255,0 0 0 114061,4 12347,112 6171,5 499,2 350,0 330,0 572,7 374500,4 11281,3

3

11 7447,6 499,2 0 0 0 63418,37 8900,012 6767,6 499,2 350,0 330,0 0 202690,4

8514,222 6767,6 499,2 350,0 330,0 0 202690,423 6067,6 499,2 350,0 330,0 700,0 412690,4 8148,4

4

11 7785,7 360,2 0 0 0 45759,81 400012 7018,3 499,2 350,0 278,4 0 191782,1

400022 7399,9 499,2 246,8 0 0 112432,823 6632,5 499,2 246,8 0 0 302950,4

400033 7034,1 499,2 350,0 262,6 0 188442,034 6266,7 499,2 350,0 330,0 700,0 412690,4 4000

Documents

Representação da Incerteza em Modelos de …...UNIVERSIDADE ESTADUAL DE CAMPINAS aculFdade de Engenharia Mecânica ernandaF Nakano Kazama Representação da Incerteza em Modelos