OTIMIZAÇÃO COM EVOLUÇÃO DIFERENCIADA ADAPTATIVA COM DUPLA
SUPERFÍCIE DE RESPOSTA DO VPL DE UMA UTE COMPLEXA E FLEXÍVEL
Bruno Itagyba Paravidino
Dissertação
de
Mestrado
apresentada
ao
Programa de Pós-graduação em Engenharia
Mecânica, COPPE, da Universidade Federal do
Rio de Janeiro, como parte dos requisitos
necessários à obtenção do título de Mestre em
Engenharia Mecânica.
Orientador: Marcelo José Colaço
Rio de Janeiro
Setembro de 2013
OTIMIZAÇÃO COM EVOLUÇÃO DIFERENCIADA ADAPTATIVA COM DUPLA
SUPERFÍCIE DE RESPOSTA DO VPL DE UMA UTE COMPLEXA E FLEXÍVEL
Bruno Itagyba Paravidino
DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DO INSTITUTO ALBERTO
LUIZ COIMBRA DE PÓS-GRADUAÇÃO E PESQUISA DE ENGENHARIA
(COPPE) DA UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE
DOS REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE
EM CIÊNCIAS EM ENGENHARIA MECÂNICA.
Examinada por:
________________________________________________
Prof. Marcelo José Colaço, Ph.D.
________________________________________________
Prof. Manuel Ernani de Carvalho Cruz, Ph.D.
________________________________________________
Dr. Leonardo dos Santos Reis Vieira, D.Sc.
RIO DE JANEIRO, RJ - BRASIL
SETEMBRO DE 2013
iii
Paravidino, Bruno Itagyba
Otimização com Evolução Diferenciada Adaptativa
com Dupla Superfície de Resposta do VPL de uma UTE
Complexa e Flexível/ Bruno Itagyba Paravidino. – Rio de
Janeiro: UFRJ/COPPE, 2013.
XII, 161 p.: il.; 29,7 cm.
Orientador: Marcelo José Colaço
Dissertação (mestrado) – UFRJ/ COPPE/ Programa de
Engenharia Mecânica, 2013.
Referências Bibliográficas: p. 155-161.
1. Evolução Diferenciada Adaptativa. 2. Superfície de
Resposta. 3. Otimização. I. Colaço, Marcelo José. II.
Universidade Federal do Rio de Janeiro, COPPE,
Programa de Engenharia Mecânica. III. Título.
iv
As minhas meninas
Lu e Lara.
v
AGRADECIMENTOS
Agradeço aos meus pais por sempre me incentivarem aos estudos. A minha
esposa Lu e a minha linda filha Lara.
Agradeço aos meus colegas de trabalho por me ajudarem a conclusão do
trabalho, principalmente ao Thiago, Franciele e Alessandro.
vi
Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos
necessários para a obtenção do grau de Mestre em Ciências (M.Sc.)
OTIMIZAÇÃO COM EVOLUÇÃO DIFERENCIADA ADAPTATIVA COM DUPLA
SUPERFÍCIE DE RESPOSTA DO VPL DE UMA UTE COMPLEXA E FLEXÍVEL
Bruno Itagyba Paravidino
Setembro/2013
Orientador: Marcelo José Colaço
Programa: Engenharia Mecânica
Este trabalho apresenta um algoritmo de otimização, o ED2SR, que emprega
uma função interpolada no lugar de uma função objetivo a ser maximizada na
otimização de uma Usina Termelétrica complexa, a partir do método heurístico
Evolução Diferenciada. A função objetivo é representada pelo Valor Presente Líquido
de um empreendimento térmico, que leva em consideração um risco associado ao
investimento de apenas 5%, bem como as diversas flexibilidades operacionais existentes
num projeto de investimento em usinas termelétricas no Brasil. Para auxiliar as soluções
dos balanços de massa e energia de cada projeto e evitar a obtenção dos modelos físico
e termodinâmico de forma explícita, este trabalho faz uso do simulador de processos
Thermoflow. Durante o processo de maximização, para o cálculo da função objetivo, o
simulador de processo é chamado diversas vezes até que se obtenha um valor ótimo. É
constatado neste trabalho que a aplicação do modelo de Superfície de Resposta
formulado em funções de base radial, como método de interpolação da função objetivo
original, implica a redução do número de avaliações da função original, sem prejuízo
significante quando da obtenção do valor ótimo. Este trabalho mostra o
desenvolvimento de uma metodologia de investimento, que ao otimizar os parâmetros
do processo de geração de uma térmica, garante a lucratividade do empreendedor e
elimina o risco do seu Valor Presente Líquido ser negativo.
vii
Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the
requirements for the degree of Master of Science (M.Sc.)
ADAPTIVE DIFFERENTIAL EVOLUTION WITH DOUBLE RESPONSE
SURFACE OPTIMIZATION OF A COMPLEX AND FLEXIBLE POWER PLANT’S
NPV
Bruno Itagyba Paravidino
September/2013
Advisor: Marcelo José Colaço
Department: Mechanical Engineering
This work presents an optimization algorithm, the ED2SR, which makes use of
an interpolation function in order to replace the objective function to be maximized in
the optimization of a complex power plant system. A heuristic method, called
Differential Evolution, was used for this task. The objective function is represented by
the Net Present Value of the project, which takes into account an investment risk of 5%,
as well as other investment flexibilities in power plants usually found in Brazil. To
solve the mass and energy conservation equations for each project and to avoid explicit
physical and thermodynamic models, this work makes use of the Thermoflow processes
simulator software. During the maximization process of the objective function, the
process simulator is called repeatedly until it obtains an optimal value of it. It is found
in this study that the application of a response surface model formulated in terms of
radial basis functions, as an approximation mathematical function of the original
physical problem, implies in a reduction of the number of evaluations of the original
objective function, without significant impairment of the optimal value obtained at the
end of the optimization task. This work shows the development of an investment
methodologies that, by optimizing the power plant’s parameters, eliminates the risk of
Net Present Value be negative and ensure the profitability of the investor.
viii
Sumário
1 INTRODUÇÃO ........................................................................................................ 1
1.1
MOTIVAÇÃO................................................................................................... 1
1.2
OBJETIVO ........................................................................................................ 3
1.3
FERRAMENTAS .............................................................................................. 4
1.3.1
Simulador Termodinâmico Thermoflow .................................................... 4
1.3.2
Plataforma de Cálculo MATLAB............................................................... 7
1.4
2
MÉTODO DE OTIMIZAÇÃO: UMA INTRODUÇÃO .................................. 8
REVISÃO BIBLIOGRÁFICA ............................................................................... 11
2.1
SETOR ELÉTRICO BRASILEIRO ............................................................... 11
2.2
ALGORITMOS EVOLUTIVOS..................................................................... 17
2.2.1
2.3
3
MÉTODO DA SUPERFÍCIE DE RESPOSTA .............................................. 24
COMERCIALIZAÇÃO DE ENERGIA ELÉTRICA NO SEB.............................. 27
3.1
4
Evolução Diferenciada ............................................................................. 19
COMERCIALIZAÇÃO NO ACR .................................................................. 27
3.1.1
Contratos de Disponibilidade ................................................................... 29
3.1.2
Leilões De Energia Nova.......................................................................... 32
3.1.3
O Índice de Custo Benefício ..................................................................... 33
3.1.4
Metodologia de Cálculo do ICB ............................................................... 36
3.1.4.1
Disponibilidade ................................................................................. 36
3.1.4.2
Inflexibilidade ................................................................................... 36
3.1.4.3
Custo Variável Unitário - CVU......................................................... 37
3.1.4.4
Custos Marginais de Operação - CMO ............................................. 39
3.1.4.5
COP e CEC........................................................................................ 41
3.1.4.6
Receita Fixa ....................................................................................... 42
3.2
COMERCIALIZAÇÃO NO ACL................................................................... 43
3.3
MERCADO DE CURTO PRAZO .................................................................. 45
USINA TERMELÉTRICA ..................................................................................... 46
4.1
PLANTA DE GERAÇÃO: CICLO COMBINADO ....................................... 46
4.2
UTE-BASE ...................................................................................................... 48
4.2.1
4.3
UTE-Base no Thermoflow ....................................................................... 49
INTERFACE THERMOFLOW E MATLAB ................................................ 53
4.3.1
Convertendo a UTE-Base Em U-Link ..................................................... 54
4.3.2
Executando a UTE-Base no Matlab ......................................................... 55
ix
5
FUNÇÃO OBJETIVO ............................................................................................ 56
5.1
6
VPL DE UMA UTE FLEXÍVEL .................................................................... 56
5.1.1
Investimento ............................................................................................. 60
5.1.2
Receita do CCEAR ................................................................................... 61
5.1.3
Receita do Contrato PPA .......................................................................... 63
5.1.4
Outros Custos ........................................................................................... 65
5.1.5
VPL da UTE ............................................................................................. 66
5.2
ANÁLISE DE RISCO ..................................................................................... 68
5.3
SIMULAÇÃO DE MONTE CARLO ............................................................. 70
5.4
EXPRESSÃO FINAL PARA A FUNÇÃO OBJETIVO ................................ 72
ALGORITMO DE OTIMIZAÇÃO ........................................................................ 76
6.1
MÉTODO DA EVOLUÇÃO DIFERENCIADA............................................ 76
6.1.1
Método Tradicional .................................................................................. 76
6.1.1.1
Mutação ............................................................................................. 77
6.1.1.2
Cruzamento ....................................................................................... 78
6.1.1.3
Seleção .............................................................................................. 80
6.1.2
Método Adotado: ED_Mod ...................................................................... 81
6.1.2.1
Adaptação dos Fatores de Cruzamento e Mutação ........................... 82
6.1.2.2
Pool de Estratégias de Mutação ........................................................ 87
6.1.2.3
Adaptação das Estratégias ................................................................. 90
6.1.2.4
ED_mod ............................................................................................. 93
6.2
MÉTODO DE SUPERFÍCIE DE RESPOSTA ............................................... 95
6.3
ALGORITMO PROPOSTO: ED2SR ............................................................ 100
6.3.1
H3_mod .................................................................................................. 100
6.3.1.1
6.3.2
6.4
7
SR .................................................................................................... 102
ED2SR .................................................................................................... 104
TESTE DE OTIMIZAÇÃO COM ED2SR: .................................................. 108
UTE ÓTIMA: UM ESTUDO DE CASO ............................................................. 117
7.1
MODELAGEM DO VALOR PRESENTE LÍQUIDO ................................. 117
7.1.1
Avaliação das Variáveis Comuns ........................................................... 118
7.1.1.1
Vida Útil .......................................................................................... 118
7.1.1.2
Taxa Mínima de Atratividade Anual e Mensal ............................... 119
7.1.1.3
Garantia Física da UTE ................................................................... 119
7.1.1.4
Disponibilidade da UTE .................................................................. 120
x
7.1.1.5
Inflexibilidade da UTE .................................................................... 121
7.1.1.6
Custo Variável Unitário – CVU ...................................................... 121
7.1.1.6.1 Custo Variável Unitário de Caráter Estrutural ............................. 123
7.1.1.6.2 Custo Variável Unitário de Caráter Operacional ......................... 123
7.1.1.7
Custo Marginal de Operação – CMO.............................................. 124
7.1.1.8
Preço de Liquidação das Diferenças – PLD .................................... 124
7.1.2
Receita da Comercialização de Energia no Ambiente Regulado ........... 126
7.1.2.1
Valor do ICB ................................................................................... 126
7.1.2.2
Valor do COP e CEC ...................................................................... 127
7.1.2.2.1 Custo de Operação ....................................................................... 127
7.1.2.2.2 Custo Econômico de Curto Prazo ................................................ 128
7.1.3
7.1.3.1
Preço da Energia Contratada ........................................................... 129
7.1.3.2
Custo Variável Unitário no mês t .................................................... 129
7.1.3.3
Número de Horas no Mês................................................................ 130
7.1.4
8
Receita da Comercialização de Energia no Ambiente Livre .................. 128
Custos Diversos ...................................................................................... 130
7.1.4.1
Investimento .................................................................................... 130
7.1.4.2
Outros Custos .................................................................................. 131
7.1.4.3
Impostos Diretos e de Renda........................................................... 132
7.1.4.4
Depreciação ..................................................................................... 132
7.2
VALOR PRESENTE LÍQUIDO DA UTE-BASE ........................................ 133
7.3
OTIMIZANDO A UTE-BASE ..................................................................... 137
7.3.1
Otimizando Com o Algoritmo da ED_mod ............................................ 139
7.3.2
Otimizando Com o Algoritmo da ED2SR .............................................. 147
CONCLUSÃO ...................................................................................................... 152
REFERÊNCIAS BIBLIOGRÁFICAS ......................................................................... 155
xi
NOMENCLATURAS
ACL
Ambiente de Contratação Livre.
ACO
Ant Colony Optimization (Método da Colônia de Formigas).
ACR
Ambiente de Contratação Regulada.
AdapSS
Adaptive Strategy Selection.
AEO
Annual Energy Outlook.
AG
Algoritmos Genéticos (Genetic Algorithms).
ANEEL
Agência Nacional de Energia Elétrica.
BP
Busca em Padrões.
CCEAR
Contratos de Comercialização de Energia Elétrica no Ambiente
Regulado.
CCEE
Câmara de Comercialização de Energia Elétrica.
CEC
Custo Econômico de Curto Prazo.
CMO
Custo Marginal de Operação.
COP
Custo de Operação.
CVU
Custo Variável Unitário.
ED
Evolução Diferenciada (Differential Evolution).
EIA
Energy Information Administration.
EP
Enxame de Partículas.
EPE
Empresa de Pesquisa Energética.
GCPS
Grupo Coordenador do Planejamento dos Sistemas Elétricos.
GEATbx
Genetic and Evolutionary Algorithm Toolbox.
GF
Garantia Física.
ICB
Índice de Custo Benefício.
LEN
Leilão de Energia Nova.
MAE
Mercado Atacadista de Energia Elétrica.
MCP
Mercado de Curto Prazo.
MCSD
Mecanismo de Compensação de Sobras e Déficits para as Distribuidoras.
MME
Ministério de Minas e Energia.
ONS
Operador Nacional do Sistema.
PE
Programação Evolutiva (Evolutionary Programming).
PEACE
Plant Engineering And Construction Estimator.
PLD
Preço de Liquidação das Diferenças.
xii
PNE
Plano Nacional de Energia.
PPA
Power Purchase Agreement.
RE-SEB
Reestruturação do Setor Elétrico Brasileiro.
PSO
Particle Swarm Otimization (Método do Enxame de Partículas).
RS
Recozimento Simulado.
SCL
Sistema de Contabilização e Liquidação.
SEB
Setor Elétrico Brasileiro.
SMS
Simulação de Monte Carlo
TMA
Taxa Mínima de Atratividade
UTE
Usina Termelétrica.
VPL
Valor Presente Líquido.
1
1
1.1
INTRODUÇÃO
MOTIVAÇÃO
Diferente da grande maioria dos países, o Brasil tem sua matriz energética com
participação majoritária de fontes renováveis. O setor elétrico utiliza como base,
principalmente, a energia proveniente de usinas hidrelétricas e as demais fontes de
energia servem como complementares no caso de hidrologia desfavorável ou elevação
de demanda. As usinas termelétricas se enquadram neste segundo caso, sendo chamadas
a despachar em períodos secos, em que os reservatórios das usinas hidrelétricas estão
com seus níveis baixos, e também para atendimento da demanda nos horários de ponta.
Desde o início da década de 90 o Setor Elétrico Brasileiro (SEB) está sendo
profundamente reestruturado, com o objetivo de introduzir a livre competição nos
segmentos de geração e de comercialização, assim como possibilitar a inserção de
novos agentes na prestação dos serviços de energia elétrica (conduzida pelo Ministério
de Minas e Energia), preservando a lucratividade e o interesse destes em investir em
geração.
O modelo vigente de comercialização de energia elétrica no Brasil estabelece a
contratação da energia para atendimento da expansão do mercado consumidor das
distribuidoras por meio dos Leilões de Energia Nova, sendo o menor preço ofertado o
critério para a seleção dos empreendimentos vencedores.
Entre as vantagens da geração termelétrica, quando comparada com uma
hidrelétrica, estão o menor tempo para instalação, o menor investimento unitário
instalado (quantificado em R$/KW) e a flexibilidade para o atendimento da demanda.
Porém, investimentos em usinas termelétricas são classificados como de alto risco
devido ao alto custo da geração térmica, incerteza da taxa de câmbio, incerteza do preço
2
do gás natural, restrição de oferta de gás natural e a baixa probabilidade de ocorrer uma
hidrologia desfavorável, dentre outros.
A flexibilidade operacional imposta pela nova estrutura do setor, onde a
termelétrica flexível somente será despachada se o preço da energia no Mercado de
Curto Prazo (preço spot) estiver acima do seu custo operacional, agrega valor ao
investimento. Com isso, é imprescindível a obtenção de um modelo de avaliação de
projetos que capture as reais vantagens provenientes da flexibilidade operacional de
uma planta de geração termelétrica [1], que transforme esse tipo de operação flexível
em valor (receita), para que seja adicionado ao Valor Presente Líquido (VPL) do
empreendimento.
As características supracitadas tornam o Sistema Elétrico Brasileiro bastante
peculiar em sua configuração. Nesse sentido, o projeto de usinas termelétricas deve
atender as exigências de um mercado baseado em contrato bilateral na modalidade
Disponibilidade de Energia Elétrica (de acordo com o Decreto n° 5.163 de 30/7/2004 e
Lei n° 10.848 de 15/3/2004, art. 27 e 28, [2]) para geração de energia e ainda manter a
atratividade do negócio de geração para o investidor, dando a possibilidade de atuar no
mercado spot. Dado este contexto, a grande vantagem para a produção termelétrica está
no fato da térmica se declarar flexível perante o Operador Nacional do Sistema, ONS,
ao invés de atuar de forma inflexível1 [1].
Diante disso, é necessário desenvolver metodologias que garantam projetos
termelétricos otimizados e com baixos custos, levando em consideração as reais
vantagens provenientes da flexibilidade operacional de uma planta de geração
termelétrica instalada no Brasil.
1
Inflexibilidade corresponde ao montante de potência mínima que deve ser obrigatoriamente despendido
e que não está sujeito, portanto, à regra de despacho do ONS.
3
1.2
OBJETIVO
Este trabalho teve como objetivo otimizar um ciclo combinado com três níveis
de pressão de vapor, reaquecimento, extração e indução, partir da avaliação das
possibilidades de ganho (receitas) provenientes de variáveis incertas como: atuar no
mercado spot (Mercado de Curto Prazo) ou no mercado de contratos bilaterais. Ou seja,
a otimização do ciclo de geração teve como função objetivo o Valor Presente Líquido
de uma térmica segundo as regras impostas pelo Setor Elétrico Brasileiro, levando em
consideração um risco associado ao investimento de apenas 5% e as diversas
flexibilidades existentes num projeto de investimento em usinas termelétricas no Brasil
[1-3].
Para o cálculo da função objetivo não foram desenvolvidos modelos físicos e
termodinâmicos de forma explícita para o ciclo de geração. Neste trabalho o simulador
de processo termodinâmico Thermoflow1 foi utilizado com o propósito de resolver os
balanços de massa e de energia para que o VPL fosse então calculado.
O principal objetivo da escolha do Thermoflow foi porque este apresenta um
banco de dados de custos completo para empreendimentos de usinas termelétricas. Nele,
são contemplados individualmente os principais sistemas e equipamentos, materiais,
custos de engenharia, custos de construção e montagem, todos compatíveis com os
preços praticados nos EUA. Um banco de dados dessa magnitude foi uma ferramenta de
grande valia para os cálculos de investimento do empreendimento, parte integrante da
função objetivo.
Para otimizar a planta termelétrica foi desenvolvida uma ferramenta
computacional que integra o Thermoflow com uma plataforma de cálculo
computacional, o MATLAB.
1
Pacote de programas de engenharia desenvolvido para a indústria de geração elétrica e cogeração.
4
Foi desenvolvido um código computacional que integra o algoritmo de
otimização baseado no método heurístico Evolução Diferenciada (ED) [4] com a função
objetivo, conforme visto na referência [5]. Como este trabalho utilizou um simulador de
processo termodinâmico para o cálculo da função objetivo, foi aplicado o método de
Superfície de Resposta na otimização da função objetivo com o intuito em diminuir o
número de simulações executadas pelo Thermoflow, ou seja, diminuir o número de
avaliações da função objetivo, reduzindo o custo computacional.
A grande diferença em relação ao que a referência [5] já aplicou em termos de
otimização de sistemas térmicos, foi a proposta de uma função objetivo, que leva em
consideração os riscos de um projeto de investimento termelétrico no Brasil; a
abordagem de modelos recentes do método da Evolução Diferenciada, considerando a
adaptação de seus parâmetros; a consideração de outros tipos de funções de base radial
no método de interpolação; a consideração de um critério baseado no erro dos mínimos
quadrados para a escolha do parâmetro de forma presente no método de Superfície de
Resposta, com intuito de melhorar a aproximação e contornar o mal condicionamento
da matriz de interpolação, conforme visto na referência [6]; a proposta de um algoritmo
de otimização que integra o método da Evolução Diferenciada Adaptativa com o
modelo de Superfície de Resposta, o ED2SR; e um teste de otimização utilizando o
ED2SR na otimização de cinco funções conhecidas: Beale; Branin; Easom; Sphere; e
Griewank.
1.3
1.3.1
FERRAMENTAS
Simulador Termodinâmico Thermoflow
O Thermoflow é um pacote de programas (software) de engenharia térmica
5
desenvolvido para a indústria de energia e cogeração, sendo comercializado desde 1987.
Com os programas do Thermoflow é possível modelar diferentes plantas de geração
termelétrica. Dentre os programas do pacote destacam-se: GTpro; GTMaster; SteamPro;
SteamMaster; e Thermoflex.
• GTpro: Este foi o primeiro programa desenvolvido para o pacote, sendo um
programa de projeto de plantas de geração e cogeração utilizando turbogeradores a gás.
Baseado em critérios de projeto termodinâmicos definidos pelo usuário, o programa
gera balanços de massa e energia e o projeto conceitual dos equipamentos mais
relevantes na planta. Aliado ao módulo PEACE (Plant Engineering And Construction
Estimator) gera ainda como dado de saída uma estimativa do custo total do
empreendimento, podendo, portanto, ser utilizado para estudos de viabilidade técnicoeconômicos;
• GTMaster: Este programa também simula modelos de plantas de cogeração e
ciclo combinado. No entanto, diferentemente do GTPro, seu objetivo principal é estudar
o comportamento de uma determinada planta quando operando em condições diferentes
das condições de projeto (off-design). Dessa forma, o GTMaster permite então estudar
diferentes cenários operacionais tanto para novos projetos quanto em plantas já
existentes;
• SteamPro: Neste programa é possível criar plantas de geração termelétrica que
sejam baseadas no ciclo Rankine. Assim como o GTPro, o usuário é responsável por
definir os principais critérios de projeto, tendo o auxílio do banco de dados do próprio
programa para decisão de parâmetros ainda não conhecidos na fase de projeto. Também
é utilizado junto com o módulo PEACE para estudos de viabilidade técnico-econômica
deste tipo de unidade;
6
• SteamMaster: Este o é equivalente do programa GTMaster para o SteamPro.
Permite estudar o comportamento de unidades novas já definidas ou unidades existentes
em condições off-design;
• Thermoflex: Diferente dos programas descritos anteriormente, que foram
desenvolvidos para atender projetos de plantas específicas (ciclos Rankine nos casos do
SteamPro e SteamMaster e ciclos Brayton+Rankine nos casos do GTPro e GTMaster), o
Thermoflex segue a filosofia fully-flexible. Com ele o usuário pode construir o modelo
de sua unidade com todas as particularidades necessárias, que de outra forma não
poderiam ser modeladas nos outros programas do pacote Thermoflow. No entanto, sua
maior flexibilidade é conseguida em detrimento da facilidade da construção dos
modelos que os outros programas possuem, justamente por serem destinados a
aplicações específicas. Portanto, o Thermoflex exige do usuário um maior
conhecimento técnico de ciclos e equipamentos e maior cuidado na elaboração dos
modelos, e certamente um maior dispêndio de tempo. Outra grande diferença entre os
outros programas e o Thermoflex é que este último não dispõe do módulo PEACE.
Além destes principais programas, o Thermoflow possui utilitários compatíveis
com a plataforma Excel, que permitem ao usuário explorar ao máximo o banco de dados
de plantas de geração termelétrica que o programa oferece, tanto no estudo de
elaboração de novos projetos quanto no estudo de alterações de plantas existentes.
Alguns dos utilitários mais conhecidos, compatíveis com a plataforma Excel, são o ELink e o TOPS. Entretanto, estes não serão abordados neste trabalho.
Outro utilitário menos conhecido é o U-Link. U-Link é a interface ActiveX do
software Thermoflow. Ele funciona em ambientes de programação que suportam
ActiveX dynamics link libraries (dll), como por exemplo: MATLAB; MS Excel; MS
7
Visual Base C; dentre outros. Ou seja, o U-Link permite o usuário criar uma interface
com modelos de plantas previamente construídas a partir dos programas do Thermoflow
(GTpro, GTMaster, SteamPro, SteamMaster, e Thermoflex) e manipulá-las em outro
ambiente de programação.
Com isso, o usuário tem a liberdade de utilizar os resultados previstos pelo
modelo do Thermoflow com outras informações, como por exemplo: informações do
desempenho previsto; de banco de dados financeiros; informações da disponibilidade do
equipamento mecânico; dados meteorológicos; previsão de demanda; etc.
Neste trabalho foram utilizados os seguintes programas: GTpro e U-Link.
1.3.2
Plataforma de Cálculo MATLAB
O MATLAB é uma plataforma da empresa MathWorks com linguagem de
programação de alto nível. Esta foi a plataforma utilizada para a criação das rotinas de
otimização através da implementação de algoritmos e com auxílio de procedimentos
prontos pertencentes à biblioteca do programa.
Através do MATLAB foram desenvolvidas as rotinas para o cálculo da função
objetivo, do algoritmo de otimização e para construção da superfície de resposta. É
durante a execução da função objetivo que o simulador de processo Thermoflow é
requisitado. Um conjunto de parâmetros de operação (variáveis de decisão) gerado pelo
algoritmo de otimização é fornecido ao Thermoflow através de linhas de comando
escritas no MATLAB, o qual também ordena que a simulação seja executada e por fim
extrai todas as informações necessárias para o cálculo do VPL de um empreendimento
termelétrico.
8
1.4
MÉTODO DE OTIMIZAÇÃO: UMA INTRODUÇÃO
O algoritmo de otimização desenvolvido neste trabalho foi baseado no Método
da Evolução Diferenciada (ED). Este método foi proposto por [4] em 1995 como uma
nova abordagem heurística baseada segundo a teoria de Darwin, onde os membros mais
fortes de uma população são mais capazes de sobreviver em uma determinada condição
ambiental.
O algoritmo da ED tradicional consiste em gerar aleatoriamente (a partir de uma
distribuição de probabilidade uniforme) uma população com NP indivíduos dentro do
domínio delimitado, sendo cada um dos NP indivíduos representado por um vetor de
variáveis de otimização. Para cada vetor (indivíduo) é calculada a função objetivo,
sendo este procedimento chamado de avaliação (evaluation).
Assim, esta população inicial sofre modificações a partir dos processos de
mutação (mutation), cruzamento (crossover) e seleção (selection), construindo-se uma
nova população (ou nova geração) a cada iteração.
Na realidade, a intenção dos processos de mutação e cruzamento é gerar
perturbações nos indivíduos da população corrente, para que assim, seja feita uma busca
pelo espaço de soluções.
O primeiro processo consiste em gerar um vetor mutante (mutant vector) para
cada vetor da população corrente. Este vetor é formado a partir de uma combinação
linear entre dois ou mais vetores escolhidos aleatoriamente da população corrente. Esta
combinação é ponderada por um fator de mutação F, sendo este fator mantido constante
durante todo o processo de otimização.
O segundo processo, o cruzamento, consiste em formar outro vetor, chamado de
vetor experimental (trial vector). Este vetor é composto pelas componentes (variáveis
de otimização) do vetor da população corrente e pelas componentes do vetor mutante
9
(criando no processo de mutação), ou seja, o vetor experimental é a mistura das
componentes dos outros dois vetores citados. A escolha de qual componente irá compor
o vetor experimental é feita a partir da comparação de um número aleatório (entre 0 e 1)
com o fator de cruzamento CR, sendo este fator, assim com o fator de mutação,
mantido constante durante todo o processo iterativo.
Por fim, a função objetivo é avaliada para o vetor experimental. No caso de uma
minimização, se o valor da função objetivo aplicada ao vetor experimental for menor do
que o valor dela aplicada ao vetor da população atual, então o vetor experimental
substitui o vetor da população na geração seguinte. Caso contrário, o vetor da população
atual é mantido na próxima geração. Esta última operação é chamada seleção e o
procedimento de otimização continua até que algum critério de parada seja cumprido.
Portanto, dada uma população inicial e definidos os valores dos fatores de
mutação e cruzamento (F e CR), os processos de mutação, cruzamento e seleção são
repetidos sucessivamente até que o procedimento de otimização seja encerrado (veja
Figura 1.1 a seguir).
Iniciação
Avaliação
Definir F e CR
Repetir
Mutação
Cruzamento
Avaliação
Seleção
Até (critério de parada)
Figura 1.1 – Passos da ED
Na realidade o fator de mutação é responsável pelo mapeamento do domínio de
busca, ou seja, o fator F permite controlar uma convergência prematura. Já o fator de
cruzamento (CR) é responsável pela velocidade da convergência: quanto maior, mais
rápida será a convergência, segundo [7].
10
Diante do que foi exposto, o algoritmo da ED tradicional baseia-se numa ideia
muito simples, ou seja, uma busca global por meio da adição de vetores, que cria um
“sobrevivente” a partir do simples processo de seleção. Além disso, a ED é muito fácil
de ser implementada (poucas linhas de programação) em qualquer plataforma de
cálculo, contendo uma quantidade limitada de parâmetros a ser sintonizados (somente
NP, F e CR). Neste sentido, este método foi escolhido para servir de base para o
algoritmo de otimização proposto neste trabalho.
11
2
REVISÃO BIBLIOGRÁFICA
Este trabalho foi fundamentado a partir dos conceitos encontrados nas
referências bibliográficas que serão vistas a seguir: o atual modelo de comercialização
de energia elétrica no Brasil; método de otimização da Evolução Diferenciada; e o
método de aproximação da Superfície de Resposta.
2.1
SETOR ELÉTRICO BRASILEIRO
Até o ano de 1995, o Setor Elétrico Brasileiro (SEB) era conduzido por
empresas federais e estaduais, ou seja, o Estado mantinha o monopólio do setor. Desde
o início da década de 90, ocorreram dois grandes processos de reforma do SEB que
contribuíram para o formato do modulo atual do setor.
O início do primeiro processo de reforma do Setor Elétrico Brasileiro ocorreu
em 1993 com a Lei nº 8.631, que extinguiu a equalização tarifária vigente e criou os
contratos de suprimento entre geradores e distribuidores. Em 1995, a reforma foi
marcada pela criação do Produtor Independente de Energia e foi criado o conceito de
Consumidor Livre, a partir da promulgação da Lei nº 9.074 de 1995 [2, 3].
Em 1996, coordenado pelo Ministério de Minas e Energia (MME), foi
implantado o Projeto de Reestruturação do Setor Elétrico Brasileiro (Projeto RE-SEB).
Nele eram previstos: a necessidade de implantar a desverticalização das empresas de
energia elétrica, a partir das divisões dos segmentos de geração, transmissão e
distribuição; incentivar a competição nos segmentos de geração e comercialização; e
manter sob regulação os setores de distribuição e transmissão de energia elétrica,
considerados como monopólio naturais, sob regulação do Estado [8]. Este projeto foi
concluído em agosto de 1998.
12
Concomitante à criação do Projeto RE-SEB, foram implementadas mudanças de
cunho institucional para viabilizar alterações no setor elétrico, destacando-se: criação de
um órgão regulador (a Agência Nacional de Energia Elétrica - ANEEL); criação de um
operador para o sistema elétrico nacional (Operador Nacional do Sistema Elétrico ONS); e criação de um ambiente para a realização das transações de compra e venda de
energia elétrica (o Mercado Atacadista de Energia Elétrica - MAE).
Todas essas medidas eram condizentes com o objetivo de aumentar a eficiência
dos serviços elétricos, restaurar o equilíbrio financeiro e sanear a crise no setor (gerada
com a crise da década de 80) atraindo novos investimentos. Porém, essas medidas não
foram eficazes, pois a reforma não foi eficiente em atrair investimentos para a
necessária expansão do sistema, o que, aliado ao crescimento da demanda e ao
deplecionamento dos reservatórios, levaram ao racionamento de energia elétrica em
2001 [9].
A crise de abastecimento de energia elétrica ocorrida em 2001 gerou uma série
de questionamentos em relação ao rumo do SEB. Com isso, visando adequar o primeiro
modelo (anteriormente citado), foi implantado em 2002 o Comitê de Revitalização do
Modelo do Setor Elétrico, com o intuito de obter propostas de alterações no SEB [8].
Com a entrada do novo governo em 2003, houve o início do segundo processo
de reforma do Setor Elétrico Brasileiro, ou seja, uma nova reestruturação onde foram
focados, sob a ótica do governo, os seguintes objetivos: garantir a segurança de
suprimento de energia elétrica; promover a modicidade tarifária, através da contratação
eficiente de energia para os consumidores regulados; e promover a inserção social no
setor elétrico, em particular através dos programas de universalização de atendimento.
Com isso, sustentado pelas Leis nº 10.847 e 10.848, de 15 de março de 2004 e
pelo Decreto nº 5.163, de 30 de julho de 2004, o Governo Federal lançou as bases do
13
novo modelo para o SEB, em que foram definidas: a criação de uma instituição
responsável pelo planejamento do setor elétrico a longo prazo (a Empresa de Pesquisa
Energética - EPE); a criação de uma instituição com a função de avaliar
permanentemente a segurança do suprimento de energia elétrica (o Comitê de
Monitoramento do Setor Elétrico - CMSE); e a criação de uma instituição para dar
continuidade às atividades do MAE, relativas à comercialização de energia elétrica no
sistema interligado (a Câmara de Comercialização de Energia Elétrica - CCEE) [2, 8].
Na Tabela 2.1 estão representadas as diversas alterações que o Setor Elétrico
Brasileiro sofreu até chegar o modelo vigente.
Uma das grandes mudanças do novo modelo passa ser em relação à
comercialização de energia, a partir da instituição de dois tipos de ambientes para
celebração de contratos de compra e venda de energia: o Ambiente de Contratação
Regulada (ACR); e o Ambiente de Contratação Livre (ACL).
No ACR, as distribuidoras compram energia a longo prazo para o atendimento
dos seus consumidores cativos por meio de licitação na modalidade de leilões
conduzidos pela Agência Nacional de Energia Elétrica (ANEEL). Neste ambiente
participam os Agentes Vendedores (comercializadores,
geradores,
produtores
independentes ou autoprodutores) e os de Distribuição, que são os compradores de
energia elétrica. A contratação neste ambiente é formalizada através de contratos
bilaterais regulados, denominados Contratos de Comercialização de Energia Elétrica no
Ambiente Regulado (CCEAR).
14
Tabela 2.1 - Transformações dos modelos do SEB [8]
Modelo Antigo
Modelo de Livre Mercado
Novo Mode lo
(até 1995)
(1995 a 2003)
(a partir de 2004)
Financiamento através de recursos
públicos
Financiamento através de recursos
públicos e privados
Consumidores Cativos
Consumidores Livres e Cativos
Financiamento através de recursos
públicos e privados
Empresas divididas por atividade:
Empresas divididas por atividade:
geração, transmissão, distribuição,
Empresas verticalizadas
geração, transmissão, distribuição e
comercialização, importação e
comercialização
exportação.
Empresas predominantemente
Abertura e ênfase na privatização das Convivência entre Empresas Estatais
Estatais
Empresas
e Privadas
Competição na geração e
Competição na geração e
Monopólios - Competição inexistente
comercialização
comercialização
T arifas reguladas em todos os
segmentos
Preços livremente negociados na
geração e comercialização
Mercado Regulado
Mercado Livre
Planejamento Determinativo - Grupo
Coordenador do Planejamento dos
Sistemas Elétricos (GCPS)
Planejamento Indicativo pelo
Conselho Nacional de Política
Energética (CNPE)
Contratação : 85% do mercado (até
agosto/2003) e 95% mercado (até
dez./2004)
Contratação: 100% do Mercado
Sobras/déficits do balanço energético
rateados entre compradores
Sobras/déficits do balanço energético
liquidados no MAE
Consumidores Livres e Cativos
No ambiente livre: Preços livremente
negociados na geração e
comercialização. No ambiente
regulado: leilão e licitação pela menor
tarifa
Convivência entre Mercados Livre e
Regulado
Planejamento pela Empresa de
Pesquisa Energética (EPE)
Contratação: 100% do mercado +
reserva
Sobras/déficits do balanço energético
liquidados na CCEE. Mecanismo de
Compensação de Sobras e Déficits
(MCSD) para as Distribuidoras.
No ACL participam os Agentes de Geração, os de Comercialização, os
Importadores e os Exportadores de energia, e os Consumidores Livres, onde são
firmados contratos bilaterais entre os mesmos e a negociação é livre. Os envolvidos
neste ambiente de contratação celebram contratos do tipo PPA (Power Purchase
Agreement), que é um contrato de compra e venda de energia por um período
determinado com condições pré-estabelecidas de preços e volumes, firmadas entre as
partes [10].
Os Agentes de Geração, sejam concessionários de serviço público de Geração,
Produtores
Independentes
de
energia
ou
Autoprodutores,
assim
como
os
Comercializadores, podem vender energia elétrica nos dois ambientes, mantendo o
caráter competitivo da geração, e todos os contratos, sejam do ACR ou do ACL, são
registrados na CCEE e servem de base para a contabilização e liquidação das diferenças
15
no Mercado de Curto Prazo [8]. Uma visão geral da comercialização de energia é
apresentada Figura 2.1, a seguir.
Figura 2.1 - Ambientes de Comercialização [8]
Além dos dois ambientes de contratação mencionados anteriormente, existe o
Mercado de Curto Prazo (MCP). O mercado se resume na apuração das diferenças entre
a energia contratada e a energia verificada e posterior valoração dessas diferenças ao
Preço de Liquidação das Diferenças (PLD) ou preço spot. Exemplificando, para um
Agente Distribuidor apura-se a diferença entre o montante total contratado e o seu
consumo, e para um Agente Gerador apura-se a diferença entre o montante total
contratado e a geração, essas diferenças (podendo ser positivas ou negativas) são
valoradas pelo Preço de Liquidação das Diferenças que irão compor os pagamentos ou
recebimentos do Agente no âmbito da CCEE [8].
O PLD é determinado para cada patamar de carga e para cada submercado (Sul,
Sudeste/Centro-Oeste, Norte e Nordeste, totalizando quatro), tendo como base o custo
marginal de operação do sistema elétrico, limitado por um preço mínimo (PLD mínimo)
e um preço máximo (PLD máximo). Esses limites existem de acordo com legislação da
ANEEL, com validade entre a primeira e a última semana operativa de preços do ano.
Conforme descrito pela Câmara de Comercialização de Energia Elétrica, a
16
metodologia para determinação do PLD é operacionalizada através do programa
NEWAVE. Este programa é utilizado para modelar o planejamento da operação
hidrotérmica, a partir do cálculo do Custo Marginal de Operação (CMO) do setor
elétrico nacional que, no caso brasileiro, serve de base para a definição do preço spot ou
Preço de Liquidação das Diferenças (PLD).
O modelo de planejamento de operação hidrotérmica (NEWAVE) representa as
energias afluentes modeladas por um processo auto regressivo periódico de ordem p,
PAR(p), aproximando o comportamento estocástico das vazões naturais. Os custos
marginais de operação apresentam uma volatilidade decorrente da aleatoriedade das
vazões afluentes. Como os preços da energia no curto prazo são calculados pelo CMO,
essa característica influencia os resultados econômicos da geração e comercialização de
energia elétrica [3].
Não foi escopo deste trabalho a abordagem da metodologia de cálculo utilizado
para o planejamento de operação, tampouco como o programa NEWAVE gera os
valores dos CMO. Entretanto, o CMO mencionado tem uma relação proporcionalmente
inversa com o volume de água contida nos reservatórios. Com isso, há maior
probabilidade deste custo ser pequeno na maioria do tempo, desfavorecendo a atuação
de qualquer empreendimento somente no Mercado de Curto Prazo. Ou seja, haverá
maior probabilidade do custo operacional do empreendimento ser maior do que o PLD
(ou preço spot) durante a vida útil do projeto [2].
Para efeito de cálculo e conclusão desta dissertação, foi utilizado o banco de
dados disponível no próprio site da CCEE [8], que contém os últimos valores de CMO
com previsões equivalentes a cinco anos.
Resumindo, a nova estrutura do Setor Elétrico Brasileiro classifica as usinas
termelétricas como Agentes Geradores, permitindo, com isso, às mesmas atuarem tanto
17
no Ambiente de Contratação Regulada, quanto no Ambiente de Contratação Livre e no
Mercado de Curto Prazo. Dado este contexto de flexibilidade operacional é
imprescindível a obtenção de um modelo de avaliação de projetos que capture as reais
vantagens provenientes da flexibilidade operacional de uma planta de geração
termelétrica, transformando esse tipo de operação flexível em valor (receita) para o
empreendimento termelétrico [1].
2.2
ALGORITMOS EVOLUTIVOS
Os algoritmos evolutivos, segundo a referência [11], são baseados nos princípios
da Teoria da Evolução e Genética, onde só há sobrevivência dos indivíduos mais aptos.
Os primeiros trabalhos envolvendo métodos evolutivos são datados da década de 30,
quando sistemas evolutivos naturais passaram a ser investigados como algoritmos de
exploração de função objetivo multimodal.
Porém, somente com o crescimento do acesso a computadores, a partir da
década de 60, que se intensificaram os desenvolvimentos de Algoritmo Evolutivos.
Acreditava-se na evolução biológica como uma potencial ferramenta para o processo de
otimização, possibilitando a sua utilização em problemas de engenharia.
Os métodos evolutivos são da classe dos métodos heurísticos (classe de
algoritmos de otimização estocástica), que utilizam estratégias de busca não
determinísticas. Esses métodos não estão fundamentados em complexos conceitos
matemáticos e não utilizam o gradiente da função objetivo para determinação da direção
de procura no processo iterativo de otimização (características marcantes dos métodos
determinísticos).
De acordo com a referência [12], os algoritmos evolutivos têm uma vantagem
importante sobre outros tipos de métodos numéricos, dentre os quais são citados os mais
18
importantes:
• Eles só exigem informações sobre a função objetivo (ou multiobjetivo) em si,
que pode ser explícita ou implícita, dando flexibilidade para lidar com um amplo
espectro de problemas;
• Apresentam simplicidade, ou seja, podem ser modelados com poucas e simples
linhas de código;
• Apresentam adaptação relativamente fácil para problemas das mais diversas
áreas;
• Eles podem ser aplicados a problema que apresentam funções objetivo (ou
multiobjetivo) descontínuas e não diferenciáveis;
• Podem facilmente escapar de ótimos locais.
Alguns dos Algoritmos Evolutivos mais comuns são: Algoritmos Genéticos [14]
(Genetic Algorithms, AG); Programação Evolutiva [15] (Evolutionary Programming,
PE) e a Evolução Diferenciada [4] (Differential Evolution, ED), dentre outros. Além
desses, existem alguns métodos baseados na analogia com o comportamento social de
espécies, como por exemplo, o Método do Enxame de Partículas [16] (Particle Swarm
Optimization, PSO) e o Método da Colônia de Formigas [17] (Ant Colony Optimization,
ACO), dentre outros.
Importantes trabalhos mostram as diferença entre os principais métodos
heurísticos aqui citados. A referência [13] mostra que o método da ED tem excelente
performance, quando comparado com outros algoritmos testados (PSO e AG). O mesmo
apresentou ser robusto, ter maior velocidade de convergência e encontrou o ótimo
global em quase todas as execuções do método, quando aplicado a funções numéricas
de benchmark. A referência [18] mostrou que o método da ED tem uma melhor
19
convergência do que o método AG, quando estes foram aplicados a problemas de
engenharia de diferentes áreas.
Atualmente o algoritmo da ED tem sido utilizado em muitos casos práticos de
engenharia [11] e vem se tornando progressivamente mais popular (quando comparado
com os outros métodos heurísticos), pois apresenta boas propriedades de convergência,
sendo considerado de fácil compreensão. Na realidade, segundo [7], esta popularidade
se deve ao fato do método cumprir com cinco requisitos que normalmente os usuários
exigem:
• Habilidade em lidar com funções não lineares, não diferenciáveis e multimodais;
• Permite cálculos paralelos de funções objetivo, quando esta tem um alto custo
computacional;
• Facilidade de uso e compreensão, ou seja, poucas variáveis de controle para
orientar a minimização ou maximização;
• Variáveis de controle robustas e de fácil de escolha;
• Convergência consistente para o mínimo global em tentativas independentes e
consecutivas.
2.2.1
Evolução Diferenciada
O algoritmo da Evolução Diferenciada (ED) foi proposto por [4] em 1995 como
uma nova abordagem heurística. A principal motivação para o desenvolvimento deste
algoritmo foi a lenta taxa de convergência e a dificuldade na determinação dos
parâmetros do polinômio de Chebychev exibida pelo algoritmo híbrido denominado
Recozimento Genético (Genectic Annealing). Foi durante a resolução deste problema
que [4] decidiu modificar o algoritmo híbrido de forma a trabalhar com codificação de
ponto flutuante e com operações aritméticas. Neste sentido, foi desenvolvido o operador
20
de mutação diferencial, o qual fundamenta o algoritmo da ED.
Em maio de 1996, o sucesso da ED foi demonstrado no Primeiro Concurso
Internacional em Otimização Evolutiva (First International Contest on Evolutionary
Optimization, 1st ICEO) e na Conferência Internacional sobre Computação Evolutiva
(International Conference on Evolutionary Computation - CEC) [19]. Nestes eventos o
algoritmo da ED terminou em terceiro lugar geral e em primeiro lugar entre os
algoritmos evolutivos.
Em 1997, a algoritmo da ED tradicional também foi apresentado no Segundo
Concurso Internacional em Otimização Evolutiva (2nd ICEO) [20] e acabou como um
dos melhores entre os algoritmos concorrentes.
Apesar da ED utilizar conceitos emprestados da ampla classe dos Algoritmos
Evolutivos, os resultados experimentais apresentados mostraram que o desempenho da
ED é melhor do que muitos dos Algoritmos Evolutivos mais comuns [7, 21]. Neste
sentido a ED ganhou popularidade por demostrar ser um algoritmo robusto, de fácil
compreensão e com poucos parâmetros de controle a serem definidos, como o tamanho
da população (NP), o fator de mutação (F) e o fator de cruzamento (CR) [7].
Entretanto, segundo a referência [22], embora o algoritmo da ED tradicional
proposto por [4] (apresentado brevemente na introdução deste trabalho, seção 1.4)
pareça ser, à primeira vista, muito eficiente, o mesmo esconde uma limitação. Se por
algum motivo o algoritmo não tiver sucesso na geração de indivíduos descendentes que
superem a população atual correspondente, o processo de busca será repetido com os
mesmos parâmetros de mutação (F) e cruzamento (CR) e provavelmente falhará por
cair em uma condição de estagnação indesejada. Em outras palavras, a principal
desvantagem da ED tradicional ocorre quando as ações inerentes ao método não forem
suficientes para gerar novas soluções promissoras, ficando a busca pelo valor ótimo
21
fortemente comprometida.
Na última década ficou evidente que o sucesso dos processos de otimização do
método da ED depende da correta definição dos seus parâmetros NP (tamanho da
população), F (fator de mutação) e CR (fator de cruzamento). Esta observação é
enfatizada em vários trabalhos, como por exemplo, nas referências [23, 24].
O tamanho da população NP corresponde à quantidade de vetores (ou
indivíduos) presentes em cada geração, sendo que cada indivíduo passa pelos processos
de mutação, cruzamento e seleção.
Para um determinado indivíduo estes processos podem ser benéficos, enquanto
para outros podem resultar em um desperdício de esforço computacional. Portanto, uma
população muito pequena pode promover uma quantidade muito limitada de processos
bem sucedidos, enquanto uma população muito grande pode conter um número elevado
de processos ineficazes. Em outras palavras, se o valor de NP for muito pequeno pode
causar convergência prematura e se for muito grande pode causar estagnação.
Segundo [15], a definição do valor de NP é semelhante ao que ocorre em outros
algoritmos evolutivos, quando se considera a dimensão (variáveis de otimização) do
problema para estimar a quantidade de indivíduos.
A referência [7] propõe um valor para NP igual a dez vezes a dimensão do
problema. No entanto, um estudo feito por [25], além de não confirmar esta proposta,
foi mostrado que uma população de tamanho menor do que a dimensionalidade do
problema pode ser o ideal em muitos casos.
A definição de valores ideais para os fatores de mutação (F) e cruzamento (CR)
é considerada uma tarefa desafiadora, por não ser simples e nem intuitiva, entretanto, é
crucial para garantir o sucesso do algoritmo.
Neste sentido, vários estudos têm sido propostos na literatura. A partir de uma
22
análise empírica, [22] chegou à conclusão que o uso de F igual a 1 não é recomendado,
pois provoca uma diminuição significativa do poder de explorar o domínio de busca.
Analogamente, o mesmo autor recomenda não utilizar um CR igual a 1, uma vez que
poderá diminuir drasticamente a possibilidade de criar vetores experimentais diferentes,
causando uma estagnação.
As referências [7, 26] recomendam utilizar valores entre [0,5; 1] e [0,8; 1] para
os fatores de mutação e cruzamento, respectivamente. Nos trabalhos publicados por [27,
28], os valores de F e CR foram fixados em 0,9, obtendo sucesso nos testes realizados.
A análise empírica relatada em [28], mostra que, em muitos casos, adotar F e CR
maiores do que 0,6 conduzem a resultados com melhor desempenho.
É fácil observar que na literatura existem várias propostas de valores ideais para
os parâmetros F e CR, confundindo assim, os usuários do algoritmo da ED, por ser
muitas vezes conflitantes.
Nos últimos anos, tornou-se claro para a comunidade científica que, apesar do
método da ED ser um bom algoritmo, existem margens consideráveis de melhoria na
estrutura do algoritmo, especialmente quando aplicado a problemas reais e de
engenharia. Com base nestas considerações, recentemente várias propostas de
modificações as estruturas da ED foram realizadas, destacando-se as adaptações dos
parâmetros F e CR durante o processo evolutivo [11-13, 23, 24, 29-36].
Segundo a referência [38], destacam-se na literatura como estado da arte os
algoritmos jDE [11], SaDE [23], JADE [35], EPSDE [33], CoDE [31], MDE_pBX [34].
Cada um destes algoritmos são modificações do método da ED tradicional proposto por
[4], onde adotam diferentes estratégias para adaptar os parâmetros F e CR e ainda
propõe outras estratégias de mutação, diferentes daquelas indicadas inicialmente por [4].
Nos últimos três anos alguns autores propuseram modificações ao algoritmo
23
JADE [35], destacando-se os algoritmos wJADE-4 [36] e AdapSS-JADE [37], propostos
nos anos de 2012 e 2011, respectivamente. Estes algoritmos obtiveram performance
melhor do que os considerados como estado da arte. Neste sentido, o algoritmo de
otimização utilizado nesta dissertação foi baseado nas propostas indicadas nos trabalhos
[36, 37], as quais serão abordadas com maiores detalhes na seção 6.1 do capítulo 6.
Quanto ao Uso da Evolução Diferenciada na Otimização de Sistemas Térmicos
Por se mostrar simples e robusto, o algoritmo da Evolução Diferenciada vem
sendo aplicado na solução de diversos problemas de engenharia na última década, tais
como: projeto de sistemas [39], sistemas lineares [40], transferência de calor [41],
projeto manipulador de robô [42], etc.
O método da ED também vem sendo aplicado aos problemas de otimização
ligados a sistema de geração e cogeração. No trabalho desenvolvido por [43] a ED foi
utilizada para otimizar um sistema de cogeração segundo o menor custo de
implementação, utilizando o conceito de termoeconomia, obtendo grande sucesso.
A referência [44] aplicou o algoritmo evolutivo na otimização termoeconômica
paramétrica do sistema benchmark de cogeração denominado CGAM. Segundo o autor,
o algoritmo se mostrou poderoso, proporcionando resultados de acordo com estudos
anteriores.
Em [45] foi apresentado a utilização de um algoritmo evolutivo para otimizar o
custo de implementação de uma UTE em ciclo combinado, segundo o conceito de
termoeconomia. Este trabalho fez uso do software comercial GateCycle para resolver o
balanço de massa e energia da planta térmica. No que diz respeito ao algoritmo de
otimização, foi utilizado o pacote de programação GEATbx (Genetic and Evolutionary
Algorithm Toolbox) disponível no MATLAB, obtendo grande sucesso.
24
A referência [46] equacionou e definiu uma estratégia de compra de energia
elétrica no leilão público para uma Distribuidora, considerando-se as restrições impostas
pelo Setor Elétrico Brasileiro (SEB). Para encontrar a melhor estratégia foi utilizado o
algoritmo da ED tradicional, obtendo o sucesso esperado.
Quanto a Otimização do VPL de uma UTE, Segundo as Regras do SEB
Foram encontradas na literatura várias propostas para equacionar o VPL de um
empreendimento termelétrico, segundo as regras impostas pelo Setor Elétrico Brasileiro.
Entretanto, nestes estudos as características da térmica eram conhecidas e o objetivo era
identificar qual a melhor forma de investir em um empreendimento previamente
conhecido [1-3, 47-53].
A proposta desta dissertação foi exatamente inversa, ou seja, identificada uma
forma de investir (podendo ser a melhor ou não) segundo as regras e restrições impostas
pelo SEB, o objetivo foi determinar os melhores parâmetros do processo de geração de
uma UTE para o tipo de investimento escolhido. Este tipo abordagem é nova e singular,
onde as referências citadas serviram apenas de base para que o Valor Presente Líquido
deste estudo fosse modelado.
Nos capítulos 3 e 5 será abordado com maior detalhe a formulação do VPL de
um empreendimento termelétrico instalado no Brasil adotada nesta dissertação,
atribuindo um risco de apenas 5%. Esta formulação foi utilizada como função objetivo
para otimizar o ciclo de geração, que será apresentado no capítulo 4.
2.3
MÉTODO DA SUPERFÍCIE DE RESPOSTA
Modelos de superfície de resposta são frequentemente usados na substituição de
modelos físicos complexos, quando o cálculo da função objetivo é dispendioso e
25
quando se tem um grande número de variáveis de decisão. Este modelo é usado no lugar
da formulação original, a fim de reduzir o custo computacional em problemas de
otimização [55].
A técnica de aproximação é construída a partir de uma formulação analítica
adequada, usando dados disponíveis sobre a função objetivo que se deseja substituir.
Uma das formulações mais populares são as funções de base radial (RBF1) [56].
Inicialmente,
as RBF foram
desenvolvidas
para determinar soluções
aproximadas de equações diferenciais parciais. A referência [57] verificou que as
funções de base radial as RBF eram capazes de construir um esquema de interpolação
com propriedades favoráveis, tais como elevada eficiência, qualidade e capacidade de
lidar com dados dispersos, especialmente para os problemas de dimensão superior. A
referência [57] cita que as RBF apresentam qualidade de interpolação superior, quando
comparadas com outras técnicas de interpolação como, por exemplo, a transformada
Wavelet, a Expansão de Taylor e simulação sequencial gaussiana.
Várias aplicações em problemas de Engenharia Mecânica têm sido feitas nos
últimos anos, por exemplo, [5, 58, 59]. A referência [5] aplicou o modelo de superfície
de resposta num problema de otimização termoeconômica paramétrica do sistema
benchmark de cogeração CGAM utilizando um simulador termodinâmico e conclui que
o emprego deste método pode reduzir o número de avaliações da função objetivo em
muitos casos.
A grande vantagem do emprego desta técnica em funções objetivo implícitas
(como a deste trabalho) consiste em obter um modelo de aproximação analítico,
permitindo que a busca pelo ponto ótimo (durante um processo de otimização) seja mais
rápida do que no modelo original. Contudo, se um modelo de superfície de resposta for
1
Radial Basis Function
26
utilizado, o grau de fidelidade diminuirá, segundo a referência [5].
Com o intuito de diminuir o número de funções avaliadas (ou seja, diminuindo o
número de simulações do Thermoflow) durante o processo de otimização, neste trabalho
foi empregado o método de superfície de resposta formulado em RBF em substituição a
formulação original da função objetivo. Para contornar a perda de fidelidade por utilizar
uma função aproximada, foi utilizado o algoritmo H3 proposto por [55].
Assim, na seção 6.3.2 do capítulo 6 será apresentado o algoritmo de otimização
proposto neste trabalho, o ED2SR. Este algoritmo inicia-se fazendo uso da metodologia
de otimização proposta por [55], através do algoritmo H3, entretanto com algumas
modificações. E finaliza a otimização fazendo um intercâmbio entre o algoritmo H3
modificado e o método da Evolução Diferenciada Adaptativa.
27
3
COMERCIALIZAÇÃO DE ENERGIA ELÉTRICA NO SEB
Como visto no capítulo 2, seção 2.1, no atual modelo do SEB as relações
comerciais se estabelecem no Ambiente de Contratação Regulada e no Ambiente de
Contratação Livre. O primeiro inclui o mercado cativo de energia das distribuidoras,
onde estas são obrigadas a comprar energia elétrica de todas as geradoras (incluindo
termelétricas) participantes dos leilões com contratos de longo prazo. O segundo inclui
a livre negociação, onde os consumidores livres e os comercializadores têm liberdade
para negociar a compra de energia elétrica, estabelecendo volumes, preços e prazos de
suprimento [8].
Os
Agentes
de
Geração,
Produtores
Independentes
de
energia
ou
Autoprodutores, assim como os Comercializadores, podem vender energia elétrica nos
dois ambientes, mantendo o caráter competitivo da geração e todos os contratos, sejam
no ACR ou do ACL, são registrados na CCEE e servem de base para a contabilização e
liquidação das diferenças no MCP [8].
3.1
COMERCIALIZAÇÃO NO ACR
A comercialização de energia elétrica neste ambiente ocorre por meio de leilões.
Todos os agentes vencedores (geradores, importadores e comercializadores) podem
vender a todos os distribuidores, onde a estes compradores há obrigação de adquirir
energia sempre através de leilão.
O art. 27 do decreto n° 5.163, de 30 de julho de 2004 [2], estabelece que os
vencedores dos leilões de energia proveniente de empreendimentos de geração novos ou
existentes deverão formalizar um contrato bilateral denominado Contrato de
Comercialização de Energia Elétrica no Ambiente Regulado (CCEAR), celebrado entre
28
cada agente vendedor e todos os agentes de distribuição.
Os CCEAR são especificados por meio dos editais publicados para cada leilão,
contendo cláusulas e condições fixas, como preço da energia, submercado de registro do
contrato e vigência de suprimento, os quais não são passíveis de alterações bilaterais por
parte dos agentes envolvidos.
Após a assinatura pelos agentes vendedores e compradores, os CCEAR são
registrados pela CCEE no Sistema de Contabilização e Liquidação (SCL), para que
possam ser considerados no processo de contabilização e liquidação financeira.
Existem duas modalidades de CCEAR: por quantidade de energia elétrica; e por
disponibilidade de energia elétrica.
No CCEAR por quantidade, os riscos hidrológicos da operação energética são
assumidos integralmente pelos geradores, cabendo a eles todos os custos referentes ao
fornecimento da energia contratada, devendo existir mecanismos específicos para o
rateio dos riscos financeiros decorrentes de diferenças de preços entre submercados e
eventualmente impostos aos agentes de distribuição que celebrarem contratos nessa
modalidade.
No CCEAR por disponibilidade, os custos decorrentes dos riscos hidrológicos e
eventuais exposições financeiras no Mercado de Curto Prazo, positivas ou negativas,
serão assumidos pelas distribuidoras, com repasse ao consumidor final, conforme
mecanismo definido pela ANEEL.
No Brasil, a contratação de energia por quantidade se aplica somente aos
empreendimentos de geração hidrelétrica e tem duração de 30 anos [47]. As usinas
termelétricas celebram contratos CCEAR na modalidade disponibilidade.
Como o espoco deste trabalho foi verificar o Valor Presente Líquido de uma
UTE flexível, não foi abordada a metodologia da receita (remuneração) do gerador
29
quando este celebra contrato na modalidade por quantidade.
A seguir, será visto a metodologia de cálculo da remuneração da UTE adotada
neste trabalho, quando esta celebra contratos na modalidade disponibilidade.
3.1.1
Contratos de Disponibilidade
Quando da participação nos leilões, os empreendimentos termelétricos (de ciclo
aberto ou combinado a gás natural) celebram contratos com os agentes distribuidores na
modalidade disponibilidade de energia elétrica. Neste caso, as distribuidoras “alugam” a
UTE, que deve estar disponível sempre que for chamada a despachar. Em troca, a usina
recebe uma receita fixa que deve ser suficiente para cobrir seus custos fixos, remunerar
seu investimento e garantir a lucratividade desejada.
Nesta modalidade de contratação, a liquidação das diferenças entre a energia
produzida (pela UTE) e contratada (pela distribuidora) é de responsabilidade das
distribuidoras, que passam a ser responsáveis por qualquer transação no MCP. Se a
usina produz mais energia do que a quantidade que foi contratada, este excedente
pertence ao distribuidor, que tem a prerrogativa de revender esta diferença no MCP. Por
outro lado, se a energia entregue pelo gerador for menor do que a energia contratada, o
distribuidor deverá comprar esta diferença no mercado spot. Ou seja, como descrito na
seção 3.1, o risco hidrológico é assumido pelas distribuidoras.
Em suma, a remuneração da UTE, num contrato por disponibilidade, é composta
por uma receita fixa e uma parcela variável, realizada em pagamentos mensais. A
receita fixa cobre os custos de combustível e de operação e manutenção associados à
inflexibilidade da UTE e demais custos, como: remuneração do investimento, custo de
conexão e uso dos sistemas de transmissão e distribuição, seguros e garantias da usina,
tributos e encargos necessários ao cumprimento do contrato. A parcela variável mensal
30
refere-se ao custo de produção da geração variável (aquela além da inflexibilidade
operativa), dada pelo produto do Custo Variável Unitário (CVU), em R$/MWh, pela
geração variável, em MWh. Este CVU deve cobrir todos os custos de operação da
termelétrica, exceto aqueles cobertos pela receita fixa.
Com tudo isso exposto, a remuneração de uma UTE e o pagamento da
distribuidora à primeira são dados da seguinte forma, respectivamente [47]:
 RF ⋅ E ⋅ ht + (CVU OP − CVU R ) ⋅ GT ⋅ ht , se PLDt ≥ CVU OP
RCCEAR = 
 RF ⋅ E ⋅ ht , se PLDt < CVU OP
(3.1).
RF ⋅ E ⋅ ht + CVU OP ⋅ GT ⋅ ht − PLDt ⋅ (GT − E ) ⋅ ht , se PLDt ≥ CVU OP
PCCEAR = 
RF ⋅ E ⋅ ht + PLDt ⋅ E ⋅ ht , se PLDt < CVU OP
(3.2).
onde,
RF : receita fixa do contrato, em R$/MWh;
E: montante contratado no leilão, em MW;
GT : geração total da UTE, em MW;
PLDt : preço spot no período t, em R$/MWh;
CVU OP : custo variável unitário de caráter operacional1, em R$/MWh;
CVU R : custo variável unitário real da UTE, em R$/MWh;
ht : número de horas no período t.
A equação (3.1) mostra a receita da UTE quando esta celebra contratos no
ambiente regulado. Quando o custo operacional for menor do que o Preço de
Liquidação das Diferenças em um determinado período, a termelétrica entra em
operação e fornece toda a energia gerada neste período ( GT ) ao distribuidor. Ainda
nesta situação, a térmica assume o risco do custo do combustível no período t não ser
1
Corresponde ao Custo Variável Unitário declarado a EPE, o qual é utilizado pelo ONS para definir
quando do despacho da térmica. Maiores informações, ver seção 3.1.4.3 desta dissertação.
31
exatamente aquele declaro ao EPE, desta forma, o custo real de operação da termelétrica
( CVU R ) pode ser maior do que o previsto ( CVU OP ) e onerar a receita da UTE.
A equação (3.2) mostra o desembolso de uma distribuidora para adquirir a
energia contratada E. Conforme descrito anteriormente, a distribuidora sempre pagará
uma parcela fixa referente ao “aluguel” da térmica ( RF ⋅ E ⋅ ht ). Entretanto, caso o custo
operacional da UTE seja menor do que o Preço de Liquidação das Diferenças em um
período t, então a distribuidora paga o custo operacional ( CVU OP ⋅ GT ⋅ ht ) a térmica,
adquirindo toda a energia gerada. Nesta mesma equação, nota-se que a distribuidora tem
a prerrogativa de liquidar a diferença entre o que foi comprado ( GT ) com o que foi
contratado (E) no Mercado de Curto Prazo ao preço do PLD. Assumindo, desta forma, o
ônus se esta diferença for negativa.
No caso em que o custo operacional for maior do que o Preço de Liquidação das
Diferenças em um determinado período, o UTE não opera, recebe o seu “aluguel” e a
distribuidora adquire a energia contratada E no MCP ao preço do PLD.
Assim, a equação (3.1) corresponde à receita da UTE quando esta celebra
contratos no ambiente regulado. Uma das parcelas mais importante desta equação é a
Receita Fixa ( RF ), sendo o seu valor aquele registrado no momento em que os agentes
envolvidos assinam e registram o contrato no âmbito do ACR após o leilão. Se esta
parcela receber um valor alto, beneficiará o vendedor, entretanto, o mesmo perderia
competitividade no leilão, pois estaria onerando o custo do comprador. Desta forma,
para quantificar esta Receita Fixa ( RF ) foi levado em consideração as regras impostas a
uma térmica, quando da sua participação em Leilões de Energia Nova.
32
3.1.2
Leilões De Energia Nova
Os leilões são a principal forma de contratação de energia no Brasil. Por meio
desse mecanismo, concessionárias, permissionárias e autorizadas de serviço público de
distribuição (distribuidores) de energia elétrica do SIN garantem o atendimento à
totalidade de seu mercado no ACR.
O Leilão de Energia Nova (LEN), realizado pela CCEE, por delegação da
ANEEL, é o que possibilita a contratação a longo prazo da energia de futuros
empreendimentos de geração, sendo utilizado o critério de menor tarifa para definir os
vencedores do certame [8].
O LEN tem como finalidade atender ao aumento de carga das distribuidoras.
Este leilão pode ser de dois tipos: A-5 e A-3 (usinas que entram em operação comercial
em até cinco e três anos, respectivamente).
Os leilões de contratos de disponibilidade de energia elétrica (foco deste estudo)
funcionam para as distribuidoras como um contrato de opção. Como visto no item
anterior, a distribuidora paga um “aluguel” (receita fixa) e tem a opção, e não a
obrigação, de comprar a energia da UTE pelo custo variável de operação declarado pelo
gerador. Ou seja, se o preço spot for menor que o custo declarado, a distribuidora não
exerce a opção de compra e adquire a energia no MCP. Caso contrário, exerce a opção
reembolsando o gerador pelo preço de exercício, podendo ainda vender no MCP a
diferença entre o montante produzido pela usina e o montante contratado [47], conforme
a equação (3.2).
Para o distribuidor, o custo da energia contratado por disponibilidade é função
do preço de venda do gerador (receita fixa) e do preço de exercício declarado pelo
empreendedor (CVU). Quanto maior o valor de CVU, menor é a probabilidade de
despacho da usina, implicando maiores custos de compra de energia no MCP. Esta
33
observação há de ser levada em consideração no momento do leilão, para que os
empreendimentos de geração (vendedores de energia) sejam comparados em uma
mesma base [60].
No Brasil, o critério para priorização de projetos de geração em LEN é o
tradicional método do Índice de Custo Benefício (ICB). Ou seja, o ICB é utilizado para
a
ordenação
econômica
de
empreendimentos
de
geração
termelétrica
e,
consequentemente, como critério de contratação por meio de contratos de
disponibilidade de energia elétrica, sob o regime de autorização por 15 anos.
3.1.3
O Índice de Custo Benefício
Em um sistema de geração predominantemente hidroelétrico como o do Brasil, o
benefício energético da operação integrada de um empreendimento de geração
(hidroelétrica ou termelétrica) pode ser avaliado pela garantia física da usina, [60].
No caso das térmicas, a garantia física poderia ser calculada pela máxima
potência instalada que a UTE consegue gerar continuamente (valor conhecido como
disponibilidade da usina, que será visto mais adiante). Todavia, térmicas com diferentes
custos operativos (CVU) tem diferentes frequências de despacho, contribuindo de
maneira diferente para a confiabilidade do sistema elétrico [47]. Com isso, quanto maior
for o custo variável menor será a garantia física da usina atribuída pela EPE como
proporção de sua potência total instalada, uma vez que quanto maior for custo variável
menor será a probabilidade desta usina vir a ser despachada pelo ONS.
Para térmicas totalmente inflexíveis (sempre operando, independente do seu
custo operacional) a garantia física é igual a sua disponibilidade. À medida que a UTE
se torna flexível, a geração desta parte fica sujeita à opção do distribuidor, como visto
nos itens anteriores.
34
Sendo assim, a garantia física é definida como um percentual da disponibilidade,
podendo ser igual ou menor que a mesma. A magnitude deste percentual varia com o
nível de inflexibilidade e CVU declarados.
Na prática a garantia física é calculada à época do leilão, aplicando-se a
metodologia definida na Portaria MME 258, de 28 de julho de 2008, e o critério
definido na Resolução n° 9, de 28 de julho de 2008, do Conselho Nacional de Política
Energética (CNPE).
Assim, o ICB é definido pela razão entre o custo global do empreendimento e o
seu benefício energético. Como visto, o benefício energético é a garantia física da UTE
e o custo global compreende três parcelas: receita fixa anual; do custo de operação
(COP); e custo econômico de curto prazo (CEC).
A parcela da receita fixa anual é aquela que corresponde ao pagamento fixo da
distribuidora, devendo ser suficiente para remunerar o investimento da térmica,
incluindo todos os custos fixos inerentes ao empreendimento e o lucro desejado. Como
visto nos itens anteriores, esta parcela corresponde ao “aluguel” da usina por parte da
distribuidora.
O COP, expressos em R$/ano é o valor esperado anual do reembolso do custo de
operação acima da inflexibilidade1 operativa. Esta parcela (assim como a garantia física)
é também função do nível de inflexibilidade no despacho da UTE e do CVU declarado à
época do LEN, os quais determinam sua condição de despacho em função dos CMO
futuros observados no SIN. Portanto, trata-se de uma variável aleatória cujo valor
esperado é calculado com base em uma amostra de valores de CMO divulgados pelo
EPE [61, 62].
O CEC, expressos em R$/ano, corresponde ao valor esperado anual das
1
Inflexibilidade corresponde ao montante de potência mínima que deve ser obrigatoriamente despendido
e que não está sujeito, portanto, à regra de despacho do ONS.
35
liquidações das diferenças no mercado de curto prazo das distribuidoras. Esta parcela
contabiliza os gastos e os ganhos da distribuidora no MCP, feita com base no CMO,
este último limitado aos preços de PLD mínimo e máximo. O Valor do CEC também é
função do nível de inflexibilidade no despacho da UTE e do CVU declarado pela usina
à época do LEN. Trata-se, portanto, de uma variável aleatória, cujo valor esperado é
calculado com base na mesma amostra de valores de CMO utilizada no cálculo da
parcela COP [60, 62].
Portanto, o ICB é calculado de acordo com a seguinte expressão:
ICB =
R F ,ano
COP + CEC
+
8760 ⋅ QL
8760 ⋅ GF
(3.3).
Sendo,
QL = x ⋅ GF
(3.4).
onde,
RF ,ano : receita fixa anual, em R$/ano;
QL: quantidade de lotes1 ofertados no leilão, MW médios;
COP: custo de operação, em R$/ano;
CEC: custo econômico de curto prazo, em R$/ano;
GF: garantia física, em MW médios;
8760: número de horas no ano;
x: fração da GF destinada ao ACR, valor entre [ xmin , 1].
É importante observar, que no caso de um empreendimento que deseja destinar
apenas uma fração (x) de sua garantia física ao ACR, sendo o restante para uso próprio
ou comercialização no ACL, o ICB, conforme equação (3.4), será calculado com base
nesta fração. Com isso, se a receita fixa não for diminuída na mesma proporção que a
1
Nos Leilões de Energia Nova o montante negociado é dividido em lotes, podendo esta divisão ser em
kW ou MW, dependendo do leilão.
36
garantia física a usina perderá competitividade durante o leilão.
3.1.4
Metodologia de Cálculo do ICB
Nesta seção será abordada a metodologia de cálculo de cada uma destas parcelas
que compõe o ICB, equação (3.3), e outros conceitos e variáveis inerentes a UTE,
como: disponibilidade; inflexibilidade; custo variável unitário (custo de operação acima
da inflexibilidade); custos marginais de operação; e custos marginais de operação
limitados.
3.1.4.1 Disponibilidade
A disponibilidade é a máxima geração média mensal da UTE, em MW. Ou seja,
corresponde à potência instalada da usina, descontando as taxas de paradas por falha ou
manutenção, conforme equação a seguir:
Disp = Pot líq ⋅ FCap max ⋅ (1 − TEIF ) ⋅ (1 − IP )
(3.5).
onde,
Potlíq : potência líquida instalada da usina em MW;
FCapmax : percentual da potência instalada que a usina gera continuamente;
TEIF: taxa média de indisponibilidade forçada;
IP: taxa de indisponibilidade programada.
3.1.4.2 Inflexibilidade
A inflexibilidade da usina, em MW, é o montante de potência mínima que deve
ser obrigatoriamente despendido e que não está sujeito, portanto, à regra de despacho do
37
ONS. A inflexibilidade da usina pode ser oriunda de razões tecnológicas como, por
exemplo, quando um equipamento específico exigir uma carga mínima de potência para
seu funcionamento adequado, ou proveniente de motivos econômicos como, por
exemplo, quando o contrato de fornecimento de gás natural for do tipo “take or pay”,
sendo economicamente inviável mantê-la desligada. Desta maneira, o nível de
inflexibilidade é calculado da seguinte forma:
Inflex = Inflex% ⋅ Disp
(3.6).
onde,
Inflex% : proporção da potência que é inflexível, em %;
Disp: Disponibilidade da usina, MW, equação (3.5).
O investidor deve declarar à EPE a proporção de sua potência que é inflexível
( Inflex% ) e deve cobrar os gastos relativos à inflexibilidade na sua receita fixa anual.
3.1.4.3 Custo Variável Unitário - CVU
O CVU é o valor do custo variável, para cada MWh gerado pela usina, expresso
em R$/MWh, informados pelo gerador, necessários para cobrir todos os custos de
operação da usina, exceto os já cobertos pela receita fixa. Este custo é composto de duas
parcelas: a primeira vinculada ao custo do combustível, e a segunda vinculada aos
demais custos variáveis [63].
CVU = C COMB + C OeM
(3.7).
onde,
CCOMB : custo do combustível destinada à geração flexível, em R$/MWh;
COeM : demais custos variáveis da usina, em R$/MWh.
De acordo com o disposto nas Portarias MME nº. 42, de 1º de março 2007 e nº.
38
46, de 9 de março de 2007 (com a redação dada na Portaria nº 175, de 16 de abril de
2009) destaca-se o cálculo de dois valores de CVU:
• O primeiro é de caráter estrutural (denominado CVU E ), detalhado no art. 5º da
Portaria MME nº 46/2007, sendo destinado ao cálculo dos parâmetros energéticos
como: Garantia Física (GF), Custo de Operação (COP) e Custo Econômico de Curto
Prazo (CEC), ou seja, o ICB, equação (3.3);
• O segundo tipo de CVU é aquele de caráter operacional (denominado CVU OP ),
definido no art. 3º da Portaria MME nº. 42, de 1º de março 2007, o qual definirá o
despacho do empreendimento, pelo ONS, acima da sua inflexibilidade operativa.
Ambos os CVU ( CVU E e CVU OP ) se diferenciam quanto ao cálculo do custo
do combustível, CCOMB .
O custo do combustível para a formação do CVU estrutural ( CVU E ) é calculado
da seguinte forma:
C COMB = i ⋅ e0 ⋅ PREF
(3.8).
onde,
i: fator de conversão, transforma o preço do combustível, em 106Btu/MWh;
e0 : média anual da taxa de câmbio divulgada pelo Banco Central, em R$/US$;
PREF : preço de referência, expectativa de preço futuro do combustível, em
US$/106Btu.
O fator de conversão, i, informado pelo agente no AEGE, constará no CCEAR,
permanecendo invariável por toda a vigência do contrato (embora este fator seja
dimensionalmente homogêneo ao consumo específico, o “heat rate”, ele não o
39
representa guardando apenas uma relação com o mesmo). Este fator deve ser declarado
com quatro casas decimais.
O preço de referência do combustível corresponde a expectativa de preço futuro
para o período de dez anos, no qual se inclui o ano de realização do leilão. Este valor é
estimado com base em projeções de combustíveis equivalentes, no cenário de referência
publicado pela Energy Information Administration (EIA) [64] no Annual Energy
Outlook (AEO), conforme metodologia descrita em Nota Técnica da EPE [63].
O custo do combustível para a formação do CVU operacional ( CVU OP ), o qual
está vinculado ao despacho do empreendimento pelo ONS, será determinado pelo preço
spot do gás natural Henry Hub, segundo a EIA [62], sendo os valores do fator de
conversão (i) e do custo de operação e manutenção variável ( COeM ) os mesmos
adotados para o CVU de caráter estrutural.
Resumindo, os custos variáveis unitários CVU E e CVU OP se diferenciam tanto
no valor do preço de combustível, em US$/106Btu, quando no valor da taxa de câmbio,
em R$/US$.
3.1.4.4 Custos Marginais de Operação - CMO
O CMO representa o custo de operação do sistema hidrotérmico brasileiro, para
gerar 1 MWh. Ou seja, se o custo de operação de uma UTE for menor do que o CMO,
então a térmica é despachada, caso contrário, não.
Para realizar o planejamento da operação do SIN é necessário realizar projeções
dos valores de CMO. Estes valores são obtidos de resultados de uma simulação da
operação mensal do SIN, com auxílio do modelo NEWAVE1. Com o resultado desta
1
Programa utilizado para modelar o planejamento de operação hidrotérmica. Este programa aproxima o
comportamento estocástico das vazões naturais e gera distintos cenários de custo operacional do SEB.
40
simulação, obtém-se uma planilha de valores de CMO com 2.000 simulações para cada
mês ao longo de 5 anos (totalizando 60 meses), para cada um dos submercados1
considerados, a qual é publicada pela EPE [62].
Para efeito do cálculo tanto do ICB [60], equação (3.3), quanto do VPL foi
adotado apenas o critério de despacho por ordem de mérito (por razões energéticas) das
usinas termelétricas. Lembrando que há outros critérios como: razões elétricas e
segurança energética.
Assim, a regra de despacho mensal ocorre da seguinte maneira: quando o CVU
for inferior ao CMO, a UTE fica despachada no limite de sua disponibilidade, equação
(3.5); caso contrário, a usina gera o equivalente à sua inflexibilidade, equação (3.6). Em
termos matemáticos, escreve-se que:
CMOs ,c ,m ≥ CVU E ∴ Gc ,m = Dispm

CMOs ,c ,m < CVU E ∴ Gc ,m = Inflexm
(3.9).
onde,
s: corresponde ao submercado (1 a 4);
c: corresponde ao índice de cada cenário hidrológico (1 a 2000);
m: corresponde ao índice de cada mês (1 a 60);
CMOs ,c ,m : custo marginal de operação, em R$/MWh;
CVU E : custo variável unitário de caráter estrutural, em R$/MWh;
Gc ,m : geração da térmica, para cada cenário e mês, MW;
Inflexm : inflexibilidade de despacho da UTE, para cada mês, em MW;
Dispm : disponibilidade (geração máxima mensal) da usina, em MW.
Resumindo para cada submercado existe uma projeção de 60 meses (cinco
1
Sul, Sudeste/Centro-Oeste, Norte ou Nordeste. Havendo, desta forma, quatro submercados.
41
anos), sendo que para cada mês existem 2.000 possíveis cenários hidrológicos. Ou seja,
para cada um dos 60 meses existem 2.000 valores de CMO.
3.1.4.5 COP e CEC
O COP leva em conta o gasto adicional da térmica, considerada como um todo,
quando gera acima de sua inflexibilidade declarada. Este gasto compreende o custo
adicional do combustível propriamente dito e os custos incrementais de operação e
manutenção.
A média destes gastos em cada mês e para cada cenário de CMO, quando
multiplicado pelo número de horas em 12 meses, fornece o valor anual do COP
(R$/ano) da seguinte forma:
∑∑ CVU ⋅ (G
m
COP =
c
i =1 j =1
E
c ,m
− Inflexm )
m⋅c
⋅ 8760
(3.10).
onde,
c: corresponde ao índice de cada cenário hidrológico (1 a 2000);
m: corresponde ao índice de cada mês (1 a 60);
CVU E : custo variável unitário de caráter estrutural, em R$/MWh;
Gc ,m : geração da térmica, MW, equação (3.9);
Inflexm : inflexibilidade de despacho da UTE, para cada mês, em MW;
8760: número de horas no ano.
O CEC, como já mencionado na seção anterior, reflete os “ganhos” e as “perdas”
obtidos no MCP do distribuidor.
Assim, como o COP, o CEC é calculado para a térmica como um todo, para cada
mês e para cada uns dos possíveis cenário hidrológicos. Independentemente do valor do
42
seu CVU, a diferença entre a garantia física e a geração despachada da usina é valorada
pelo CMO limitado ao PLD mínimo e PLD máximo. Com isso, é calculada uma média
mensal de todos os possíveis “ganhos” e “perdas”, que quando multiplicado pelo
número de horas em 12 meses, fornece o valor anual do CEC, como segue:
m
c
∑∑ CMO
CEC =
i =1 j =1
*
s ,c ,m
⋅ (GF − Gc ,m )
m⋅c
⋅ 8760
(3.11).
onde,
c: corresponde ao índice de cada cenário hidrológico (1 a 2000);
m: corresponde ao índice de cada mês (1 a 60);
CMO s*,c ,m : custo marginal de operação limitado, em R$/MWh;
GF: garantia física da UTE, em MW;
Gc ,m : geração da térmica, MW, equação (3.9);
8760: número de horas no ano.
3.1.4.6 Receita Fixa
A receita fixa é parte integrante tanto remuneração de uma UTE, equação (3.1),
quanto da formulação do ICB, equação (3.3). Maximizar a receita fixa sem levar em
consideração as restrições impostas no ACR, tornará o projeto termelétrico menos
competitivo durante o LEN. Neste sentido, foi adotado para este trabalho como receita
fixa o equacionamento do ICB. Com isso, substituindo a equação (3.4) em (3.3) e
colocando a receita fixa anual em evidência, obtém-se:
RF ,ano = x ⋅ (8760 ⋅ GF ⋅ ICB − COP − CEC )
onde,
ICB: índice de custo benefício, em R$/MWh;
(3.12).
43
COP: custo de operação, em R$/ano;
CEC: custo econômico de curto prazo, em R$/ano;
GF: garantia física, em MW médios;
8760: número de horas no ano;
x: fração da Garantia Física destinada ao ACR, valor entre [ xmin , 1].
Dividindo a equação (3.12) pelo valor da garantia física (GF) e por 8.760 horas
no ano, obtém a receita fixa em R$/MWh. Ao substituir este resultado na equação (3.1)
é obtido à remuneração da UTE no âmbito do ACR.
3.2
COMERCIALIZAÇÃO NO ACL
Nesse ambiente, a comercialização de energia elétrica é feita com contratos
bilaterais ou com negociação de curto prazo (através do mercado spot). Esses contratos
são resultados da negociação entre agentes do setor elétrico, que são grandes
consumidores de energia, distribuidores e geradores.
No ACL podem participar os agentes de comercialização, geração, de
exportação, de importação, consumidores livres e consumidores especiais. Dessa
maneira, ficam de fora as distribuidoras, que apenas comercializam os seus excedentes
no MCP.
Neste trabalho, o montante de energia comercializada no ACL foi o excedente
do que foi ofertado no ACR. Ou seja, de acordo com as equações (3.4), parte integrante
da equação (3.3), foi comercializado no ACR um montante x da Garantia Física da
UTE. Desta forma, o excedente correspondeu a (1 - x) da Garantia Física.
Desta forma, se, por exemplo, um gerador tem 90% ( x = 0,9) da sua garantia
física direcionada ao mercado cativo (ACR), esse empreendedor comercializa os 10%
(1 − x = 0,1) restantes no mercado livre (ACL).
44
Uma UTE, quando atua no ACL, faz um papel parecido com os dos geradores,
ou seja, caso a mesma gere acima da fração (1 - x) negociada, esta pode comercializar o
excedente no MCP, sendo que a UTE não tem a obrigação (tem a opção) de gerar. Ou
seja, se não for economicamente vantajoso gerar a energia contratada, a térmica pode
comprar a fração (1 - x) da energia contratada no mercado spot (Mercado de Curto
Prazo) e honrar seu contrato bilateral no ambiente livre.
Esse contrato é do tipo PPA (Power Purchase Agreement), que é um contrato de
compra e venda de energia por um período determinado com condições préestabelecidas de preços e volumes (limitado ao montante da fração), firmadas entre as
partes envolvidas [10].
Diante desta flexibilidade operacional a remuneração de uma UTE neste
ambiente e os pagamentos dos consumidores livres são dados da seguinte forma,
respectivamente:
R PPA
 PPPA ⋅ E PPA ⋅ ht + PLDt ⋅ (G PPA − E PPA ) ⋅ ht −

=
G PPA ⋅ CVU R ⋅ ht , se PLDt ≥ CVU OP
 P ⋅ E ⋅ h − E ⋅ PLD ⋅ h , se PLD < CVU
PPA
t
t
t
OP
 PPA PPA t
Pg PPA = PPPA ⋅ E PPA ⋅ ht
onde,
PPPA : preço da energia contratada no PPA, em R$/MWh;
E PPA : energia contratada no PPA, em MW;
GPPA : geração total no PPA da UTE, em MW;
PLDt : preço spot de energia no período t, em R$/MWh;
CVU OP : custo variável unitário de caráter operacional, em R$/MWh;
CVU R : custo variável unitário real da UTE, em R$/MWh;
ht : número de horas no período t.
(3.13).
(3.14).
45
Neste ambiente os empreendimentos termelétricos assumem os riscos
hidrológicos. Na equação (3.13) pode notar que a UTE sempre honrará seu contratado
PPA comprando no MCP, correndo o risco de obter uma remuneração negativa quando
o valor de PLD for maior do que o preço da energia contratada ( PPPA ).
Com isso, é necessário averiguar se um projeto exposto ao risco gera maiores
receitas do que aquele que somente negocia no ACR, além de verificar qual é o projeto
termelétrico que melhor se adequa ao risco hidrológico.
3.3
MERCADO DE CURTO PRAZO
Como visto nas seções (3.1) e (3.2), todos os contratos de compra e venda de
energia, quando celebrados nos ACR e ACL, devem ser registrados na CCEE, que
realiza a medição dos montantes efetivamente produzidos e consumidos por cada
agente. As diferenças apuradas, positivas ou negativas, são contabilizadas para posterior
liquidação financeira no mercado de curto prazo (MCP) e valoradas ao preço de
liquidação das diferenças (PLD), ou preço spot.
Assim, o MCP pode ser definido como o segmento da CCEE onde são
contabilizadas as diferenças entre os montantes de energia elétrica contratada pelos
agentes e os montantes de geração e de consumo efetivamente verificados e atribuídos
aos respectivos agentes.
No MCP não existem contratos, ocorrendo à contratação multilateral, conforme
as regras de comercialização [8].
46
4
USINA TERMELÉTRICA
Neste capítulo será apresentada a planta termelétrica que foi otimizada neste
trabalho, segundo a nova estrutura do Setor Elétrico Brasileiro. Assim, será detalhado
como foi feita a comunicação entre o simulador termodinâmico (pacotes de programa
Thermoflow) e a plataforma de calculo escolhido (Matlab) durante o processo iterativo
de otimização.
4.1
PLANTA DE GERAÇÃO: CICLO COMBINADO
O ciclo combinado é o processo de produção de energia elétrica utilizando
turbinas a gás e turbinas a vapor associadas em uma única planta. O combustível é
queimado em um uma turbina a gás e o calor contido nos gases de exaustão produz
vapor em uma caldeira de recuperação. Este vapor aciona uma turbina a vapor de
condensação. Tanto a turbina a gás quanto a turbina a vapor acionam geradores para
produção de energia elétrica.
A Figura 4.1 mostra um ciclo combinado típico. A parte tracejada corresponde à
turbina a gás, também conhecido como ciclo Brayton. Neste ciclo, a expansão dos gases
resultantes da queima do combustível (gás natural, por exemplo) aciona a turbina, que
está diretamente acoplada ao gerador e, desta forma, a potência mecânica é
transformada em potência elétrica.
A outra parte da Figura 4.1 (fora do tracejado) é conhecida como ciclo Rankine.
Neste ciclo, vapor é gerado a partir da água que circula por uma extensa rede de tubos
que revestem a parede da caldeira recuperadora, absorvendo o calor dos exaustos do
ciclo Brayton. O vapor movimenta as pás da turbina a vapor, fazendo com que seu rotor
gire, gerando potência mecânica. Através de um eixo, esta potência é transferida ao
47
gerador, que, transforma a potência mecânica em potência elétrica. Ao sair da turbina, o
vapor é resfriado em um condensador e convertido outra vez em água, retornando a
caldeira recuperadora, e dando início a um novo ciclo.
Figura 4.1 - Ciclo Combinado
O ciclo combinado corresponde ao acoplamento do ciclo Rankine ao ciclo
Brayton. Esta combinação prioriza a eficiência de conversão da energia do combustível
para a energia elétrica.
Desta forma as termelétricas a gás natural de ciclo combinado vêm sendo
adotadas em todo o mundo, desde a década de oitenta [66], e vem sendo a solução
escolhida para as termelétricas brasileiras a gás natural nos últimos anos, como por
exemplo: as duas térmicas vencedoras do LEN A-3 de 2011, UTE Baixada Fluminense
[67] e UTE Maranhão III [68]; e a UTE Mauá 3 [69], todas com previsão para entrar em
operação em 2014/1015.
48
Com isso, o tipo de planta termelétrica escolhida para ser otimizada neste
trabalho foi aquela que opera em ciclo combinado. Uma UTE foi modelada no
simulador termodinâmico Thermoflow com as características das últimas térmicas
vencedoras do LEN de 2011, sendo chamada de UTE-Base. A UTE-Base, como indica
o nome, serviu de base durante o processo iterativo de otimização, em busca do melhor
projeto, segundo as restrições impostas pelo Sistema Elétrico Brasileiro. Na próxima
seção será detalhada a construção da UTE-Base.
4.2
UTE-BASE
Como mencionado nas seções anteriores, para o cálculo e otimização da função
objetivo não foram desenvolvidos modelos físicos e termodinâmicos de forma explícita
para o ciclo de geração. Neste trabalho o simulador de processos Thermoflow foi
utilizado com o propósito de resolver os balanços de massa e de energia para que o VPL
seja então calculado.
A UTE-Base foi modelada de acordo com as principais características das
últimas térmicas que ganharam o LEN de 2011: UTE Baixada Fluminense, com
capacidade instalada de 530 MW e garantia física de 430,2 MW [67] ; e UTE Maranhão
III, com capacidade instalada de 499,22 MW e garantia física de 470,7 MW [68].
Ambas as térmicas são de ciclo combinado do tipo 2-2-1, ou seja, faz uso dois
turbogeradores a gás (TGG, conjunto turbina a gás e gerador dedicado), duas caldeiras
recuperadoras de calor (HRSG) e um turbogerador a vapor (TGV, conjunto turbina a
vapor e gerador dedicado).
Nas Figura 4.2 e Figura 4.3 seguem o diagrama esquemático e o infográfico
simplificado do tipo da térmica adotada para UTE-Base.
49
Figura 4.2 - Diagrama da UTE-Base [70]
Figura 4.3 – Infográfico da UTE-Base [70]
4.2.1
UTE-Base no Thermoflow
50
Dentre os programas que o pacote Thermoflow oferece, foi utilizado o GTpro
para simular a UTE-Base.
O GTpro se destina principalmente à obtenção do balanço de massa e energia de
plantas termelétricas em estado estacionário e seus respectivos resultados econômicos a
partir dos custos de equipamentos e de implementação do projeto. Para tanto, o GTpro
dispõe de um vasto banco de dados obtido diretamente de fabricantes de equipamentos,
o qual é renovado periodicamente, tanto no aspecto dos preços praticados no mercado,
quanto no tocante às características técnicas, como performance de máquinas, perdas de
carga, melhorias em modelos existentes e lançamentos. Por todos esses aspectos, os
programas Thermoflow são mundialmente utilizados por empresas em projetos de
usinas termelétricas, plantas de cogeração e casas de força de refinarias [70].
O projeto da UTE-Base seguiu um projeto standard default que o GTpro
oferece. Entretanto, para que ficasse de acordo com a realidade praticada na Baixada
Fluminense
[67] e UTE Maranhão III [68], alguns detalhes de projeto foram
modificados. Assim, a UTE-Base foi projetada com as seguintes principais
características:
1- Geração em ciclo combinado do tipo 2-2-1: dois TGG; duas HRSG; e um TGV;
2- Dois TGG da General Eletric, modelo 7FA (modelo utilizado nas térmicas
vencedoras do LEN 2011);
3- Duas HRSG gerando vapor em três níveis de pressão, projetadas segundo a
menor temperatura da chaminé, definida em 90°C;
4- Turbina a vapor com três níveis de pressão, com reaquecimento e indução de
vapor de baixa pressão;
5- Condensador com desaerador integrado;
6- Torre de Resfriamento com resfriamento a ar.
51
Ao simular este projeto no Thermoflow, segundo máxima eficiência térmica
(correspondendo ao máximo aproveitamento do calor contido nos gases exaustos do
ciclo Brayton pelo ciclo Rankine), foram obtidos os seguintes resultados:
7- Potência Instalada Bruta, correspondente aos dois TGG e a Turbina a Vapor:
522 MW;
8- Potência Instalada Líquida, correspondente aos dois TGG e a Turbina a Vapor:
(descontado os consumos internos e respectivas perdas): 506 MW;
9- Eficiência Bruta do ciclo térmico: 56,82%;
10- Eficiência Líquida do ciclo térmico: 55,08%;
11- Heat Rate Líquido: 6,1991 106Btu/MWh;
12- Pressão, temperatura e vazão mássica do vapor na entrada da turbina a vapor:
124 bar, 566 °C e 362,4 t/h, respectivamente;
13- Pressão, temperatura e vazão mássica do vapor no reaquecimento: 27,6 bar, 566
°C e 433,3 t/h, respectivamente;
14- Pressão, temperatura e vazão mássica do vapor na indução: 3,582 bar, 208,6 °C
e 59,78 t/h respectivamente;
15- Pressão, temperatura e vazão mássica do vapor no condensador: 0,1263 bar,
50,47 °C e 503,4 t/h;
16- Vazão mássica da água de reposição: 4,236 t/h;
17- Vazão mássica de gás natural (combustível): 66,06 t/h;
18- Temperatura dos exaustos na entrada de cada HRSG: 615,9 °C;
19- Custo do Investimento: US$ 355.650.528,00.
A Figura 4.4, a seguir, mostra o resultado da UTE-Base simulado no
Thermoflow.
52
Figura 4.4 - Fluxograma da UTE-Base
53
Algumas das características listadas anteriormente (por exemplo, itens 12, 13 e
14) são variáveis do processo de geração que podem ser prescritas pelo usuário do
Thermoflow. Como descrito anteriormente, a ideia foi variar os seus valores (dentro de
um domínio definido) e, para cada conjunto de novos valores, para que fosse projetada
uma nova térmica com o software Thermoflow. Durante o processo de busca pelo
projeto ótimo, a iteração entre o Thermoflow e o algoritmo de otimização proposto (que
será visto nos próximos capítulos) ocorreu da seguinte forma:
1- O algoritmo fornece os novos valores para as variáveis do processo de geração
escolhidas através dos inputs do Thermoflow;
2- O algoritmo aciona o simulador, para que este projete a nova térmica, segundo
os novos valores das variáveis do processo de geração;
3- O algoritmo faz a leitura do novo projeto simulado através dos outputs do
Thermoflow.
Para que a primeira e a terceira etapa ocorra, é necessário definir as variáveis do
processo de geração como inputs e as variáveis necessárias para o cálculo da função
objetivo como output na UTE-Base, respectivamente.
Assim, na seção a seguir será abordado como foi criado o arquivo da UTE-Base
com seus respectivos inputs e outputs escolhidos pelo usuário, e como foi obtido a
comunicação entre o Thermoflow e Matlab.
4.3
INTERFACE THERMOFLOW E MATLAB
A integração entre o simulador (Thermoflow) e o programa de otimização
(Matlab) foi estabelecida, pois ambos possuem tecnologia COM (Component Object
Model) da Microsoft, que permite a comunicação entre estes dois softwares, [5].
O Thermoflow possui um utilitário chamado U-Link. O U-Link permite o
54
usuário criar uma interface com modelos de plantas previamente construídas a partir dos
programas do Thermoflow (neste caso, o GTpro) e manipulá-las em outro ambiente de
programação (o software Matlab, no caso desta dissertação).
Com isso, após criar a UTE-Base no Thermoflow (utilizando o pacote GTpro)
foi necessário converter este arquivo ao ambiente U-link, para que então fosse
manipulado pelo Matlab.
4.3.1
Convertendo a UTE-Base Em U-Link
Esta seção mostrará o passo a passo de como converter a UTE-Base criada no
ambiente de programação GTpro ao ambiente U-Link. Esta conversão pode ser
realizada em qualquer plataforma de cálculo, porém, neste trabalho foi utilizado o
Matlab. A seguir, na Tabela 4.1 seguem a linhas de programação utilizadas.
Tabela 4.1 - Convertendo o arquivo GTpro em U-Link
Linhas de Comando
TF = actxserver('TFULINKLxx.ULINK')
TF.ModelFileType = 0
TF.ModelFileLoad('Nome_do_arquivo.gtp')
TF.SelectInputVariables
Descrição
Estabelece a comunicação entre Matlab e
Thermoflow.
Obs.: "xx" corresponde a versão do
Thermoflow. Neste trabalho foi utilizada a
versão 23.
Define que será utilizado o pacote de programas
do Thermoflow GTpro.
Carrega o arquivo onde a planta térmica foi
modelada no GTpro.
Seleciona as variáveis de entrada, que serão as
variáveis de otimização.
TF.SelectOutputVariables
Seleciona as variáveis de saída.
TF.FileSave('Nome_do_arquivo.ulk')
Salva a planta térmica no formato U-Link.
Ao aplicar as linhas de programação indicadas na Tabela 4.1 ao arquivo da
55
UTE-Base, foi obtido o projeto base na extensão “ulk”. Assim, ficou possível manipular
as variáveis de entrada, executar um novo projeto e ler as variáveis de saída facilmente.
4.3.2
Executando a UTE-Base no Matlab
A manipulação do projeto da UTE-Base, agora no ambiente U-Link, através da
plataforma de cálculo Matlab foi estabelecida utilizando as linhas de programação
definidas na Tabela 4.2, a seguir:
Tabela 4.2 - Integração entre MATLAB e Thermoflow via U-Link
Linhas de Comando
TF = actxserver('TFULINKLxx.ULINK')
set(TF,'InputValue',i , x i)
Descrição
Estabelece a comunicação entre Matlab e
Thermoflow.
Obs.: "xx" corresponde a versão do
Thermoflow. Neste trabalho foi utilizada a
versão 23.
Fornece o valor da variável de otimização x i à
variável "i" definida como input .
TF.Compute
Executa o simulador e calcula um novo projeto.
yi = TF.OutputValue(j)
Faz a leitura da variável de sáida "j".
56
5
5.1
FUNÇÃO OBJETIVO
VPL DE UMA UTE FLEXÍVEL
O Valor Presente Líquido (VPL) é um dos indicadores financeiros mais
utilizados, por espelhar em termos presentes, a real lucratividade que um determinado
aporte financeiro retornará ao investidor. Para o seu cálculo utilizam-se métodos de
fluxos financeiros descontados, ou seja, todos os fluxos de caixa1 esperados são trazidos
ao valor presente a uma dada taxa de atratividade, que corresponde ao custo de
oportunidade do investidor. Se o cálculo do VPL chegar a um valor positivo, o
investimento tem um retorno superior à taxa de atratividade desejada, tão melhor quanto
maior for o valor do VPL.
De uma maneira bem simples, o VPL de qualquer projeto pode ser calculado da
seguinte forma:
n
VPL = ∑
t a =0
FCta
(1 + TMA)t
a
−I
(5.1).
onde,
ta : período em base anual;
FC ta : fluxo de caixa no período ta , em R$;
I: investimento inicial desprendido, em R$;
TMA: taxa mínima de atratividade (ou taxa de desconto), em %;
n; período final, normalmente correspondendo à vida útil do projeto.
No atual modelo do SEB, os empreendimentos de geração podem reservar um
percentual de sua energia assegurada (garantia física) para comercialização no mercado
cativo (ACR), podendo o restante ser negociado no mercado livre (ACL).
1
Fluxos de caixa são as entradas e as saídas de dinheiro ao longo do tempo que um determinado
investimento promove.
57
Desta forma, nas seções a seguir será abordada a metodologia de cálculo adotada
para a composição de cada variável da equação (5.1), onde os fluxos de caixa foram
compostos das receitas oriundas da negociação tanto no ACR, quanto no ACL. Também
serão identificadas as variáveis de saída do simulador de processo (outputs do
Thermoflow) necessárias para quantificar o VPL do empreendimento termelétrico em
sua totalidade.
Entretanto, para que todas as receitas fossem quantificadas de acordo com cada
projeto termelétrico simulado no Thermoflow, algumas variáveis precisaram ser
dimensionadas, como: a garantia física; a disponibilidade, equação (3.5); a
inflexibilidade, equação (3.6); o custo variável unitário (CVU), equações (3.7) e (3.8); o
custo marginal de operação (CMO); e o preço de liquidação das diferenças (PLD) ou
preço spot da energia.
• Garantia Física da UTE
Como descrito na seção 3.1.3, a garantia física foi definida de acordo com a
seguinte equação:
GF = GF% ⋅ Pot
(5.2).
onde,
Pot: potência instalada da térmica, 1° output do Thermoflow, em MW;
GF% : proporção da potência, em %.
• Disponibilidade da UTE
De acordo com a equação (3.5):
Disp = Pot líq ⋅ FCap max ⋅ (1 − TEIF ) ⋅ (1 − IP )
(5.3).
FCapmax corresponde ao Fator de Capacidade máximo. Este número está
58
relacionado quanto à térmica consegue gerar continuamente, ou seja, este valor leva em
consideração as variações das condições ambientais no decorrer de um ano. Como a
térmica foi dimensionada para uma condição ambiental (temperatura e umidade) média,
então foi assumindo um FCapmax igual a um. Com isso:
Disp = Pot líq ⋅ (1 − TEIF ) ⋅ (1 − IP )
(5.4).
onde,
Potlíq : potência líquida instalada da térmica, 1° output do Thermoflow, em MW;
TEIF: taxa média de indisponibilidade forçada;
IP: taxa de indisponibilidade programada.
• Inflexibilidade da UTE
Inflex = Inflex% ⋅ Disp
(5.5).
onde,
Inflex% : proporção da potência que é inflexível, em %, definida pelo investidor;
Disp: Disponibilidade da usina, MW, equação (5.4).
• Custo Variável Unitário - CVU
De acordo com as equações (3.7) e (3.8):
CVU = i ⋅ e0 ⋅ PREF + C OeM
(5.6).
onde,
COeM : demais custos variáveis da usina, em R$/MWh.
i: fator de conversão, em 106Btu/MWh;
e0 : média anual da taxa de câmbio divulgada pelo Banco Central, em R$/US$;
PREF : preço de referência, em US$/106Btu.
59
A única variável que faz com que o CVU se diferencie para cada projeto é o
fator de conversão i, pois este carrega uma relação com o heat rate (taxa de consumo
energético em combustível por potência elétrica gerada) da planta. Neste trabalho foi
adotado que o fator de conversão fosse igual ao valor do heat rate da térmica,
acrescentado de duas correções: uma em relação ao próprio valor do heat rate; e outra
em relação ao preço do combustível de referência. Estas duas correções serão abordadas
com maiores detalhes quando for apresentado o estudo de caso, capítulo 7.
Desta forma, o fator de conversão adotado foi definido da seguinte forma:
i = HR ⋅ CorHR ⋅ CorCOMB
(5.7).
onde,
HR: heat rate da planta de geração, 2° output do Thermoflow, em 106Btu/MWh;
CorHR : correção do heat rate, em %;
CorCOMB : correção do preço do combustível, em %;
Diante da equação (5.7) pode-se perceber que há um valor de CVU único para
cada valor de heat rate, ou seja, tanto o CVU de caráter estrutural, quanto o de caráter
operacional (ver seção 3.1.4.3) de uma UTE são únicos, de acordo com o valor do seu
heat rate.
• Custo Marginal de Operação - CMO
Os valores de CMO utilizados para calcular o VPL de um empreendimento
térmico foram aqueles disponibilizados pela Empresa de Pesquisa Energética [62]. Os
CMO são resultantes da simulação da operação realizada com o uso do modelo
NEWAVE. O NEWAVE é o modelo utilizado pelo ONS e pela CCEE no planejamento
da operação e na determinação do PLD. Este modelo gera 2.000 cenários de CMO por
mês, ao longo de 60 meses (5 anos). Ou seja, o CMO é uma variável aleatória, onde
60
para cada um dos 60 meses há 2.000 valores com distintas probabilidades de ocorrer.
• Preço de Liquidação das Diferenças - PLD
Os PLD correspondem ao preço da energia no mercado de curto prazo. Para este
estudo, os preços spot foram determinados com base nos CMO (disponibilizados pela
EPE), porém limitados ao valor máximo e mínimo, também divulgados pela EPE [62].
Como estes valores são, na maioria dos casos, o próprio CMO, então, o PLD também é
considerado uma variável aleatória.
Cumpre lembrar que as térmicas a gás podem vir a entrar em operação por
outras razões além daquela descrita acima. Para efeito de cálculo do VPL foi apenas
adotado o critério de despacho por ordem de mérito, ou seja, a UTE somente entra em
operação quando o CVU da mesma for menor do que o PLD. Ou seja, não foi levado
em consideração nos custos e nas receitas a entrada em operação da UTE por razões
elétricas (devido a uma necessidade no sistema) ou por segurança energética.
5.1.1
Investimento
Para o cálculo do investimento foi assumido que o custo total da UTE fosse
descontado ao longo dos anos de construção. Desta forma, o cálculo depende do tipo do
LEN: A-3 ou A-5. Por exemplo, se a térmica for participar de um LEN do tipo A-3,
então o custo total do investimento é dividido em três parcelas e cada uma destas é
trazida ao valor presente. Dito isso, o cálculo do investimento a ser desprendido ficou
da seguinte forma:
61
I=
tipo −1
∑
ta =0
CTF ⋅ e2013 + CCX
tipo ⋅ (1 + TMAa )
ta
(5.8).
Onde,
ta : período em base anual;
CTF : custo do investimento, 3° output do Thermoflow, em US$;
C CX : custo de conexão ao SIN, em R$;
e2013 : taxa de câmbio do ano 2013, R$/US$;
tipo: valor igual a 3 ou 5, dependendo tipo do LEN (A-3 ou A-5);
TMAa : taxa mínima de atratividade anual, em %.
5.1.2
Receita do CCEAR
A remuneração da UTE, quando esta comercializa no mercado cativo foi
valorada de acordo com a equação (3.1), reescrita abaixo:
 RF ⋅ E ⋅ ht + (CVU OP − CVU R ) ⋅ GT ⋅ ht , se PLDt ≥ CVU OP
RCCEAR = 
 RF ⋅ E ⋅ ht , se PLDt < CVU OP
(5.9).
Como não foi escopo deste trabalho analisar o comportamento das receitas (tanto
no ACR, quanto no ACL) segundo as flutuações dos preços spot do combustível, então
na equação acima foi considerado que o valor do CVU OP (declarado no sistema AEGE,
quando da participação no LEN) fosse equivalente ao real custo operacional variável (
CVU R ) da térmica. Ou seja, qualquer ganho ou perda em relação ao preço do
combustível não foi considerado. Com isso, a equação (5.9), tomou o seguinte formato:
RCCEAR = RF ⋅ E ⋅ ht
onde,
RF : receita fixa do contrato, em R$/MWh;
(5.10).
62
E: montante contratado no leilão, em MW;
ht : número de horas no período t.
Substituindo na equação (5.10) o valor da receita fixa anual (convertida para
R$/MWh) calculada na equação (3.12), obtém-se
 x ⋅ (8760 ⋅ GF ⋅ ICB − COP − CEC )
RCCEAR = 
 ⋅ E ⋅ ht
8760 ⋅ x ⋅ GF

(5.11).
Esta receita corresponde ao “aluguel” da térmica ao longo da vida útil do
projeto. Para obter a receita do CCEAR em base anual, basta substituir as horas do
período (ht ) por 8.760 (que corresponde ao número de horas no ano) e o montante a ser
contratada no leilão (E) pela a fração da garantia física destinada ao mercado cativo
(x.GF). Desta forma:
RCCEAR ,ta = x ⋅ (8760 ⋅ GF ⋅ ICB − COP − CEC )
(5.12).
Onde,
ta : período em base anual;
x: fração da GF destinada ao ACR;
8760: número de horas no ano;
GF: garantia física, em MW médios;
ICB: índice de custo benefício, em R$/MWh;
COP: custo de operação, em R$/ano;
CEC: custo econômico de curto prazo, em R$/ano.
Com isso, fazendo uso dos respectivos valores de garantia física, de
disponibilidade, de inflexibilidade, CVU de caráter estrutural e CMO (todos calculado
na seção 5.1), cada uma das parcelas da equação (5.12) pode ser calculada conforme
descrito na seção 3.1.4.6.
63
5.1.3
Receita do Contrato PPA
A remuneração da UTE, quando esta comercializa no mercado livre foi valorada
de acordo com a equação (3.13), reescrita abaixo:
R PPA
 PPPA ⋅ E PPA ⋅ ht + PLDt ⋅ (G PPA − E PPA ) ⋅ ht −

=
G PPA ⋅ CVU R ⋅ ht , se PLDt ≥ CVU OP
 P ⋅ E ⋅ h − E ⋅ PLD ⋅ h , se PLD < CVU
PPA
t
t
t
OP
 PPA PPA t
(5.13).
Onde,
PPPA : preço da energia contratada no PPA, em R$/MWh;
E PPA : energia contratada no PPA, em MW;
GPPA : geração total no PPA da UTE, em MW;
PLDt : preço spot de energia no período t, em R$/MWh;
CVU OP : custo variável unitário de caráter operacional1, em R$/MWh;
CVU R : custo variável unitário real da UTE, em R$/MWh;
ht : número de horas no período t.
A seguir será abordada a forma de cálculo de cada parcela da equação (5.13).
• Energia contratada, EPPA
A energia a ser contratada no ACL correspondeu à fração da garantia física que
não foi destinada ao ACR, sendo assim:
E PPA = (1 − x ) ⋅ GF
(5.14).
• Geração total, GPPA
Analogamente, a energia total gerada no ACL correspondeu à fração da
1
Corresponde ao Custo Variável Unitário declarado a EPE, o qual é utilizado pelo ONS para definir
quando do despacho da térmica. Maiores informações, ver seção 3.1.4.3 desta dissertação.
64
disponibilidade total da térmica, conforme equação a seguir:
GPPA = (1 − x ) ⋅ Disp
(5.15).
• Número de horas no período, ht
Como os valores de PLD são disponibilizados ao mês, então, o número de horas
no período correspondeu ao número de horas no mês.
Conforme descrito na seção 5.1.2, não foi escopo deste trabalho analisar o
comportamento das receitas (tanto no ACR, quanto no ACL) segundo as flutuações dos
preços spot do combustível, então na equação (5.13) foi considerado que o valor do
CVU OP fosse equivalente ao real custo operacional variável ( CVU R ) da térmica.
Diante de tudo que foi exposto, a equação (3.13), em R$/mês, ficou da seguinte
forma:
RPPA,t = (PPPA − PLDt ) ⋅ E PPA ⋅ ht + max(PLD t −CVU OP ; 0) ⋅ GPPA ⋅ ht
(5.16).
Sendo,
GPPA = (1 − x ) ⋅ Disp
(5.17).
E PPA = (1 − x ) ⋅ GF
(5.18).
Onde,
Disp: disponibilidade da UTE, em MW, equação (5.4);
GF: garantia física da UTE, em MW, equação (5.2);
PPPA : preço da energia contratada no PPA, em R$/MWh;
E PPA : energia contratada no PPA, em MW;
GPPA : geração total no PPA da UTE, em MW;
PLDt : preço spot de energia no período t, em R$/MWh;
65
CVU OP : custo variável unitário de caráter operacional, em R$/MWh;
ht : número de horas no período t.
Para converter a equação (5.16) em R$/ano, aplicou-se a metodologia do valor
futuro num período de 12 meses, ou seja:
12
RPPA ,ta = ∑ RPPA ,t ⋅(1 + TMAm )
t
t =1
(5.19).
onde,
t: é o referido mês;
ta : período em base anual;
RPPA,t : remuneração mensal do contrato PPA, equação (5.16);
TMAm : taxa mínima de atratividade mensal, em %.
Como neste trabalho o valor do PLD foi determinado com base no valor do
CMO, então, para cada mês “t” houve várias possibilidades de PLD, ou seja, houve
várias possibilidades de receita mensal do contrato PPA. Para cada mês existiram
cenários de PLD com valores altos e muitos outros com valores baixos. A existência de
vários cenários é o reflexo das incertezas dos afluentes, que foram transmitidas ao
CMO. Quando os volumes estão baixos cresce a probabilidade das térmicas entrarem
em operação e nestes casos os preços spot podem atingir valores muito altos. Diante
disso, houve a necessidade de ser analisado o risco inerente a estas incertezas. Nas
próximas seções será demostrado a metodologia adotada para lidar com estes vários
cenários de PLD.
5.1.4
Outros Custos
Para formalizar a receita líquida de uma térmica há necessidade de considerar
66
outros custos anuais, como por exemplo: o custo do uso do sistema de distribuição e
transmissão; o custo fixo de operação e manutenção; e o custo referente à garantia da
planta. O cálculo dos outros custos anuais ficou:
Coutros ,ta = C fixo + Ctarifa + C garantia
(5.20).
C fixo = Potbruta ⋅ CO&M , fixo ⋅ 8760
(5.21).
Ctarifa = Pot bruta ⋅ Ctusd / tust ⋅ 8760
(5.22).
C garatia = 1% ⋅ (CTF ⋅ e2013 + CCX )
(5.23).
sendo,
onde,
ta : período em base anual;
Potbruta : potência bruta instalada da térmica, 4° output do Thermoflow, em MW;
CO&M , fixo : custo fixo de operação e manutenção, em R$/MWh;
Ctusd / tust : tarifa do uso dos sistemas de distribuição/transmissão, em R$/MWh;
CTF : custo do investimento, 3° output do Thermoflow, em US$;
C CX : custo de conexão ao SIN, em R$;
8.760: número de horas durante um ano.
5.1.5
VPL da UTE
O Valor Presente Líquido da planta de geração pode ser calculado a partir da
definição do fluxo de caixa ( FC ta ) de cada período, conforme apresentado na equação
(5.1).
Para o cálculo do fluxo de caixa de cada período, este trabalho levou em
consideração os descontos referentes à aplicação dos impostos aplicados diretamente a
67
receita bruta e do imposto de renda. Com isso, o fluxo de caixa ( FC ta ) foi definido da
seguinte forma:
[(
]
)
FCta = RCCEAR ,ta + R PPA ,t a ⋅ (1 − ID ) − Coutros ,ta − Deta ⋅ (1 − IR ) + Deta
(5.24).
Onde,
ta : período em base anual;
RCCEAR ,ta : receita no ACR (mercado cativo), em R$/ano, equação (5.12);
RPPA ,ta : receito no ACL (mercado livre), em R$/ano, equação (5.19);
C outros ,ta : outros custos (O&M, tarifas, garantia), em R$/ano, equação (5.20);
ID: Impostos Diretos, em %;
IR: Imposto de Renda, em %;
Deta : depreciação, em R$/ano.
Em posse dos fluxos de caixa, equação (5.24), o Valor Presente foi calculado da
seguinte forma:
T
VP = ∑
[(R
CCEAR ,t a
)
]
+ RPPA ,t a ⋅ (1 − ID ) − C outros ,ta − Deta .(1 − IR ) + Deta
t a =1
(1 + TMAa )
T
(5.25).
Onde,
ta : período em base anual;
RCCEAR ,ta : receita no ACR (mercado cativo), em R$/ano, equação (5.12);
RPPA ,ta : receita no ACL (mercado livre), em R$/ano, equação (5.19);
C outros ,ta : outros custos (O&M, tarifas, garantia), em R$/ano, equação (5.20);
ID: Impostos Diretos, em %;
IR: Imposto de Renda, em %;
Deta : depreciação, em R$/ano;
68
TMAa : taxa mínima de atratividade anual, em %;
T: vida útil do projeto, período total.
A função da equação (5.25) foi trazer os fluxos de caixa para o marco zero, que
corresponde ao início de operação da térmica. Entretanto, para analisar o investimento
corretamente foi necessário trazer este montante ao marco atual, ou seja, antes da
construção ser realizada. Desta forma o Valor Presente Líquido foi calculado da
seguinte forma:
VPL =
VP
−I
(1 + TMAa )tipo
(5.26).
Onde,
tipo: valor igual a 3 ou 5, dependendo tipo do LEN (A-3 ou A-5);
VP: Valor Presente dos fluxos de caixa, equação (5.25), em R$;
I: custo do investimento, equação (5.8), em R$;
TMAa : taxa mínima de atratividade anual, em %;
T: vida útil do projeto, período total.
5.2
ANÁLISE DE RISCO
No Sistema Elétrico Brasileiro as afluências naturais são aleatórias em razão da
incidência de chuvas. O uso de séries (ou cenários) históricas (dada pelo registro de
afluências observadas no passado) é insuficiente para compor uma amostra de tamanho
necessário para estimar índices de risco com incertezas aceitáveis. Entretanto, as
características básicas da série histórica podem ser capturadas por modelos estocásticos
capazes de produzir séries sintéticas de afluências, diferentes da série histórica, mas
igualmente prováveis. Como já mencionado nas seções anteriores, o modelo NEWAVE
realiza este tipo de procedimento, sendo este utilizado para o planejamento da operação
69
hidrotérmica no Brasil [3].
O NEWAVE gera, para cada mês, 2.000 cenários de custo marginal de operação
(representando o custo de operação do sistema hidrotérmico). Os CMO apresentam uma
volatilidade decorrente da aleatoriedade das vazões afluentes. Como o preço da energia
no mercado de curto prazo (o preço spot ou o PLD) é calculado pelo CMO, essa
característica influencia os resultados econômicos da geração e comercialização de
energia elétrica no mercado livre.
Com isso, a proposta da função objetivo neste trabalho foi capturar o risco de um
empreendimento termelétrico flexível através da distribuição de probabilidade de cada
série de cenários de cada um dos 60 meses, sendo este risco medido segundo a
probabilidade de ocorrência de um retorno ser inferior ao desejado.
De acordo com a planilha disponibilizada pela EPE [62], os valores de CMO são
apresentados da seguinte forma:
Tabela 5.1 - Amostra de valores de CMO [62]
mês1
mês2
mês3
...
mês59
mês60
cenário1
148,04
43,17
34,20
...
19,50
19,48
cenário2
88,94
66,30
72,56
...
41,02
51,73
cenário3
95,19
92,23
112,13
...
10,63
1,62
...
...
...
...
...
...
cenário1999
65,13
0
0
...
20,73
37,06
cenário2000
112,53
207,46
210,10
...
19,18
40,68
Para o cálculo do VPL os CMO acima foram utilizados como preço spot, porém
limitados a um valor máximo e mínimo, sendo estes limites também disponibilizados
pela EPE.
Para cada cenário há respectivamente um valor de CMO para cada um dos 60
meses. Com isso, verificou-se que, ao utilizar estas previsões de CMO no cálculo da
70
remuneração da energia negociada no mercado livre, equação (5.19), diferentes cenários
de VPL seriam gerados, sendo alguns vantajosos e outros até negativos.
Diante disso, foi definida como função objetivo de um projeto termelétrico o
valor do VPL, no qual apenas 5% de todos os possíveis VPL (gerados a partir de cada
cenário de CMO) fossem menores do que o próprio. Ou seja, 95% de todos os possíveis
VPL teriam valores maiores do que o VPL escolhido para representar a função objetivo.
Assim, a função objetivo assumiu o valor do VPL com apenas 5% de risco de não
ocorrer.
Desta maneira, os cenários de CMO foram utilizados para determinar as
diferentes trajetórias dos preços spot, o que caracterizou um procedimento de simulação
de Monte Carlo. Cada conjunto de cenários (composto pelo respectivo cenário de cada
mês) foi considerado como uma trajetória dos PLD. Em uma Simulação de Monte
Carlo, cada conjunto de cenários de PLD é visto como uma amostra aleatória de todos
os possíveis conjuntos de cenários.
5.3
SIMULAÇÃO DE MONTE CARLO
O método de Monte Carlo é um procedimento numérico que se utiliza de
números aleatórios, ou pseudoaleatórios, para computar algumas quantidades não
necessariamente aleatórias, com base na Lei dos Grandes Números e no Teorema do
Limite Central. Devido à simplicidade das ideias envolvidas na concepção do método e
ao grande avanço dos computadores pessoais, o método de Monte Carlo se tornou uma
poderosa e atrativa ferramenta para lidar com problemas da teoria financeira [72].
A técnica de Monte Carlo consiste em gerar valores aleatórios para cada
distribuição de probabilidades dentro de um modelo com o objetivo de produzir
centenas ou milhares de cenários. A distribuição dos valores calculados (para cada caso)
71
deve refletir a probabilidade de ocorrência dos mesmos.
O método da Simulação de Monte Carlo (SMC) consiste basicamente nas
seguintes etapas [72, 73]:
1- Selecionar as variáveis aleatórias em estudo;
2- Identificar a função de probabilidade (ou distribuição) acumulada de cada
variável aleatória selecionada. Esta função pode ser uma distribuição teórica
(por exemplo, Uniforme, Triangular, Normal, Beta, Weibull, etc.) ou uma
distribuição empírica qualquer. A função de distribuição acumulada F(x) é a
qual fornece a probabilidade, P, em que X é menor ou igual a x. Assim:
F ( x ) = P( X ≤ x ), onde 0 ≤ F ( x ) ≤ 1
(5.27).
3- Obter a função inversa G(F(x)), que calcule o valor de x para o valor de F(x)
correspondente. Então:
G (F (x )) = x
(5.28).
4- Validar o modelo estatístico pressuposto (passos 2 e 3), pela utilização de um
teste de aderência. Para a elaboração do teste de aderência, utiliza-se, por
exemplo, os testes Qui-Quadrado, Kolmogorov-Smirnov, ou o teste de
Anderson-Darling;
5- Gerar um número aleatório, r, a partir de uma distribuição uniforme com
intervalo [0,1]. Onde o valor de r representará o valor da função de distribuição
acumulada, F(x).
6- Gerar o valor para a variável de risco em estudo através da função inversa, a
partir do número aleatório r. Ou seja:
G(r ) = x
(5.29).
O número r é gerado de uma distribuição uniforme [0,1] para permitir igual
72
oportunidade de que o valor x seja gerado em cada percentil do domínio com as
probabilidades equivalentes.
Apesar da técnica de Monte Carlo satisfazer o objetivo de se obter uma
amostragem aleatória, muitas partes da distribuição podem vir a ser subestimadas ou
superestimadas. Ainda assim, a distribuição pode ser melhor replicada se o número de
interações for muito grande. Ou seja, repetir os passos 5 e 6 em centenas ou milhares de
vezes [73].
5.4
EXPRESSÃO FINAL PARA A FUNÇÃO OBJETIVO
A princípio talvez não fosse necessário utilizar a Simulação de Monte Carlo para
gerar os cenários de PLD, uma vez que, conforme mencionado anteriormente, a EPE
[62] disponibiliza 2.000 cenários de CMO para cada um dos 60 meses. Entretanto há
dois grandes motivos para que fosse usado a SMC: primeiro porque a EPE apenas
disponibiliza cenários de CMO para os cinco primeiros anos de operação, e neste
trabalho, foi analisado um horizonte de 15 anos; segundo, para ter uma estimativa mais
precisa da função objetivo proposta, é considerado como boa prática um grande número
de simulações, por exemplo, 5.000, 10.000, ou maior [48, 72, 73].
Para contornar o primeiro problema citado, foi adotado o mesmo método visto
nas referências [1, 2], onde foi assumido que as funções de probabilidade acumuladas
dos cenários dos últimos 120 meses (aqueles não fornecidos pela EPE), teriam o mesmo
perfil de distribuição de probabilidade dos últimos 12 meses daqueles disponibilizados
pela EPE, respectivamente. Em relação ao segundo problema, foi aplicado a SMC com
10.000 iterações.
Desta forma, para determinar os cenários do PLD, foram criados 10.000 cenários
de CMO para cada um dos 180 meses (15 anos), onde cada um destes meses
73
correspondeu a uma variável aleatória no método SMC. Para simular o conjunto de
10.000 cenários de CMO, que representam o conjunto total de possibilidades, foram
aplicados os passos descritos na seção anterior da seguinte forma:
1- Selecionar as variáveis aleatórias: as variáveis aleatórias em estudo foram os
preços spot de cada um dos 180 meses;
2- Identificar a função de probabilidade acumulada de cada variável aleatória: neste
trabalho foi utilizada uma função empírica.
3- Obter a função inversa de F(x): matematicamente a inversa generalizada de uma
função de distribuição G(F(x)) é chamada função quantil da função de
distribuição F(x) [74]. Assim, a simulação de um novo cenário de uma variável
aleatória foi obtida utilizando uma função interna do MATLAB chamada
quantile.m. Esta função interpola um conjunto de cenários originais e gera,
empiricamente, um novo cenário, a partir do valor da F(x) fornecido. Ou seja:
x = quantile( X , F ( x )) = G(F ( x ))
(5.30).
4- Validar o modelo estatístico pressuposto: para validar o modelo admitido no
passo 3, foi utilizado o teste de hipótese Kolmogorov-Smirnov (KS) para duas
amostras. Este tipo de teste é utilizado para determinar se duas distribuições de
probabilidade subjacentes diferem uma da outra. Neste teste, a hipótese nula diz
que as duas amostras são oriundas de uma mesma distribuição de probabilidade
e a hipótese alternativa diz o contrário. O teste KS neste trabalho foi realizado
com uso de outra função interna do MATLAB: kstest2.m. Na prática, o teste foi
feito utilizando a equação (5.30), por exemplo, 10.000 vezes. Em seguida, foram
comparados esses 10.000 valores com os 2.000 valores originais de CMO,
conforme equação abaixo.
h = kstest 2( X , X 10.000 )
(5.31).
74
Nesta dissertação obteve-se h igual à zero (ou seja, hipótese nula) para as 180
variáveis aleatórias consideradas. Desta forma, pode ser assumido que os valores
presentes em X (que representam os 2.000 valores de CMO da variável aleatória
analisada) e os valores gerados ( X 10.000 ) são oriundos da mesma distribuição de
probabilidade acumulada.
5- Gerar um número aleatório, r, a partir de uma distribuição uniforme com
intervalo [0,1]: Neste caso foi utilizado a função rand() do MATLAB.
6- Gerar o valor para a variável de risco em estudo através da função inversa: após
o número aleatório r ter sido gerado no passo 5, o mesmo foi substituído na
equação (5.30), assim:
cmoi ,c = quantile(CMOi , ri ,c )
(5.32).
Onde,
i: corresponde ao mês analisado, i = 1 a 180;
c: cenário criado, sendo neste caso, c = 1 a 10.000;
CMOi : conjunto dos 2.000 valores disponibilizado pela EPE no mês i.
Com isso, fazendo uso da equação (5.32) foi possível gerar 10.000 cenários para
cada um dos 180 meses, calculando assim, 10.000 possíveis VPL. A partir destes
10.000 prováveis valores de VPL foi possível definir a função objetivo da seguinte
forma:
VPL5% = quantile( VPL10.000 , 5%)
(5.33).
Onde,
VPL10.000 : conjunto dos 10.000 possíveis VPL calculados (um vetor coluna);
Nesta dissertação, o VPL5% , equação (5.33), foi escolhido como função objetivo
75
de um projeto termelétrico. Assim, durante o processo de otimização, foi assumido um
risco de apenas 5% do valor escolhido para representar a função objetivo não ocorrer.
76
6
ALGORITMO DE OTIMIZAÇÃO
No capítulo anterior foi descrita e formalizada a função objetivo a ser otimizada.
Nas seções a seguir deste capítulo será discutido a metodologia adotada como algoritmo
de otimização, sendo baseado no Método da Evolução Diferenciada fazendo uso de
funções interpoladas pelo Método Superfície de Resposta.
6.1
MÉTODO DA EVOLUÇÃO DIFERENCIADA
No algoritmo da Evolução Diferenciada tradicional os parâmetros de mutação
(F), de cruzamento (CR) e a estratégia de mutação são definidos no início do método e
permanecem imutáveis durante todo o processo de otimização.
Como mencionado na seção 2.2.1, o método adaptativo permite que estes
parâmetros (F e CR) e a estratégia se modifiquem a cada iteração, para cada um dos
indivíduos da população, melhorando o desempenho da otimização.
Nas seções, a seguir, será detalhado o algoritmo tradicional e apresentado o
algoritmo da ED utilizado neste trabalho (o ED_mod), levando em consideração as
adaptações dos fatores de cruzamento, mutação e da estratégia.
6.1.1
Método Tradicional
Como mencionado na referência bibliográfica (capítulo 2), o algoritmo da ED
tradicional [4, 7] inicia-se a partir da geração aleatória (com distribuição de
probabilidade uniforme) de uma população com NP indivíduos (vetores) dentro do
domínio delimitado, conforme a equação abaixo:
popi , j = LI j + Ri , j ⋅ (LS j − LI j )
Onde,
(6.1).
77
i = 1, 2, 3, ... , NP;
j = 1, 2, 3, ... , Nvar;
NP: tamanho da população;
Nvar: número de variáveis, ou seja, dimensão de cada vetor (indivíduo);
LI j : limite mínimo da variável j;
LS j : limite máximo da variável j;
R j : número aleatório com distribuição uniforme entre 0,1 ;
pop : matriz referente a população inicial, contendo os NP indivíduos.
No método da evolução diferenciada tradicional há basicamente três processos.
Ou seja, para cada vetor da geração atual (G) ocorrem três modificações: mutação;
cruzamento; e seleção.
6.1.1.1 Mutação
No processo de mutação, para cada vetor da população corrente é gerado um
vetor mutante, também denominado de vetor diferencial. Este vetor pode ser criado a
partir de vários tipos de estratégias. A seguir, segue a lista das estratégias mais comuns
encontradas na literatura [7, 32, 36, 38, 75]:
1- DE/rand/1
(
viG = pop rG1 + F ⋅ pop rG2 − pop rG3
)
(6.2).
2- DE/rand/2
(
)
(
viG = pop rG1 + F ⋅ pop rG2 − pop rG3 + F ⋅ pop rG4 − pop rG5
3- DE/best/1
)
(6.3).
78
(
G
viG = pop otimo
+ F ⋅ pop rG1 − pop rG2
)
(6.4).
4- DE/best/2
)
(
)
(6.5).
(
)
(
)
(6.6).
(
)
(
)
(6.7).
)
(
)
(6.8).
(
G
viG = pop otimo
+ F ⋅ pop rG1 − pop rG2 + F ⋅ pop rG3 − pop rG4
5- DE/rand-to-best/1
G
viG = pop rG1 + F ⋅ pop otimo
− pop rG1 + F ⋅ pop rG2 − pop rG3
6- DE/current-to-best/1
G
v iG = pop iG + F ⋅ pop otimo
− pop iG + F ⋅ pop rG1 − pop rG2
7- DE/current-to-rand/1
(
viG = pop iG + rand ⋅ pop rG1 − pop iG + F ⋅ pop rG2 − pop rG3
Onde,
G: geração corrente;
i = 1, 2, 3, ... , NP, posição do vetor corrente;
i ≠ r1 ≠ r2 ≠ r3 ≠ r4 ≠ r5 : índices (ou posições) escolhidos aleatoriamente;
F: fator de mutação;
pop G : matriz referente a população da geração G, contendo os NP indivíduos;
G
popotimo
: corresponde ao melhor indivíduo entre todos os existentes em pop G .
De acordo com a referência [7], o método da ED tradicional utiliza a estratégia
de mutação DE/rand/1, equação (6.2), sendo F o fator de mutação de valor real e
constante no intervalo [0, 2]. Este fator controla a amplificação da variação da diferença
entre os vetores mutantes com a finalidade de controlar uma convergência prematura.
6.1.1.2 Cruzamento
79
Tipicamente no Método da Evolução Diferenciada são utilizados dois tipos de
cruzamento: o binomial e o exponencial [4, 7]. Na última década, o cruzamento
binomial foi ganhando popularidade frente ao exponencial [38]. Desta forma, nesta
dissertação foi adotado o cruzamento do tipo binomial.
No processo de cruzamento do método da ED tradicional, para cada vetor da
população atual é criado um vetor experimental. Este vetor é composto pelas
componentes do vetor da população atual ( popiG ) e pelas componentes do vetor
mutante ( viG ), obtendo o seguinte formato:
{
uiG = uiG,1 , uiG, 2 ,⋅ ⋅ ⋅, uiG, N var −1 , uiG, N var
}
(6.9).
Cada componente do vetor experimental é escolhida a partir da comparação de
um número aleatório (entre 0 e 1) com o fator de cruzamento, CR, conforme equação a
seguir.
viG, j , se (rand j < CR ) ou (randi = j )
u iG, j = 
G
 popi , j , se (rand j ≥ CR ) e (randi ≠ j )
(6.10).
Onde,
i = 1, 2, 3, ... , NP;
j = 1, 2, 3, ... , Nvar;
rand j : número aleatório entre [0, 1] para cada j;
randi: escolha aleatória de um dos índices {1, 2, ..., Nvar}.
Pode ser observado na equação (6.10) que pelo menos uma componente do vetor
mutante estará presente no vetor experimental garantindo que o vetor experimental seja
sempre diferente do vetor da geração atual.
Uma vez que o vetor experimental esteja estruturado, avalia-se o valor da função
80
objetivo do mesmo. Note que a função objetivo pode ser representada tanto por uma
função implícita, quanto explícita, ou até por uma função interpolada.
6.1.1.3 Seleção
No processo de seleção ocorre a comparação do valor da função objetivo do
vetor atual ( popiG ) como o valor da função objetivo calculado para o vetor
experimental ( uiG ) de acordo com a seguinte equação:
( )
(
 popiG , se f u iG for pior do que f popiG
popiG +1 =  G
u i , caso contrário
)
(6.11).
Onde,
i = 1, 2, 3, ... , NP;
uiG : vetor experimental;
f (uiG ) : valor da função objetivo do vetor experimental;
popiG : vetor i da população atual;
(
)
f popiG : valor da função objetivo do vetor i da população atual;
popiG +1 : vetor i da população da próxima geração.
De acordo com a equação (6.11), se o valor da função objetivo aplicada ao vetor
experimental for igual ou melhor do que o valor dela aplicada ao vetor da população
corrente, então o vetor experimental substitui o vetor da população na geração seguinte.
Caso contrário, o vetor da população atual é mantido na próxima geração.
Por fim, estes procedimentos continuam até que algum critério de parada seja
cumprido.
Na Figura 6.1 pode ser visto um pseudo-algoritmo do Método da Evolução
81
Diferenciada tradicional.
01 Início
02
Definir Nvar , o número de variáveis
03
Definir NP , o número da população
04
Definir F e CR ;
05
Criar população inicial aleatória {pop i,0 |i = 1 a NP }
06
Avaliar a função objetivo para cada indivíduo.
07
Iniciar a contagem das gerações: g = 0;
08
While (algum critério de parada não é atingido) do
09
g = g + 1;
10
For i = 1 a NP
11
Escolher aleatoriamente r 1 ≠ r 2 ≠ r 3 ≠ i
12
v i = pop r1,g + F .(pop r2,g - pop r3,g )
13
Gerar j rand = randint(1, Nvar )
14
For j = 1 a Nvar
15
if j = j rand or rand < CR i
16
u i,j,g = v i,j,g
17
else
u i,j,g = pop i,j,g
18
19
end if
20
end for
21
Avaliar a função objetivo de u i,g
22
if f (pop i,g ) melhor f (u i,g )
23
pop i,g+1 = pop i,g
24
else
pop i,g+1 = u i,g
25
26
end if
27
end for
28
end while
29 end
Figura 6.1 - Pseudo-Algoritmo: ED tradicional [7]
6.1.2
Método Adotado: ED_Mod
O método da Evolução Diferenciada adotado neste trabalho foi baseado no
método tradicional visto na seção anterior, porém com três grandes modificações.
As duas primeiras modificações foram em relação à adaptação dos fatores de
cruzamento e mutação, associados a cada estratégia.
O algoritmo adotado faz uso de um pool de estratégias. Em cada iteração, cada
indivíduo da população é associado a uma das estratégias presentes no pool, de acordo
com a probabilidade que cada estratégia tem de ser escolhida.
82
Nas seções a seguir serão apresentados como ocorrem as adaptações dos fatores
de cruzamento, mutação e das estratégias, respectivamente, bem como, o pool de
estratégias.
Por fim, será apresentada a integração destas modificações ao método da ED
tradicional, formalizando assim, a Evolução Diferenciada adotada: ED_mod.
6.1.2.1 Adaptação dos Fatores de Cruzamento e Mutação
As primeiras modificações ao método da ED tradicional estão ligadas as
gerações independentes de novos fatores de mutação (F) e cruzamento (CR) para cada
indivíduo, em cada geração. Neste trabalho foi adotada a metodologia apresentada na
referência [36], que é uma modificação ao que foi utilizado no algoritmo JADE [35].
No algoritmo JADE o fator de cruzamento é gerado segundo uma distribuição
normal, com média igual a µ CR e desvio padrão igual a 0,1. Depois de gerado, o fator é
truncado em [0, 1].
CRi = randni (µ CR , 0,1)
(6.12).
A média µ CR inicia com valor igual a 0,5, entretanto, a cada geração o seu valor
é atualizado da seguinte forma:
µ CR = (1 − c ) ⋅ µ CR + c ⋅ meanA (S CR )
(6.13).
Onde,
NP: número total de indivíduos da população;
c: constante com valor igual a 0,1;
S CR : conjunto composto por todos os fatores de cruzamento que obtiveram
sucesso na geração anterior;
mean A (*) : operador matemático responsável pela média aritmética;
83
De forma a manter a diversidade da população, analogamente, o fator de
mutação é gerado independentemente para cada indivíduo segundo uma distribuição de
Cauchy, com parâmetro de locação e de escala iguais a µ F e 0,1, respectivamente.
Fi = randci (µ F , 0,1)
(6.14).
i = 1, 2, 3, ... , NP;
NP: número total de indivíduos da população;
Caso o valor Fi seja maior do que 1, o mesmo é truncado em 1, ou gerado
novamente se for menor ou igual a zero.
O parâmetro de locação µ F inicia com valor igual a 0,5. Entretanto, a cada
geração, este parâmetro é atualizado da seguinte forma:
µ F = (1 − c ) ⋅ µ F + c ⋅ meanL (S F )
(6.15).
Onde mean L (*) corresponde ao operador matemático responsável pela média de
Lehmer, assim:
SF
mean L (S F ) =
∑F
j =1
2
j
(6.16).
SF
∑F
j =1
j
Onde,
c: constante com valor igual a 0,1;
S F : conjunto composto por todos os fatores de mutação que obtiveram sucesso
na geração anterior.
Neste trabalho foi utilizada esta metodologia de adaptação dos parâmetros,
84
porém as equações (6.13) e (6.15) foram calculadas de acordo com o que foi proposto
na referência [36].
As equações (6.13) e (6.15) mostram que a média µ CR e o parâmetro de locação
µ F são atualizados de acordo com os valores de CRi e Fi que obtiveram sucesso na
geração anterior, respectivamente. Esta forma de atualização pressupõe que todos os
valores de sucesso têm o mesmo peso. No entanto, intuitivamente, se, por exemplo, os
fatores CR1 e F1 são capazes de obter maior recompensa1 do que CR2 e F2 , então seria
mais adequado que fosse atribuído um peso mais elevado a estes, quando comparados
com os CR2 e F2 .
Com base na análise acima, a referência [36] propôs uma abordagem, onde os
valores médios ponderados são utilizados nas equações (6.13) e (6.15). Desta maneira,
os operadores mean A (*) e mean L (*) são substituídos respectivamente por:
mean wA (S CR ) =
S CR
∑w
j =1
j
⋅ CR j
j
⋅ F j2
SF
meanwL (S F ) =
∑w
j =1
SF
∑w
j =1
j
⋅ Fj
(6.17).
(6.18).
Onde,
S CR : conjunto composto por todos os fatores de cruzamento que obtiveram
sucesso na geração anterior;
S F : conjunto composto por todos os fatores de mutação que obtiveram sucesso
na geração anterior;
w j : peso atribuído a cada fator que obteve sucesso.
1
Quanto maior o sucesso na busca pelo o ponto ótimo, maior deve ser a recompensa atribuída aos
parâmetros que foram responsáveis por achar um ponto melhor.
85
Desta forma, os novos operadores meanwA (* ) e meanwL (* ) , equações (6.17) e
(6.18), levam em consideração as respectivas recompensas, de acordo com o tamanho
do sucesso que cada fator proporcionou, sendo estas recompensas quantificadas de
acordo com a metodologia do fitness improvement, proposta na referência [76].

f

ηj =
f

δ
(u ) ⋅ f ( pop ) − f (u ) , se for minimizar
(u ) ⋅ f ( pop ) − f (u ) , se for maximizar
δ
G
j
G
j
G
j
G
j
G
j
G
j
(6.19).
Onde,
j = 1, 2, ..., n, onde n é o número total de CRi e Fi que obtiveram sucesso.
Repare que, n = S CR = S F ;
δ : melhor valor de função objetivo encontrada até o momento;
(
)
f pop Gj : valor da função objetivo do vetor j da população corrente;
( )
f u Gj : valor da função objetivo do vetor experimental.
Na referência [37] a equação (6.19) é utilizada para quantificar as recompensas
de cada par de fatores ( CRi e Fi ). Entretanto, nesta dissertação não foi utilizada a
formulação proposta na equação (6.19) por três grandes motivos:
1- Caso δ assumisse um valor igual a zero durante uma minimização ou f (u Gj )
fosse zero durante um processo de maximização, então, todos os valores η j
também seriam nulos, o que não faria nenhum sentido uma vez que se deseja
mensurar as recompensas, ou seja, η j > 0 ;
2- Caso f (u Gj ) fosse igual a zero durante uma minimização ou δ fosse zero
86
durante um processo de maximização, o valor η j correspondente assumiria um
valor infinito, o que iria trazer problemas numéricos no momento em que os
operadores meanwA (* ) e meanwL (* ) , equações (6.17) e (6.18), fossem aplicados;
3- Caso o sinal de δ fosse diferente de f (u Gj ) , o valor η j correspondente
assumiria um valor negativo, o que também não iria fazer nenhum sentido, pois
haveria uma recompensa negativa.
Em vista disso, antes de utilizar a equação (6.19), os valores das funções
objetivo f (u Gj ) e f ( pop Gj ) , partes integrantes desta equação, foram normalizados
dentro de uma escala linear entre 1 e 100, conforme equação a seguir:
 f − f min
f norm = 
 f max − f min

 ⋅ (100 − 1) + 1

(6.20).
[(
) ( )]
(6.21).
[(
) ( )]
(6.22).
Sendo
f min = min f popiG ; f u iG
f max = max f popiG ; f u iG
Onde
i = 1, 2, 3, ... , NP;
NP: número total de indivíduos da população;
Assim, substituindo a equação (6.20) na equação (6.19), obtêm-se:
 δ norm
⋅ f norm pop Gj − f norm u Gj , se for minimizar
G
f
 norm u j
ηj = 
G
 f norm u j ⋅ f
G
G
, se for maximizar
norm pop j − f norm u j
 δ
norm

( )
( )
(
)
( )
(
)
( )
(6.23).
j = 1, 2, ..., n, onde n é o número total de CRi e Fi que obtiveram sucesso;
87
δ norm : melhor valor de função objetivo encontrada até o momento, normalizada
entre [1, 100];
(
)
f norm pop Gj : valor da função objetivo vetor j da população corrente,
normalizada entre [1, 100];
( )
f u Gj : valor da função objetivo do vetor experimental, normalizada entre [1,
100];
É fácil perceber que ao utilizar a equação (6.23) os três problemas citados
anteriormente são contornados, pois os valores de função objetivo normalizados nunca
serão menores do 1, nunca serão maiores do que 100 e sempre serão positivos.
Com isso, em posse das recompensas quantificadas, os pesos w j utilizados nas
equações (6.17) e (6.18) são calculados da seguinte forma:
wj =
ηj
n
∑η
j =1
j
(6.24).
Onde,
n: número total de CRi e Fi que obtiveram sucesso;
i = 1, 2, 3, ... , NP;
NP: número total de indivíduos da população;
η j : recompensa calculada para cada par CRi e Fi que obtiveram sucesso, a
partir da equação (6.23).
6.1.2.2 Pool de Estratégias de Mutação
No algoritmo da ED tradicional é adotada uma única estratégia de mutação no
88
decorrer de todo o método. Nos últimos 15 anos foram sugeridas outras diferentes
estratégias no intuito de aumentar a robustez do algoritmo subjacente. As estratégias
mais comuns encontradas na literatura foram listadas na seção 6.1.1.1, equações (6.2) a
(6.8).
Apesar destas estratégias aumentarem a robustez do algoritmo, a definição de
qual delas seria mais adequada para otimizar um determinado problema pode se tornar
inviável, por se configurar um processo de tentativa e erro.
Desta forma, no algoritmo utilizado neste trabalho foi adotado um pool de
estratégias contendo três distintas estratégias: a DE/rand-to-pbest/1, a qual foi utilizada
por [34-37]; a DE/rand/2, equação (6.3); e a DE/current-to-rand/1, equação (6.8).
• DE/rand-to-pbest/1
Esta estratégia de mutação vem sendo muito utilizada por alguns autores. Ela é
oriunda da estratégia DE/rand-to-best/1, equação (6.6), entretanto, com uma diferença
no vetor “best”. O vetor “pbest” desta nova estratégia não necessariamente será o
melhor vetor de toda a população. O mesmo pode aleatoriamente ser o segundo, o
terceiro, ou até o quarto melhor. Em cada iteração (ou geração) é definido um grupo dos
melhores vetores. A quantidade de melhores indivíduos presentes neste grupo é definida
pelo usuário, porém na literatura encontram-se valores referentes de 5% a 15% da
quantidade total de indivíduos de toda a população. Neste sentido, nesta dissertação foi
assumido um grupo “pbest” contendo 15% dos melhores vetores, ficando o vetor
mutante gerado por esta estratégia da seguinte forma:
(
)
(
viG = poprG1 + Fi ⋅ pop Gp _ otimo − poprG1 + Fi ⋅ poprG2 − po~prG3
Onde,
i = 1, 2, 3, ... , NP;
)
(6.25).
89
Fi : fator de mutação correspondente a cada i da geração G;
pop Gp _ otimo : corresponde a um vetor escolhido aleatoriamente entre o grupo
100p% dos melhores vetores1 da geração G;
pop rG1 : vetor escolhido aleatoriamente da população pop, sendo r1 ≠ i ;
pop rG2 : vetor escolhido aleatoriamente da população pop, sendo r1 ≠ r2 ≠ i ;
po~prG3 : vetor escolhido aleatoriamente da população pop união com a população
A, sendo r1 ≠ r2 ≠ r3 ≠ i ;
Outra modificação proposta na referência [35] nesta estratégia está ligada a
utilização da população A corresponde ao Archive, que tem como objetivo manter a
diversidade das gerações.
Archive é o grupo de vetores que não obtiveram sucesso na geração anterior.
Inicialmente este grupo começa vazio e vai sendo preenchido com os vetores que seriam
eliminados durante o processo de seleção. A cada geração o grupo A é renovado.
Existem dois limites ao grupo A: não há vetores iguais; e o tamanho da população é
limitado a NP. Se o tamanho da população de A ultrapassar a NP, então, alguns vetores
de A são aleatoriamente removidos, mantendo no máximo NP distintos indivíduos.
• DE/rand/2
Apesar da estratégia DE/rand/1, equação (6.2), ser uma das mais utilizadas [38,
75], a DE/rand/2, equação (6.3), foi preferida a pertencer ao pool, pois esta tem chance
de gerar melhores permutações, o que permite gerar diferentes vetores experimentais.
Isso se deve ao fato da estratégia DE/rand/2 considerar dois diferenciais, ao passo que a
1
O grupo 100p% corresponde ao grupo formado por “p” melhores vetores de toda a população pop.
90
estratégia DE/rand/1 considera apenas um.
• DE/current-to-rand/1
Essa estratégia não faz uso do cruzamento, ou seja, esta gera automaticamente o
vetor experimental, sendo uma combinação do próprio vetor mutante com o vetor da
população correspondente (vetor alvo), ver equação (6.8). Esta estratégia foi muito
utilizada por vários autores que obtiveram resultados bastante competitivos [31, 32, 37].
6.1.2.3 Adaptação das Estratégias
O mecanismo de seleção das estratégias foi baseado no método AdapSS
(Adaptive Strategy Selection) adotado em [37]. O objetivo deste método é selecionar
automaticamente as estratégias de mutação presentes no pool, de acordo com o
desempenho de cada uma durante o processo de otimização.
Para isso, o método faz uso de duas avaliações. A primeira correspondendo a
recompensa de cada estratégia, a qual quantifica o impacto das estratégias durante a
otimização. A segunda referente ao mecanismo de seleção das estratégias que, com base
nas recompensas calculadas, escolhe qual estratégia deve ser aplicada.
No método AdapSS cada estratégia tem uma probabilidade definida de ser
escolhida. E durante cada iteração a probabilidade de cada estratégia é adaptada de
acordo com o seu desempenho nas iterações anteriores.
Desta forma, como neste trabalho há três estratégias presentes no pool, as
probabilidades de cada estratégia são apresentadas no seguinte formato:
3


P(g ) = {p1 ( g ), p 2 (g ), p3 ( g )}, onde  ∀ g : p min ≤ p a ≤ p max ; ∑ p a (g ) = 1
a =1


Onde,
(6.26).
91
g: geração (ou iteração) corrente;
a: estratégia “a”, onde a = 1, 2 e 3;
p1 (g ) : corresponde a probabilidade da primeira estratégia ser escolhida;
p 2 ( g ) : corresponde a probabilidade da segunda estratégia ser escolhida;
p3 ( g ) : corresponde a probabilidade da terceira estratégia ser escolhida;
p min : valor mínimo de uma probabilidade, para evitar que seja zero;
p max = 1 − ( a − 1)⋅ p min : sendo a a quantidade total de estratégias.
A cada nova geração as probabilidades de cada estratégia podem ser
recalculadas de duas formas: a primeira forma é atribuída à estratégia que obteve maior
sucesso na geração anterior, equação (6.27); e a outra forma é atribuída às estratégias
restantes, equação (6.28).
[
]
p a* ( g + 1) = p a* ( g ) + β ⋅ p max − p a* ( g )
(6.27).
∀a ≠ a * : p a ( g + 1) = p a ( g ) + β [ p min − p a ( g )]
(6.28).
Onde,
a * : estratégia que teve maior sucesso na geração G;
a: estratégia “a”, podendo ser igual a 1, 2 ou 3;
β ∈ (0, 1] : parâmetro que controla a taxa de crescimento das probabilidades.
A estratégia a * é definida a partir da quantificação numérica da qualidade de
cada estratégia, sendo que cada qualidade é calculada da seguinte maneira:
qa (g + 1) = qa ( g ) + α ⋅ [ra ( g ) − qa (g )]
Onde,
g: geração (ou iteração) corrente;
(6.29).
92
a: estratégia “a”, onde a = 1, 2 e 3;
α ∈ (0, 1] : parâmetro que controla a taxa de crescimento das qualidades;
ra ( g ) : recompensa referente ao sucesso de cada estratégia durante a geração G.
Baseado na equação (6.29), a estratégia a * é aquela com maior valor de
qa ( g + 1) , ou seja:
a* = arg maxa (q a ( g + 1))
(6.30).
Por fim, para que as probabilidades sejam calculadas no início de cada geração é
necessário definir os valores de cada recompensa, ra .
As recompensas são quantificadas fazendo uso da metodologia do fitness
improvement adotada neste trabalho, ver equações (6.20) e (6.23).
No trabalho proposto por [37], o fitness improvement de cada indivíduo é
utilizado para quantificar as recompensas de cada estratégia. Desta forma, quando o
valor da função objetivo do vetor experimental não é melhor do que o do indivíduo
correspondente, então o fitness improvement correspondente é nulo ( ηi = 0 ).
Com isso, define-se S a como um conjunto de todos os fitness improvement que
obtiveram sucesso (aqueles diferentes de zero) para cada estratégia a (a = 1, 2 e 3)
durante a geração G. Em posse de cada conjunto S a (um para cada estratégia), as
recompensas são calculadas da seguinte forma:
ra' ( g ) = media(S a )
ra ( g ) =
ra'
max r1' ; r2' ; r3'
(
(6.31).
)
(6.32).
93
6.1.2.4 ED_mod
Estas três modificações descritas nas seções acima foram incorporadas ao
método da Evolução Diferenciada tradicional com intuito de tornar o método adaptativo
em relação aos fatores de cruzamento e mutação e às estratégias de mutação. Neste
trabalho, o algoritmo da ED modificado foi chamado de ED_mod.
Desta maneira, ao início de cada iteração do ED_mod são escolhidas
aleatoriamente uma das estratégias presentes no pool para cada um dos indivíduos da
população de acordo com suas respectivas probabilidades de escolha.
Como há três tipos de estratégias de mutação presente no pool, então, após
selecionar uma estratégia para cada indivíduo, a população total é dividida em três
grupos, cada um referente a uma estratégia.
Cada grupo tem o seu próprio fator de mutação e cruzamento, onde as
adaptações destes fatores ocorrem por grupo, não havendo nenhum tipo de mistura.
Na Figura 6.2 abaixo pode ser visualizado pseudo-algoritmo do método da
ED_mod adotado neste trabalho.
94
01 Início
02
Definir Nvar , o número de variáveis
03
NP = 10.Nvar, número da população
04
Definir µ CR = 0.5; µ F = 0.5;
05
Definir o grupo Archive: A = 0;
06
Criar população inicial aleatória {popi,0 |i = 1 a NP }
07
Avalia a função objetivo para cada indivíduo.
08
Iniciando a contagem das gerações: g = 0;
09
Definir quantidade de estratégias: K = 3;
10
Definir a qualidade inicial para cada estratégia a : q a (g ) = 0;
11
Definir a probabilidade inicial para cada estratégia a : p a (g ) = 1/K ;
12
While (algum critério de parada não é atingido) do
13
g = g + 1;
14
For i = 1 a NP
15
Selecionar a estratégia 'a ' de acordo com cada p a (g )
16
Escolher aleatoriamente r 1 ≠ r 2 ≠ r 3 ≠ r 4 ≠ r 5 ≠ i
17
Gerar CR a,i = randni (µ CR,a ,0.1); F a,i = randci (µ F,a ,0.1)
18
Gerar j rand = randint(1, Nvar )
19
For j = 1 a Nvar
20
if j = j rand or rand < CRa,i
21
if a == 1
p
22
Escolher aleatoriamente pop
23
u i,j,g é gerado com a estratégia 1!
24
do grupo 100p%
else if a == 2
25
u i,j,g é gerado com a estratégia 2!
26
else if a == 3
27
u i,j,g é gerado com a estratégia 3!
28
end if
29
else
30
u i,j,g = pop i,j,g
31
end if
32
end for
33
Avaliar a função objetivo de u i,g
34
if f (pop i,g ) melhor f (u i,g )
35
pop i,g+1 = pop i,g
36
37
otimo,g
ηi = 0
else
38
pop i,g+1 = u i,g
39
Calcular η i
40
pop i,g → A; Atualizar o grupo Archive
41
CR a,i → S CRa, Guardar o valor de CR da estratégia a
42
43
F a,i → S Fa, Guardar o valor de F da estratégia a
end if
44
end for
45
Remover vetores de A , tal que, |A | ≤ NP
46
Atualizar µ CRa e µ Fa para cada estratégia: a = 1, 2 e 3
47
Calcular r a (g ) para cada estratégia: a = 1, 2 e 3
48
Atualizar q a (g ) para cada estratégia: a = 1, 2 e 3
49
50
Atualizar p a (g ) para cada estratégia: a = 1, 2 e 3
end while
51 end
Figura 6.2 - Pseudo-Algoritmo: ED_mod
95
6.2
MÉTODO DE SUPERFÍCIE DE RESPOSTA
De acordo com o que foi descrito na seção 2.3, modelos de superfície de
resposta são frequentemente usados na substituição de modelos físicos complexos, a fim
de reduzir o custo computacional em problemas de otimização.
Quando o cálculo da função objetivo é dispendioso e se tem um grande número
de variáveis de decisão, um método que pode se tornar viável é o uso de um modelo
substituto no lugar do modelo original. Estes modelos de substituição conhecidos como
modelos de superfície de resposta são construídos a partir de uma formulação analítica
adequada usando dados conhecidos sobre o modelo original.
A grande vantagem do emprego desta técnica é que uma vez construído o
modelo de aproximação, a busca pelo ponto ótimo nesta nova função é mais fácil e mais
rápida do que no modelo original para um conjunto de variáveis de decisão pertencente
ao domínio da superfície de resposta.
Assim, o primeiro passo para formular a função aproximada é elencar distintos
vetores dentro do domínio desejado e avaliar a função objetivo real de cada um destes
pontos. Estes pontos, conhecidos como pontos de interpolação, são escolhidos a partir
de uma geração aleatória uniformemente distribuída, conforme a equação abaixo:
Psr = Psri , j = LI j + Ri , j ⋅ (LS j − LI j )
Onde,
i = 1, 2, 3, ..., Nsr;
j = 1, 2, 3, ..., Nvar;
Nsr: quantidade de pontos de interpolação distribuídos dentro do domínio;
Nvar: dimensão (número de variáveis) dos pontos (vetores);
LI j : limite mínimo da variável j;
(6.33).
96
LS j : limite máximo da variável j;
Ri , j : número aleatório com distribuição uniforme entre [0, 1].
Segundo as referências [55-57] uma boa opção é utilizar no lugar de Ri , j na
equação (6.33) a sequência de Sobol, pois este é um tipo sequência de baixa
discrepância, que tende a preencher o domínio de busca de modo mais uniforme e
organizado.
Com isso, após a avaliação da função objetivo original para os vetores gerados
pela equação (6.33), as duas informações necessárias para formular analiticamente a
função aproximada estão disponíveis: a matriz composta pelos vetores gerados, Psr ; e
o vetor compostos pelas avaliações da função objetivo original, f (Psr ) .
Neste trabalho foi adotado como o modelo de aproximação o método de
superfície de resposta formulado em RBF (funções de base radial), como visto no
trabalho de [55]. As funções de base radial são assim denominadas porque apresentam
simetria radial. São funções reais que dependem somente da distância entre dois pontos.
A formulação geral de uma função de base radial é descrita da seguinte forma:
(
φ (r j ) = x j − xi
)
(6.34).
Onde,
xi : centro da função;
x j − xi : distância euclidiana entre o vetor x j e xi .
Um modelo de aproximação através de funções de base radial é conseguido
somando-se diversas funções aplicadas a cada ponto de interpolação. A equação (6.35)
mostra o modelo de aproximação s(x), adotado por [55], que foi adotado neste trabalho.
97
Neste modelo, além de fazer uso do somatório das RBF, é utilizado um polinômio.
s( x ) = ∑ α i ⋅ φ ( x − Psri ) + ∑ ∑ β j ,k ⋅ pk (x j ) + β 0 ≈ f ( x )
Nsr
M N var
i =1
k =1 j =1
(6.35).
Onde,
x = {x1 , x2 ,...,x N var −1 , x N var } ;
Nvar: dimensão do vetor x;
f(x): função objetivo original do vetor x;
α i : valor dos coeficientes constantes ligado a RBF;
β j ,k : valor dos coeficientes constantes ligado ao polinômio p k (x j ) ;
pk (x j ): polinômio de grau M, onde p k (x j ) = x kj [55].
Ao substituir a equação (6.33) na equação (6.35), uma vez que os valores de
f (Psr ) são conhecidos, os coeficientes α i e β j ,k podem ser obtidos resolvendo-se um
sistema linear de ( Nsr + M ⋅ N var + 1 ) equações, pois (segundo a referência [55]) a
equação (6.35) é sujeita as seguintes restrições:
i
⋅ p k (Psri ,1 ) = ∑ α i ⋅ (Psri ,1 ) = 0
(6.36).
i
⋅ p k (Psri ,2 ) = ∑ α i ⋅ (Psri ,2 ) = 0
(6.37).
Nsr
∑α
i =1
Nsr
Nsr
∑α
i =1
k
i =1
Nsr
k
i =1
M
Nsr
∑α
i =1
i ⋅ p k (Psri ,N var ) = ∑ α i ⋅ (Psri ,N var ) = 0
Nsr
i =1
Nsr
∑α
i =1
Onde,
k = 1, 2, ..., M;
i
=0
k
(6.38).
(6.39).
98
Nsr: quantidade de pontos de interpolação;
Nvar: dimensão dos vetores.
Assim, o sistema de equações montado fica da seguinte forma:
A
BT

B  α   f (Psr )
⋅
=
0   β   0 
(6.40).
Sendo:
 φ ( Psr1 − Psr1 ) K φ ( Psr1 − PsrNsr ) 


A=
M
O
M

φ ( PsrNsr − Psr1 ) L φ ( PsrNsr − PsrNsr )
1 (Psr1,1 )1 L (Psr1,N var )1 L (Psr1,1 )k L (Psr1,N var )k 


B = M
M
M
M
M
M
M
M

1
k
k
1(Psr )1 L (Psr
Nsr ,1
Nsr ,N var ) L (PsrNsr ,1 ) L (PsrNsr ,N var ) 

(6.43).
 β0 
 β

1


β=
 M



 β ( M ⋅ N var ) 
(6.44).
Onde,
Psr : matriz composta pelos vetores gerados pela equação (6.33);
Nsr: quantidade de pontos de interpolação;
M: grau do polinômio;
(6.42).
 α1 
α =  M 
α Nsr 
 f (Psr1 ) 

f (Psr ) = 
M

 f (PsrNsr )
Nvar: dimensão dos vetores;
(6.41).
(6.45).
99
B T : matriz transposta de B.
Ao resolver o sistema de equação, equação (6.40), obtêm-se os respectivos
coeficientes α e β . Repare que a formalização da equação (6.40) distingue a matriz
total em duas matrizes, uma correspondendo aos valores das funções de base radiais
(matriz A) e outra correspondendo apenas aos polinômios de ordem M (matriz B).
Assim, é possível facilmente obter uma função aproximada para diferentes formulações
de RBF.
Neste trabalho foram adotados três tipos de RBF: a Multiquadrática, φ M (rj ) ; a
Gaussiana, φG (r j ) ; e a Multiquadrática Cúbica, φC (r j ) . A seguir, seguem suas
respectivas formulações:
φ M (r j ) = r j2 + c 2
(6.46).
[
φG (r j ) = exp − (c ⋅ r j )2
φC (r j ) =
[r
2
j
+ c2
]
]
(6.47).
3
(6.48).
Onde,
r j : distância euclidiana entre o vetor apurado e o vetor central do método;
c: fator de forma, valor constante a ser determinado pelo usuário, sendo c > 0 .
Nesta seção foi apresentado como formular analiticamente uma função interpolada
pelo método de Superfície de Resposta. Na seção a seguir será apresentado o algoritmo de
otimização proposto por este trabalho, que fez uso do algoritmo da ED_mod (visto nas
seções anteriores) para otimizar a função aproximada.
100
6.3
ALGORITMO PROPOSTO: ED2SR
Nesta seção será apresentado o algoritmo proposto chamado ED2SR. Este
algoritmo inicia-se fazendo uso da metodologia de otimização proposta por [55], através
do algoritmo H3, porém com pequenas modificações. E finaliza a otimização fazendo
um intercâmbio entre o algoritmo H3 modificado e o método da ED_mod (visto nas
seções anteriores).
Assim, nas seções a seguir será apresentado o algoritmo H3 com as respectivas
modificações e a integração deste algoritmo (H3_mod) com o ED_mod formando o
ED2SR.
6.3.1
H3_mod
O que motivou a utilizar a metodologia de otimização proposta por [55] foi em
relação ao número de funções objetivo executadas durante cada iteração: no máximo
duas vezes.
A grande modificação ao que foi proposto por [55] foi em relação ao
procedimento utilizado para escolher a RBF ( φ M , φG ou φC , visto na seção anterior) e o
grau do polinômio (valor de M, visto na seção anterior) a serem utilizados para formular
analiticamente a função interpolada s(x), equação (6.35).
O procedimento global do algoritmo H3_mod foi baseado nos seguintes passos:
1- Carregar os pontos de interpolação responsáveis pela criação da Superfície de
Resposta e chamar de Psr. Carregar o valor da função objetivo real, f(x),
correspondente a cada vetor presente em Psr e chamar de f_Psr. Carregar os
pontos de treinamento e chamar de Ptr. Carregar o valor da função objetivo real,
f(x), correspondente a cada vetor presente em Ptr e chamar de f_Ptr;
101
2- Achar o melhor vetor (aquele com melhor valor de função objetivo real)
presente em Psr e chamar de xotimo;
3- Identificar na matriz Psr a posição do vetor mais distante de xotimo e chamar de
ifar;
4- Gerar uma função interpolada a partir do método Superfície de Resposta com a
matriz Psr e o vetor coluna f_Psr. Para escolher a RBF ( φ M , φG ou φC , visto na
seção anterior) e o grau do polinômio (valor de M, visto na seção anterior) a
serem utilizados, o algoritmo utiliza os pontos de treinamento (Ptr) para testar a
Superfície de Resposta construída para cada par ( φ , M). Desta forma, o par ( φ ,
M) que proporcionar o menor erro RMS (Root Mean Square) são os candidatos a
gerar a função interpolada chamada de g(x). Na próxima seção será detalhada a
sub-rotina chamada SR, responsável pela definição do melhor par ( φ , M);
5- Otimizar a função interpolada, g(x), com método da ED_mod (descrito na seção
6.1.2). Chamar o melhor ponto encontrado de xint. Calcular a função objetivo
real de xint: f(xint);
6- Se o valor de f(xint) for melhor do que todos os vetores presentes em Psr, então,
substituir o vetor de Psr na posição ifar por xint. Caso contrário, gerar um novo
vetor e chamar de xnovo. Avaliar a f(xnovo) e substituir xnovo na posição ifar da
matriz Psr.
Na Figura 6.3, a seguir, se encontra detalhado o algoritmo H3_mod utilizado
nesta dissertação.
102
01 Início
02 Definir Nvar , o número de variáveis
Passo 1
03
Carrega Psr , os pontos de interpolação
04
Carregar f_Psr , os valores da f (x ) de cada vetor de Psr
05
Carregar Ptr , os pontos de treinamento
06
Carregar f_Ptr , os valores da f (x ) de cada vetor de Ptr
Passo 2
07
Calcular x otimo , vetor de Psr com melhor valor de f (x )
Passo 3
Calcular i far , posição do vetor mais distante de x otimo
08
Passo 4
Gerar uma função interpolada g (x ) usando o algoritmo SR
09
g (x ) = SR(Psr , f_Psr , Ptr , f_Ptr )
Passo 5
Calcular x int , ponto ótimo da g (x ) usando o algorimo DE_mod
10
x int = DE _mod (g (x ))
11
Avalia a função objetivo real, f (x ), do vetor x int
Passo 6
12
if f (xint ) melhor do que f (xotimo )
13
Psr ifar = x int , substitui x int na posição i far
14
f_Psr ifar = f(x int )
15
16
else
Gerar i novo = SOBOL(1, Nvar ), vetor pseudo aleatório
17
Gerar x novo , novo vetor dentro dos limites de busca, utilizando i novo
18
Avalia a função objetivo real, f (x ), do vetor x novo
19
Psr ifar = x novo , substitui x novo na posição i far
20
21
f_Psr ifar = f(x novo )
end if
22 end
Figura 6.3 - Pseudo-Algoritmo: H3_mod
6.3.1.1 SR
A cada vez que as linhas de programação do algoritmo H3_mod forem
utilizadas, de acordo com toda a metodologia descrita na seção 6.2, é necessário definir
a função RBF, o grau do polinômio e o fator de forma para que a função interpolada,
g(x), seja criada.
Este trabalho não fez uso do método cross-validation1 anunciado na referência
[55] como o procedimento padrão para as escolhas do melhor conjunto RBF, grau do
polinômio e fator de forma. Este método não se mostrou robusto, pois para uma mesma
população inicial de pontos, o mesmo não garantiu um repetitividade das escolhas.
1
Para maiores informações do cross-validation consultar a referência [55].
103
Aliás, foi exatamente por este motivo que foi escolhido um método de otimização com
dupla superfície de resposta, pois para que houvesse repetitividade na escolha do melhor
conjunto (RBF, grau do polinômio e fator de forma) verificou-se a necessidade de haver
uma população independente que fizesse o papel de pontos de treinamento.
Em testes realizados foi constatado que a repetitividade depende diretamente da
quantidade e da qualidade dos pontos de treinamento. Desta forma, durante os processos
iterativos verificou-se que os pontos de treinamento também deveriam ser otimizados.
Assim, foi determinado que os pontos de treinamento fizessem o papel de uma segunda
superfície de resposta, como será descrito na próxima seção.
Para este trabalho foi criado uma sub-rotina chamada SR, que ao fornecer os
dados da superfície de resposta (Psr e f_Psr) e dos pontos de treinamento (Ptr e f_Ptr),
obtém-se a g(x). Esta rotina pode ser visualiza na Figura 6.4.
104
01 Início
02
Carrega Psr , os pontos de interpolação
03
Carregar f_Psr , os valores da f (x ) de cada vetor de Psr
04
Carregar Ptr , os pontos de treinamento
05
Carregar f_Ptr , os valores da f (x ) de cada vetor de Ptr
06
Calcular D sr , matriz contendo as distântias euclidiana dos vetor de Psr
07
Definir c = min (D sr )/Nsr , fator de forma inicial
08
For i = -1 a 5
09
grau = i
10
B = Matriz B (Psr , grau )
Montar a matriz B
11
For j = 1 a 3
12
RBF = j
13
Verificar se o fator de forma será maximizado ou minimizado!
14
Definir passo
15
RMS após = Inf
16
Do While RMS após < RMS i,j
17
RMS i,j = RMS após
18
c = c .passo , calculo do novo fator de forma
Montar a matriz A
19
A = Matriz A (D sr , RBF , c )
Construi a matriz M
20
M = [A , B | B trans , 0]
Inverter a matriz M com a método SVD
21
M inv = SVD (M )
22
α = M inv . f_Psr
Avaliar a função interpolada nos pontos de treinamento
23
Calcular D tr , distância entre os vetores de Ptr e Psr
24
A tr = Matriz A (D tr , RBF , c )
25
B tr = Matriz B (Ptr , grau )
26
M tr = [A tr , B tr ]
27
f _Ptr int = M tr . α, valor da função interpolada de Ptr
28
error = f _Ptr - f _Ptrint
Calcular o RMS do vetor error
29
RMS i,j = RMS (error )
30
α i,j = α
31
c i,j = c
32
33
34
end do while
end for
end for
Identificar as posições i e j do menor valor de RMS i,j
35
[i , j ] = min (RMS i,j )
36
RBF = j , defini a melhor Função de Base Radial
37
grau = i , defini o melhor grau do polinômio
38
α = α i,j , defini os melhores coeficientes da função interpolada
39
c = c i,j , defini o melhor fator de forma
40 end
Figura 6.4 - Pseudo-Algoritmo: SR
6.3.2
ED2SR
Esta dissertação propôs um algoritmo de otimização baseado no Método da
Evolução Diferenciada Adaptativa (o ED_mod) fazendo uso de duas superfície de
105
resposta sendo chamado de ED2SR.
O algoritmo da ED2SR é dividido em duas partes. Na primeira, este utiliza duas
Superfícies de Respostas ao mesmo tempo (Psr1 e Psr2), enquanto na segunda parte é
criada uma única superfície, a partir dos melhores vetores presentes nas superfícies
anteriores.
O objetivo da primeira parte é que ambas as superfícies sejam otimizadas pelo
algoritmo H3_mod, onde em um determinado momento a superfície Psr2 é substituída
como pontos de treinamento, testando a interpolação criada com a superfície Psr1, sendo
este papel invertido em outro momento. Ainda nesta primeira parte, as superfícies
também trocam informações, caso seja encontrado algum um ponto ótimo.
A segunda parte começa após um determinado número de funções avaliadas,
onde os vetores presentes nas matrizes Psr1 e Psr2 formam um único grupo, sendo este
grupo divido em duas matrizes. Uma nova Superfície de Resposta é criada a partir dos
melhores vetores do grupo anteriormente formado, sendo esta superfície chamada Psr.
Os vetores restantes formam o grupo dos pontos de treinamento, Ptr. Assim, nesta
segunda parte apenas a superfície Psr é otimizada com o algoritmo H3_mod e os pontos
de treinamento somente irão testar as superfícies criadas.
Ainda nesta segunda parte, quando o método não for capaz de encontrar um
novo ótimo durante 10 iterações consecutivas, o algoritmo da ED_mod é aplicado nos
pontos Psr com o objetivo de diminuir a área da superfície, deixando os pontos mais
próximos e, com isso, melhorando a eficácia do algoritmo H3_mod.
Assim, o procedimento global do algoritmo da ED2SR seguem os seguintes
passos:
1- Gerar duas populações usando a função objetivo real, f(x). Sendo estas duas
populações chamadas de Psr1 e Psr2, respectivamente;
106
2- Aplicar o algoritmo H3_mod em cada população criada no passo anterior. Na
geração de cada função interpolada, g(x), utilizar a outra população como pontos
de treinamento, respeitando os limites de interpolação, para não haver
extrapolação;
3- Verificar qual população obteve sucesso na busca de um novo ponto ótimo
durante o segundo passo. Se Psr1 obtiver sucesso, incluir o melhor vetor de Psr1
em Psr2 e aplicar o segundo passo em Psr2 e, posteriormente retirar o vetor
incluso. Se Psr2 obtiver sucesso, repetir este procedimento na população Psr1.
Repetir este terceiro passo até que nenhum ótimo seja encontrado.
4- Se o número de funções avaliadas for maior do que 200 vezes a dimensão do
problema, seguir para o próximo passo, referente à segunda parte do algoritmo
da ED2SR. Caso contrário voltar ao segundo passo.
5- Definir a população Psr com os Nsr melhores vetores de Psr1 união com Psr2.
Definir a população Ptr com os vetores de Psr1 união com Psr2, que não foram
escolhidos para Psr;
6- Aplicar o algoritmo H3_mod na população Psr, utilizando a população Ptr
apenas como pontos de treinamento. Verificar se os pontos Ptr extrapolam a
interpolação. Se houver extrapolação, gerar novos pontos de treinamento e
certificar que o número total de vetores em Ptr não seja menor do que 60% de
Nsr;
7- Caso não seja encontrado um novo ponto ótimo durante 10 iterações seguidas,
então, aplicar o algoritmo da ED_mod na população Psr, utilizando a função
real;
8- Voltar ao sexto passo até que algum critério de parada seja atingido.
107
Estes passos podem ser visualizados nos pseudo-algoritmos a seguir. A Figura
6.5 corresponde à primeira parte do ED2SR e a Figura 6.6 corresponde à segunda parte.
01 Início
02
Definir Nvar , o número de variáveis
03
Nsr 1 = 15.Nvar, número da população da 1° Superfície
04
Nsr 2 = 15.Nvar, número da população da 2° Superfície
05
Gerar r 1 = SOBOL(Nsr 1 , Nvar ), 1° matriz pseudo aleatório
06
Cria Psr 1 , 1° população inicial com a matriz r 1
07
Calcular f_Psr 1 , função objetivo de cada vetor de Psr 1
08
Gerar r 2 = rand (Nsr 2 , Nvar ), 2° matriz aleatório
09
Cria Psr 2 , 2° população inicial com a matriz r 2
10
Calcular f_Psr 2 , função objetivo de cada vetor de Psr 2
11
Calcular x otimo1 , vetor de Psr 1 com melhor valor de f (x )
12
Calcular x otimo2 , vetor de Psr 2 com melhor valor de f (x )
13
nfaval = Nsr 1 + Nsr 2 , número de funções avaliadas
14
While nfaval < 200.Nvar do
15
x otimo1_ant = x otimo1
Calcular x otimo1 , usando o algorimo H3_mod
Sendo Psr 1 os ponto de interpolação e Psr 2 os pontos de treinamento
16
17
[x otimo1 , nfaval ] = H3_mod (Psr 1 , f_Psr 1 , Psr 2 , f_Psr 2 )
x otimo2_ant = x otimo2
Calcular x otimo2 , usando o algorimo H3_mod
Sendo Psr 2 os ponto de interpolação e Psr 1 os pontos de treinamento
18
[x otimo2 , nfaval ] = H3_mod (Psr 2 , f_Psr 2 , Psr 1 , f_Psr 1 )
19
While (x otimo1 melhor x otimo1_ant ) OU (x otimo2 melhor x otimo2_ant ) do
20
Substitui x otim o1 por x otimo2 na população Psr2
21
Substitui x otim o2 por x otimo1 na população Psr1
22
if x otimo2 melhor do que x otimo2_ant
23
x otimo1_ant = x otim o1
24
[x otimo1 , nfaval ] = H3_mod (Psr 1 , f_Psr 1 , Psr 2 , f_Psr 2 )
25
end if
26
if x otimo1 melhor do que x otimo1_ant
27
x otimo2_ant = x otim o2
28
29
30
31
32
[x otimo2 , nfaval ] = H3_mod (Psr 2 , f_Psr 2 , Psr 1 , f_Psr 1 )
end if
Desfaz as trocas dos x otimo entre as populações Psr 1 e Psr 2
end while
end while
33
Nsr = 10.Nvar, número da população da Superfície Única
34
Definir Psr com os Nsr melhores vetores de [Psr 1 ; Psr 2 ]
35
Definir f _Psr com os Nsr melhores valores de [f_Psr1; f _Psr 2 ]
36
Definir Ptr com os outros vetores de [Psr 1 ; Psr2 ]
37
Definir f _Ptr com os outros valores de [f _Psr 1 ; f _Psr 2 ]
38
[x otimo , f otimo ] = 2°p_ED2SR (Psr , f_Psr , Ptr , f_Ptr , SOBOL, nfaval )
Aplicar a 2° parte do algoritmo ED2SR
39 end
Figura 6.5 - Pseudo-Algoritmo:Primeira parte do ED2SR
108
01 Início
02
Definir Nvar , o número de variáveis
03
Carrega Psr , os pontos de interpolação
04
Carregar f_Psr , os valores da f (x ) de cada vetor de Psr
05
Carregar Ptr , os pontos de treinamento
06
Carregar f_Ptr , os valores da f (x ) de cada vetor de Ptr
07
Carregar nfaval , número de função objetivo já avaliada
08
Ntr min = 6.Nvar , Quantidade mínima de pontos de treinamento
09
niter ED = 3.Nvar , número de iterações do algoritmo DE_mod
10
s_suc = 10, contagem de insucessos
11
While (algum critério de parada não é atingido) do
12
if s_suc == 10
13
s_suc = 0
14
[Psr , f_Psr , nfaval ] = DE_mod (Psr , f_Psr , niter ED )
Otimiza com DE_mod a população Psr em niter ED iterações
15
end if
16
Calcular Ntr , número de vetores de Ptr dentro do limite de interp.
17
if Ntr < Ntr min
18
Gerar (Ntr min - Ntr ) novos pontos de treinamento
19
nfaval = nfaval + (Ntr min - Ntr )
20
end if
21
x otimo_ant = x otimo
Calcular x otimo , utilizando o algorimo H3_mod
Sendo Psr os ponto de interpolação e Ptr os pontos de treinamento
22
[x otimo , f otimo , nfaval ] = H3_mod (Psr , f_Psr , Ptr , f_Ptr )
23
if x otimo melhor do que x otimo_ant
24
25
26
27
28
s_suc = 0, reinicia a contagem de insucessos
else
s_suc = s_suc + 1
end if
end while
29 end
Figura 6.6 - Pseudo-Algoritmo: segunda parte ED2SR
6.4
TESTE DE OTIMIZAÇÃO COM ED2SR:
Nesta seção serão apresentadas as comparações realizadas entre os algoritmos da
ED_mod e da ED2SR, quando aplicados a otimizações de algumas funções Benchmark
(funções de teste).
Foram escolhidas cinco funções de referência: Beale ( f1 ); Branin ( f 2 ); Easom
( f 3 ); Sphere ( f 4 ); e Griewank ( f 5 ). As três primeiras funções, Beale, Branin e Easom,
são bidimensionais, enquanto as outras podem assumir qualquer dimensão.
As formulações analíticas de cada uma destas funções estão representadas nas
equações a seguir, onde “n” corresponde a dimensão do problema:
109
(
f1 (x ) = (1,5 − x1 + x1 ⋅ x2 ) + 2 ,25 − x1 + x1 ⋅ x 22
2
) + (2,625 − x
2
1
+ x1 ⋅ x 23
)
2
(6.49).
2


5,1 ⋅ x12 5 ⋅ x1
1 

f 2 ( x ) =  x2 −
+
− 6  + 10 ⋅ 1 −
 ⋅ cos ( x1 ) + 10
2
π
4 ⋅π
 8 ⋅π 


[
f 3 ( x ) = − cos ( x1 ) ⋅ cos ( x 2 ) ⋅ exp − ( x1 − π ) − ( x 2 − π )
2
2
]
(6.50).
(6.51).
n
f 4 ( x ) = ∑ xi2
i =1
n
f 5 (x ) = ∑
i =1
n
xi2
x 
− ∏ cos i  + 1
4000 i =1
 i
(6.52).
(6.53).
A primeira função teste corresponde a função Beale com 2 variáveis de
otimização. O limite de busca foi definido dentro do intervalo − 4 ,5 ≤ x ≤ 4 ,5 , onde o
mínimo global localiza-se em para x = (3, 0,5), sendo f (x ) = 0 .
Na Figura 6.7 constam os resultados obtidos ao otimizar a função Beale com os
algoritmos ED_mod e ED2SR. Ambos os algoritmos encontraram o mínimo global,
sendo que o ED_mod necessitou calcular a função objetivo 3.620 vezes, quando o
ED2SR precisou apenas de 2.892.
O algoritmo da ED2SR se comportou muito bem no início do processo iterativo,
entretanto, necessitou fazer uso do método da Evolução Diferenciada para que o
resultado fosse refinado.
110
Figura 6.7 - Função Beale - 2D
A segunda função teste corresponde a função Branin com 2 variáveis de
otimização. O limite de busca foi definido dentro do intervalo − 5 ≤ x1 ≤ 10 e
0 ≤ x2 ≤ 15 , onde o mínimo global encontra-se em x = (-π, 12,275), (π, 2,275),
(9,42478, 2,475), sendo f ( x ) = 0,397887 .
Na Figura 6.8 constam os resultados obtidos ao otimizar a função Branin com os
algoritmos ED_mod e ED2SR. Ambos os algoritmos encontraram o mínimo global,
sendo que o ED_mod encontrou o mínimo com um número menor de função objetivo.
Entretanto, apesar do ED_mod apresentar melhor resultado, o algoritmo da ED2SR
encontrou a região muito próxima a do mínimo global com apenas 423 cálculos de
função objetivo.
111
Figura 6.8 - Função Branin - 2D
A terceira função teste analisada foi a função Easom também com 2 variáveis de
otimização. O limite de busca foi definido dentro do intervalo − 100 ≤ x ≤ 100 , onde o
mínimo global localiza-se em x = (π, π), sendo f ( x ) = −1 .
Na Figura 6.9 constam os resultados obtidos ao otimizar a função Easom com os
algoritmos ED_mod e ED2SR.
Ambos os algoritmos encontraram o mínimo global, sendo que o ED_mod
necessitou calcular a função objetivo mais de 5.000 vezes, quando o ED2SR precisou
apenas de 2.528.
112
Figura 6.9 - Função Easom - 2D
A quarta função teste analisada foi a função Sphere. Esta função foi otimizada
para 3 casos, com 2, 5 e 10 variáveis de otimização, respectivamente. O limite de busca
foi definido dentro do intervalo − 100 ≤ x ≤ 100 , onde o mínimo global localiza-se em
x = 0 , sendo f ( x ) = 0 .
Na Figura 6.10, Figura 6.11 e Figura 6.12 constam os resultados para os 3 casos
citados, respectivamente.
Nos três casos ambos os algoritmos encontraram o mínimo global. O algoritmo
da ED2SR foi extremamente superior ao ED_mod ao identificar a região do mínimo
global. Entretanto, com já observado na função Beale, o algoritmo da ED2SR necessitou
fazer uso do método da Evolução Diferenciada para que o resultado fosse refinado.
113
Figura 6.10 - Função Sphere - 2D
Figura 6.11 - Função Sphere - 5D
114
Figura 6.12 - Função Sphere - 10D
A quinta função teste analisada foi a função Griewank. Esta função também foi
otimizada para 3 casos, com 2, 5 e 10 variáveis de otimização, respectivamente. O
limite de busca foi definido dentro do intervalo − 600 ≤ x ≤ 600 , onde o mínimo global
encontra-se em x = 0 , sendo f ( x ) = 0 .
Na Figura 6.13, Figura 6.14 e Figura 6.15 constam os resultados para os 3 casos
citados, respectivamente.
Nos três casos ambos os algoritmos encontraram o mínimo global. O algoritmo
da ED2SR encontrou o mínimo global com apenas 1.018, 2.949 e 5.265 cálculos de
funções objetivo para cada um dos três casos, respectivamente.
Para esta função o algoritmo da ED2SR foi superior ao ED_mod, quando este
último necessitou calcular 6.400, 41.600 e 90.200 funções objetivo para cada um dos
três casos, respectivamente.
115
Figura 6.13 - Função Griewank - 2D
Figura 6.14 - Função Griewank - 5D
116
Figura 6.15 - Função Griewank - 10D
Diante dos resultados obtidos com as funções de teste ficou evidente o benefício
obtido quando se utilizou o algoritmo de otimização com Superfície de Resposta.
Assim, no próximo capítulo serão apresentados os resultados quando o algoritmo
proposto (o ED2SR) foi utilizado para otimizar a função objetivo proposta nesta
dissertação: o VPL5% de um empreendimento termelétrico inserido no Brasil.
117
7
UTE ÓTIMA: UM ESTUDO DE CASO
No capítulo 5 foi apresentado um modelo de avaliação econômica de uma usina
termelétrica a gás natural inserida no Sistema Interligado Nacional, segundo a
contratação pela opção de disponibilidade de energia. Este modelo foi utilizado como
função objetivo nos algoritmos de otimização descritos no capítulo 6 (o ED_mod e o
ED2SR) para otimizar a UTE-Base proposta no capítulo 4.
Nas seções a seguir será apresentado como foi modelada e quantificada a função
objetivo, bem como os resultados obtidos com a sua otimização.
7.1
MODELAGEM DO VALOR PRESENTE LÍQUIDO
Para que fosse possível realizar a simulação dos fluxos de caixa descontados
pelas estratégias propostas neste trabalho (vistas nos capítulos 3 e 5), foi necessário
quantificar as receitas obtidas nos ambientes de comercialização regulado (ACR) e livre
(ACL), as receitas obtidas no mercado de curto prazo (MCP) e os respectivos custos.
Conforme discutido no capítulo 5, segundo as equações (5.8), (5.25) e (5.26), o
VPL de um empreendimento termelétrico foi modelado em função das seguintes
variáveis:
(
VPL = f RCCEAR ,ta , RPPA ,t a ,C outros ,ta , I ,T ,tipo ,TMAa , ID , IR , Deta
)
(7.1).
Onde,
ta : período em base anual;
RCCEAR ,ta : receita no ACR (mercado cativo), em R$/ano, equação (5.12);
RPPA ,ta : receita no ACL (mercado livre), em R$/ano, equação (5.19);
C outros ,ta : outros custos (O&M, tarifas, garantia), em R$/ano, equação (5.20);
118
ID: Impostos Diretos, em %;
IR: Imposto de Renda, em %;
Deta : depreciação, em R$/ano;
I: custo do investimento, em R$, equação (5.8);
tipo: LEN tipo A-3, ou seja, valor igual a 3;
TMAa : taxa mínima de atratividade anual, em %;
T: vida útil do projeto, período total.
Os valores das principais variáveis da equação (7.1), listadas acima, precisaram
ser estimados para a elaboração do estudo de caso. Esses valores foram obtidos por
meio de coleta de dados realizada em órgãos federais, institutos de pesquisa e estudos
acadêmicos anteriores.
Entretanto, para que estas principais variáveis da equação (7.1) fossem
quantificadas de acordo com cada projeto termelétrico simulado no Thermoflow, outras
variáveis precisaram ser dimensionadas, tais como: a garantia física; a disponibilidade; a
inflexibilidade; o custo variável unitário (CVU); o custo marginal de operação (CMO);
e o preço de liquidação das diferenças (PLD) ou preço spot da energia.
Assim, nas próximas seções serão abordados as equações e valores de cada
variável presente na equação (7.1), bem como as demais variáveis auxiliares para o seu
cálculo (chamadas de variáveis comuns).
7.1.1
Avaliação das Variáveis Comuns
7.1.1.1 Vida Útil
A vida útil representa o período de tempo que o empreendedor espera utilizar os
119
ativos da UTE. Nesta dissertação foi utilizado um período de 15 anos (ou 180 meses),
conforme referência [60], sendo contados após os anos referentes à construção do
empreendimento. No estudo de caso que será apresentado, a UTE foi avaliada para um
leilão do tipo A-3, sendo os três primeiros anos reservados para a construção da térmica.
7.1.1.2 Taxa Mínima de Atratividade Anual e Mensal
A Taxa Mínima de Atratividade (TMA) é o mínimo que um investidor se propõe
a ganhar quando faz um investimento. A definição desta taxa foi baseada nas taxas
utilizadas no Plano Nacional de Energia (PNE 2030) de 10%. Assim, neste trabalho
adotou-se como taxa mínima de atratividade anual o valor de 10%.
TMAa = 10%
(7.2).
A taxa mínima de atratividade mensal foi calculada com base no valor atribuído
a taxa anual, desta forma:
TMAm =
[
12
(1 + 0 ,1) − 1]⋅100 = 0,7974 %
(7.3).
7.1.1.3 Garantia Física da UTE
Conforme equação (5.2) e identificando a primeira variável de saída do
Thermoflow como TF1 , a garantia física foi definida de acordo com a seguinte equação:
GF = GF% .TF1
(7.4).
Onde,
TF1 : potência instalada líquida da térmica, 1° output do Thermoflow, em MW;
GF% : proporção da potência, em %.
120
De acordo com a seção 3.1.3, a proporção da potência GF% é função do Custo
Variável Unitário (CVU) e da inflexibilidade operativa da UTE. Neste trabalho foi
assumido que não haveria inflexibilidade operativa, ou seja, foi considerada uma
inflexibilidade igual à zero (de acordo com a UTE Baixada Fluminense [67]).
Assim, de acordo com a metodologia de cálculo disponibilizada pela Empresa de
Pesquisa Energética (EPE) e assumindo uma inflexibilidade zero, foi possível gerar
valores de GF% em função do CVU estrutural (visto na seção 3.1.4.3). Desta forma, o
valor de GF% correspondente ao CVU da usina foi obtido a partir da interpolação linear
dos valores de GF% listados na tabela abaixo:
Tabela 7.1 - Valores de GF% versus CVU, com inflexibilidade igual a zero
CVU [R$/MWh]
5,0
10,0
15,0
20,0
25,0
30,0
35,0
40,0
45,0
50,0
%GF
100,00%
100,00%
100,00%
100,00%
100,00%
100,00%
100,00%
99,91%
98,81%
97,67%
CVU [R$/MWh]
55,0
60,0
65,0
70,0
75,0
80,0
85,0
90,0
95,0
100,0
%GF
96,26%
94,73%
93,10%
91,43%
89,53%
87,59%
85,31%
83,15%
81,20%
79,21%
CVU [R$/MWh]
105,0
110,0
115,0
120,0
125,0
130,0
135,0
140,0
145,0
150,0
%GF
77,03%
75,06%
73,03%
70,86%
69,04%
66,96%
64,88%
63,00%
60,63%
59,08%
CVU [R$/MWh]
155,0
160,0
165,0
170,0
175,0
180,0
185,0
190,0
195,0
200,0
%GF
57,32%
55,75%
54,09%
52,44%
51,01%
49,61%
47,57%
45,89%
44,67%
43,18%
CVU [R$/MWh]
205,0
210,0
215,0
220,0
225,0
230,0
235,0
240,0
245,0
250,0
%GF
42,26%
41,52%
40,67%
39,89%
39,33%
38,71%
38,16%
37,60%
37,16%
36,62%
CVU [R$/MWh]
255,0
260,0
265,0
270,0
275,0
280,0
285,0
290,0
295,0
300,0
%GF
35,76%
35,28%
34,36%
34,00%
33,61%
33,35%
33,17%
32,80%
32,46%
32,24%
7.1.1.4 Disponibilidade da UTE
A disponibilidade foi calculada de acordo com a equação (5.4), que faz uso da
potência instalada líquida (variável de saída do Thermoflow, TF1 ), assim:
Disp = TF1 ⋅ (1 − TEIF ) ⋅ (1 − IP )
(7.5).
As variáveis TEIF e IP correspondem à taxa média de indisponibilidade forçada
e a taxa de indisponibilidade programada, respectivamente. Para usinas termelétricas
121
são considerados como referência valores na ordem de 2% e 3,5%, respectivamente.
Com isso, a disponibilidade foi calculada em MW da seguinte forma:
Disp = 0,9457 ⋅ TF1
(7.6).
Onde,
TF1 : potência líquida instalada da térmica, 1° output do Thermoflow, em MW;
7.1.1.5 Inflexibilidade da UTE
Conforme já descrito, neste trabalho foi admitida uma inflexibilidade
operacional igual à zero. Assim, a equação (5.5) tomou a seguinte forma:
Inflex = 0
(7.7).
7.1.1.6 Custo Variável Unitário – CVU
De acordo com a seção 3.1.4.3 há dois tipos de CVU, sendo um de caráter
estrutural ( CVU E ) e outro para definir o despacho ( CVU OP ) da térmica. Ambos os
CVU fazem uso do mesmo fator de conversão i, equação (5.7), e do mesmo valor do
custo de operação e manutenção variável, COeM .
Conforme equação (5.7), o fator de conversão foi calculado da seguinte forma:
i = HR ⋅ CorHR ⋅ CorCOMB
(7.8).
Onde,
HR: heat rate da planta de geração, 2° output do Thermoflow, em 106Btu/MWh;
CorHR : correção do heat rate, em %;
CorCOMB : correção do preço do combustível, em %;
122
O heat rate, que corresponde à taxa de consumo energético em combustível por
potência elétrica gerada da planta, foi obtido a partir de outra variável de saída do
Thermoflow, sendo chamada de TF2 .
A correção relacionada ao valor do heat rate, CorHR , tem a função de corrigir as
perdas inerentes à transmissão e distribuição de energia elétrica. Valores típicos para
estas perdas são de 2,5% e 1%, respectivamente. Desta forma:
CorHR =
1
(1 − 0,025 ) ⋅ (1 − 0,01)
(7.9).
A correção relacionada ao preço do combustível, CorCOMB , tem a função de
diminuir o risco quando o preço do combustível na prática for maior do que o previsto.
Assim, o valor de CorCOMB correspondeu a razão entre a previsão do preço de
combustível no último ano da vida útil da UTE (no caso deste trabalho, ano 2031) e o
preço atual.
Assim, segundo a referência [64], no ano de 2031 a previsão do preço do gás
natural é de 5,630530 US$/106 Btu e o preço de referência considerado para o ano de
2013, segundo a referência [77], é de 3,94 US$/106 Btu.
CorCOMB =
5,630530
3,94
(7.10).
Com isso, substituindo as equações (7.9) e (7.10) em (7.8), obtém:
i = 1,4805 ⋅ TF2
(7.11).
Onde,
TF2 : heat rate da UTE, 2° output do Thermoflow, em 106Btu/MWh.
Segundo a referência [49] valores típicos de operação e manutenção variável são
da ordem de 10 R$/MWh. Desta forma, neste trabalho foi adotado este valor para o
custo de operação e manutenção variável, COeM .
123
7.1.1.6.1 Custo Variável Unitário de Caráter Estrutural
O CVU estrutural é aquele utilizado no cálculo do custo de operação (COP) e
custo econômico de curto prazo (CEC), partes integrantes da receita no ambiente
regulado, equação (5.12). A taxa de câmbio e o preço do combustível de referência para
calcular o CVU estrutural são valores referenciais, disponibilizados pela EPE na época
de cada leilão. Assim, de acordo com a referência [77], a taxa de câmbio e preço do
combustível foram de 1,9484 R$/US$ e 3,94 US$/106 Btu, respectivamente.
Assim, substituindo esses valores e o valor de COeM na equação (5.6), obteve-se:
CVU E = 11,3655 ⋅ TF2 + 10
(7.12).
Onde,
TF2 : heat rate da UTE, 2° output do Thermoflow, em 106Btu/MWh.
7.1.1.6.2 Custo Variável Unitário de Caráter Operacional
O CVU referente ao despacho é aquele que define se a UTE entra em operação
ou não. Para calcular o este CVU, foi levado em consideração o valor do dólar médio.
Para o preço do combustível foi feito uma análise em relação às previsões dos seus
valores durante a vida útil do projeto termelétrico.
Pela razão da sua volatilidade, o câmbio considerado foi de 2,1303 R$/US$, que
corresponde à média dos últimos 6 meses. Estes valores estão disponíveis no site do
Banco Central do Brasil [78].
O preço do combustível foi definido como a média dos valores previstos para os
anos entre 2017 e 2031 (vida útil do projeto), sendo este valor de 4,8267 US$/106Btu.
Os valores de cada previsão estão disponíveis na referência [64].
124
Assim, o CVU de caráter operacional de cada UTE foi calculado da seguinte
forma:
CVU OP = 15,2231⋅ TF2 + 10
(7.13).
Onde,
TF2 : heat rate da UTE, 2° output do Thermoflow, em 106Btu/MWh.
7.1.1.7 Custo Marginal de Operação – CMO
O CMO é aquele utilizado no cálculo do custo de operação (COP), equação
(3.10), e custo econômico de curto prazo (CEC), equação (3.11), partes integrantes da
receita no ambiente regulado, equação (5.12). Estes valores são disponibilizados pela
Empresa de Pesquisa Energética [62], na época do leilão. Neste trabalho foram
utilizadas as séries de CMO disponibilizadas na referência [79], correspondente ao
leilão ocorrido em agosto de 2013.
Para o cálculo do CEC foram utilizados valores de CMO limitados ao um Preço
de Liquidação das Diferenças máximo e outro mínimo, PLDmáx e PLDmin . De acordo
com a mesma referência [79] estes valores são de 780,03 R$/MWh e 14,13 R$/MWh,
respectivamente.
7.1.1.8 Preço de Liquidação das Diferenças – PLD
Os PLD correspondem ao preço da energia no mercado de curto prazo. Para este
estudo, os preços spot foram determinados com base nos CMO (disponibilizados pela
Empresa de Pesquisa Energética [62]), conforme descrito detalhadamente na seção 5.4.
Na realidade a EPE apenas disponibiliza valores futuros de CMO para os
primeiros 60 meses da vida útil do projeto. Como o VPL foi calculado para 180 meses
125
(15 anos), então, para os 120 meses restantes, esta dissertação utilizou o mesmo
mecanismo de cálculo visto nas referências [1, 2], quando os últimos 12 meses dos 60
iniciais foram repetidos para os 120 meses restantes.
Na realidade isso não quer dizer que as séries (ou cenários) de PLD foram
repetidas. O que realmente foi repetido foi a curva de distribuição acumulada empírica
dos CMO. Assim, as curvas de distribuição acumulada dos últimos 10 anos foram iguais
à distribuição acumulada do último ano disponibilizado pela EPE.
Desta forma, os cenários de PLD para cada um dos 180 meses foram obtidos a
partir da equação (5.32), sendo esta equação limitada a um valor máximo de 780,03
R$/MWh e um valor mínimo de 14,13 R$/MWh. Desta forma:
PLDi ,c = min(max(cmoi ,c ; 14,13 ); 780,03)
(7.14).
cmoi ,c = quantile(CMOi , ri ,c )
(7.15).
Sendo,
Onde,
i: corresponde o mês analisado, i = 1 a 180;
c: cenário criado, sendo c = 1 a 10.000;
CMOi : conjunto dos 2.000 valores disponibilizado pela EPE no mês i.
ri ,c : número aleatório a partir de uma distribuição uniforme com intervalo [0,1]
Desta maneira foram gerados aleatoriamente 10.000 valores de PLD para cada
um dos 180 meses (que corresponde ao tempo de vida útil do projeto, 15 anos), a partir
de cada uma das 180 curvas de distribuição de probabilidade acumulada empírica
(obtidas pelos 2.000 cenários de CMO de cada um dos 180 meses).
126
7.1.2
Receita da Comercialização de Energia no Ambiente Regulado
De acordo com a equação (5.12), a RCCEAR ,ta foi definida em função das seguintes
variáveis:
RCCEAR,ta = f (8760, x ,GF , ICB,COP,CEC)
(7.16)
Onde,
ICB: índice de custo benefício, em R$/MWh;
COP: custo de operação, em R$/ano;
CEC: custo econômico de curto prazo, em R$/ano;
GF: garantia física, em MW médios;
8760: número de horas no ano;
x: fração da GF destinada ao Ambiente de Contratação Regulado.
7.1.2.1 Valor do ICB
Conforme descrito nas seções 3.1.2 e 3.1.3, na prática o ICB é utilizado para a
ordenação
econômica
de
empreendimentos
de
geração
termelétrica
e,
consequentemente, como critério de contratação por meio de contratos de
disponibilidade de energia elétrica, sob o regime de autorização por 15 anos.
De acordo com a referência [62], o valor médio do ICB do Leilão de Energia
Nova (LEN) do tipo A-3 ocorrido em 2011 foi de 102,07 R$/MWh. Neste mesmo ano,
houve um LEN do tipo A-5, onde o ICB médio foi de 102,18 R$/MWh. Desta forma,
neste estudo foi adotado como valor do ICB a média destes dois valores citados:
ICB =
(102 ,07 + 102 ,18) = 102 ,125
2
(7.17)
127
7.1.2.2 Valor do COP e CEC
Antes de modelar o Custo de Operação (COP) e o Custo Econômico de Curto
Prazo (CEC), foi necessário definir a geração. Desta forma, de acordo com a equação
(3.9), o valor da geração, Gc ,m , foi obtida da seguinte forma:
CMOc ,m ≥ CVU E ∴ Gc ,m = Disp

CMOc ,m < CVU E ∴ Gc ,m = Inflex
(7.18).
Onde,
c: corresponde ao índice de cada cenário hidrológico (1 a 2000);
m: corresponde ao índice de cada mês (1 a 60);
CMOc ,m : custo marginal de operação, definido na seção 7.1.1.7, em R$/MWh;
CVU E : CVU estrutural, equação (7.12), em R$/MWh;
Inflex : inflexibilidade operativa, equação (7.7), em MW;
Disp : disponibilidade, equação (7.6), em MW.
7.1.2.2.1 Custo de Operação
De acordo com a equação (3.10), o Custo de Operação foi calculado da seguinte
forma:
∑∑ CVU ⋅ (G
m
COP =
c
i =1 j =1
E
c ,m
− Inflex )
m⋅c
⋅ 8760
Onde,
c: corresponde ao índice de cada cenário hidrológico (1 a 2000);
m: corresponde ao índice de cada mês (1 a 60);
CVU E : CVU estrutural, equação (7.12), em R$/MWh;
(7.19).
128
Gc ,m : geração da térmica, para cada cenário e mês, em MW, equação (7.18);
Inflex : inflexibilidade operativa, equação (7.7), em MW;
8760: número de horas no ano.
7.1.2.2.2 Custo Econômico de Curto Prazo
O Custo Econômico de Curto Prazo foi calculado a partir da equação (3.11),
onde:
∑∑ CMO ⋅ (GF − G )
m
CEC =
c
i =1 j =1
*
c ,m
c ,m
m⋅c
⋅ 8760
(7.20).
Conforme descrito na seção 7.1.1.7, a parcela CMOc*,m corresponde ao valor do
CMOc ,m limitado ao valor de PLDmáx e PLDmin . Assim:
CMOc*,m = min(max (CMOc ,m ; 14,13 ); 780,03)
(7.21).
Onde,
c: corresponde ao índice de cada cenário hidrológico (1 a 2000);
m: corresponde ao índice de cada mês (1 a 60);
GF: garantia física da UTE, em MW, equação (7.4);
Gc ,m : geração da térmica, para cada cenário e mês, em MW, equação (7.18);
8760: número de horas no ano.
7.1.3
Receita da Comercialização de Energia no Ambiente Livre
De acordo com a equação (5.16) e (5.19), a Receita do Contrato de
Comercialização de Energia no Ambiente Livre ( RPPA ,ta ) foi definida em função das
seguintes variáveis:
129
R PPA ,ta = f (Disp ,GF , PPPA , PLDt ,CVU t , ht ,TMAm )
(7.22).
Onde,
PPPA : preço da energia contratada no PPA, em R$/MWh;
PLDt : preço spot no mês t, em R$/MWh, equação (7.14);
CVUt : custo variável unitário da UTE no mês t, em R$/MWh;
Disp: disponibilidade da UTE, em MW, equação (7.6);
GF: garantia física da UTE, em MW, equação (7.4);
TMAm : taxa mínima de atratividade mensal, em %, equação (7.3);
ht : número de horas no mês t;
TMAm : taxa mínima de atratividade mensal, em %.
7.1.3.1 Preço da Energia Contratada
O preço da energia contratada ( PPPA ) refere-se ao valor pago por MWh em
contratos de longo prazo de compra de energia elétrica incentivada no ambiente livre.
Na prática, por se tratarem de acordos bilaterais, os valores são definidos livremente
pelas partes e os agentes envolvidos não precisam informar os preços negociados.
As referências [49, 50] indicaram um valor de 155 R$/MWh, com isso, neste
trabalho foi adotado para preço de energia negociado no ambiente livre o valor de 155
R$/MWh.
7.1.3.2 Custo Variável Unitário no mês t
Neste trabalho foi adotado para todos os meses ao longo da vida útil do projeto o
mesmo valor do Custo Variável Unitário de caráter operacional, equação (7.13).
130
7.1.3.3 Número de Horas no Mês
Para o número de horas no mês, ht , foi adotado um valor médio, que
corresponde ao número de horas em um ano divido por 12 meses, assim:
ht =
7.1.4
8760
= 730
12
(7.23).
Custos Diversos
7.1.4.1 Investimento
De acordo com a equação (5.8), o investimento necessário foi definido em
função das seguintes variáveis:
I = f (CTF ,CCX ,e2013 ,tipo ,TMAa )
(7.24).
Onde,
CTF : custo do investimento, correspondendo ao TF3 , 3° output do Thermoflow,
em US$;
C CX : custo de conexão ao SIN, em R$;
e2013 : taxa de câmbio de 2,1303, R$/US$, definido na seção 7.1.1.6.2;
tipo: LEN tipo A-3, ou seja, valor igual a 3;
TMAa : taxa mínima de atratividade anual, em %, equação (7.2).
O custo do empreendimento termelétrico disponibilizado pelo Thermoflow não
leva em consideração o valor real do custo de conexão ( C CX ) ao sistema elétrico
praticado no Brasil. Para uma termelétrica com potência instalada na ordem de 500 MW
131
a experiência operacional em empresas de energia indica um valor de R$ 40.000.000.
7.1.4.2 Outros Custos
Baseado na equação (5.20) a variável C outros ,ta , referente a outros custos, foi
definida em função das seguintes variáveis:
C outros ,ta = f (Pot bruta ,CO&M , fixo ,Ctusd / tust ,CTF ,e2013 )
(7.25).
Onde,
Potbruta : potência bruta instalada, correspondendo ao TF4 , 4° output do
Thermoflow, em MW;
CO&M , fixo : custo fixo de operação e manutenção, em R$/MWh;
Ctusd / tust : tarifa do uso dos sistemas de distribuição/transmissão, em R$/MWh;
CTF : custo do investimento, TF3 , 3° output do Thermoflow, em US$;
C CX : custo de conexão ao SIN, em R$, seção 7.1.4.1.
Os custos fixos são custos considerados constantes durante a vida útil de um
empreendimento térmico, como por exemplo, os encargos trabalhistas.
Não há na literatura valores disponíveis para o custo fixo de operação e
manutenção, entretanto a experiência operacional em empresa de energia indica um
valor em torno de um terço do valor considerado para o custo variável, que foi de 10
R$/MWh (seção 7.1.1.6.1). Neste trabalho foi assumido um valor de 3 R$/MWh para o
custo fixo de operação e manutenção, CO&M , fixo .
De acordo com as referências [49, 50] foi considerado para tarifa Ctusd / tust um
valor de 5 R$/MWh;
132
7.1.4.3 Impostos Diretos e de Renda
Segundo a referência [50], os impostos aplicados diretamente às receitas brutas
(ID) são da ordem de 9,75%. Para o imposto de renda (IR) considerou-se uma carga
tributária de 34%.
7.1.4.4 Depreciação
Depreciação é um custo ou a despesa decorrente do desgaste ou da
obsolescência dos ativos imobilizados (máquinas, veículos, móveis, imóveis e
instalações) da empresa.
Ao longo do tempo, com a obsolescência natural ou desgaste com uso na
produção, os ativos vão perdendo valor. Essa perda de valor é apropriada pela
contabilidade periodicamente até que esse ativo tenha valor reduzido à zero.
No Brasil, em termos contábeis, o cálculo da depreciação deve obedecer aos
critérios determinados pelo governo, através da Secretaria da Receita Federal [80], que
estipula o prazo de 5, 10 ou 25 anos para depreciar os ativos.
Neste trabalho adotou-se uma depreciação linear em dez anos, onde todo o valor
investido na planta de geração termelétrica foi depreciado e utilizado para abater o
imposto de renda, conforme visto na equação (5.24).
Dito isto, a depreciação anual foi calculada da seguinte forma:
Deta =
(TF3 ⋅ e2013 + CCX )
10
TF3 : custo do investimento, 3° output do Thermoflow, em US$;
C CX : custo de conexão ao SIN, em R$, seção 7.1.4.1;
e2013 : taxa de câmbio de 2,1303, R$/US$, definido na seção 7.1.1.6.2.
(7.26).
133
Repare que o valor da depreciação somente foi aplicado aos primeiros 10 anos.
Como o projeto de investimento considerado neste trabalho tem vida útil de 15 anos,
então, o valor da depreciação para os 5 anos restantes foi considerado zero.
7.2
VALOR PRESENTE LÍQUIDO DA UTE-BASE
Para que o calculo do VPL fosse realizado, inicialmente foi definido quanto da
energia total gerada pela UTE seria negociada no mercado cativo (ACR) e no mercado
livre (ACL), respectivamente.
Segundo as referências [2, 8], um empreendimento termelétrico não pode
negociar menos dos 70% de sua energia assegurada (garantia física) no mercado cativo
(Ambiente de Contratação Regulado - ACR). Com isso, para este estudo foi
estabelecida que a UTE comercializar-se no ACR 85% da sua garantia física, sendo o
restante (15%) negociado no ambiente livre (ACL).
Assim, a partir do uso de todos os dados descritos nas seções anteriores e das
quantidades de energias negociadas em cada ambiente de contratação, a quantificação
do Valor Presente Líquido ficou em função apenas das variáveis de saída do
Thermoflow, ou seja:
VPL = f (TF1 ,TF2 ,TF3 ,TF4 )
Onde,
TF1 : potência líquida instalada da térmica, em MW;
TF2 : heat rate da UTE, em 106Btu/MWh;
TF3 : custo do investimento, em US$;
TF4 : potência bruta instalada, em MW;
(7.27).
134
Com isso, foi possível calcular o VPL da UTE-Base (modelada no Thermoflow,
capítulo 4, seção 4.2). Recapitulando, a UTE-Base foi projetada com as seguintes
principais características:
1- Geração em ciclo combinado do tipo 2-2-1: dois TGG; duas HRSG; e um TGV;
2- Dois TGG da General Eletric, modelo 7FA (modelo utilizado nas térmicas
vencedoras do LEN 2011);
3- Duas HRSG gerando vapor em três níveis de pressão, projetadas segundo menor
temperatura da chaminé: 90°C;
4- Turbina a vapor com três níveis de pressão, com reaquecimento e indução de
vapor de baixa pressão;
5- Condensador com desaerador integrado;
6- Torre de Resfriamento com resfriamento a ar.
Ao simular este projeto no Thermoflow foi obtido um valor de 506 MW de
potência líquida instalada ( TF1 ), um heat rate de 6,1991 106Btu/MWh ( TF2 ), um custo
referencial de US$ 355.650.528,00 ( TF3 ) e uma potência instalada bruta de 522 MW
( TF4 ).
Substituindo os valores de TF1 , TF2 , TF3 e TF4 na equação (7.27) e levando em
consideração toda a filosofia de cálculo do VPL apresentada na seção 7.1 deste capítulo,
obteve-se 10.000 possíveis valores de VPL para a UTE-Base, um para cada cenário de
PLD, conforme apresentado na Figura 7.1, a seguir:
135
Cenário
PLD t
...
3
T
VPL
1
136,45
14,13
14,13
... 27,67
79,66
VPL 1
2
229,55
414,85
14,13
... 33,43
78,63
VPL 2
3
14,13
182,42
14,13
... 68,80
14,13
VPL 3
...
...
...
...
...
T -1
...
2
...
1
9.998
220,03
14,13
76,97
... 50,48
21,13
VPL 9.998
9.999
35,96
14,13
14,13
... 14,13
85,94
VPL 9.999
10.000
38,05
96,29
14,13
... 54,79
211,55
VPL 10.000
Figura 7.1 - Cenários de VPL
A partir destes 10.000 valores (cenários) de VPL foi possível traça a curva de
probabilidade acumulada da UTE-Base, como pode ser visto na Figura 7.2. Na Tabela
7.2 estão identificados os valores dos cenários mais otimista, mais pessimista, o ponto
de VPL igual à zero, o VPL médio e o VPL5%.
Tabela 7.2 - VPL da UTE-Base
Valor Presente Líquido da
UTE-Base
VPL Otimista
VPL Pessimista
VPL igual a zero
VPL médio
VPL5% - Função Objetivo
6
VPL [10 R$]
21,851
-57,721
0
-12,497
-28,562
Probabilidade de Ocorrer
Valores Menores
100%
0%
90,38%
49,42%
5%
136
Figura 7.2 - Probabilidade Acumulada: UTE-Base
Analisando a Figura 7.2 em conjunto com a Tabela 7.2 foi possível notar que há
90,38% de chance do VPL da UTE-Base ser menor do zero, ou seja, o risco de um
empreendimento com as caraterística da UTE-Base não trazer retorno ao investidor é de
90,38%.
A visão otimista, apesar de propiciar um VPL de R$ 21,851x106, no ponto de
vista de um investidor não pode ser levada em consideração, pois a probabilidade deste
cenário ocorrer é praticamente nula.
Diante destes resultados, um investidor não pode apenas levar em consideração
o tipo de investimento. Além de escolher o melhor investimento, segundo as regras
impostas pelo Setor Elétrico Brasileiro, é necessário averiguar se o projeto que se deseja
investir é realmente o melhor tecnicamente.
Não foi encontrado na literatura este tipo de abordagem. Desta forma, baseado
nas configurações da UTE-Base, o objetivo deste trabalho foi otimizar os parâmetros do
processo de geração desta térmica, em busca de um projeto tecnicamente mais vantajoso
137
segundo as regras impostas pelo Setor Elétrico Brasileiro.
Na realidade, o ideal é mitigar ao máximo o risco de um projeto. Por isso, foi
adotado como função objetivo o valor do VPL com o risco de apenas 5% de não
ocorrer.
7.3
OTIMIZANDO A UTE-BASE
A UTE-Base foi projetada no Thermoflow segundo a máxima eficiência térmica
possível, onde todos os parâmetros do processo de geração foram escolhidos pelo
próprio software, conforme descrito na 4.2.1 do capítulo 4.
Para otimizar o projeto térmico da UTE-Base alguns parâmetros do processo de
geração foram escolhidos para serem prescritos e, com isso, modificar o projeto de
acordo com os novos valores.
Assim, antes de iniciar o processo de otimização, foram escolhidos os
parâmetros do processo de geração que mais impactam no dimensionamento dos
principais equipamentos, tais como, a turbina a vapor e as caldeiras recuperadoras de
calor. Como na planta Base há três níveis de pressão, foram escolhidas como variáveis
de otimização as pressões, temperaturas e vazões mássicas de vapor de alta pressão e de
reaquecimento, totalizando seis variáveis.
A cada novo conjunto de valores para estas 6 variáveis o Thermoflow projeta um
distinto projeto. Cada novo projeto tem equipamentos com dimensões e características
diferentes, proporcionando, por exemplo, novos custos, nova capacidade de geração,
eficiência térmica distinta, dentre outras modificações.
Durante o processo de otimização não foi modificado o tipo de turbina a gás
escolhida para o UTE-Base (seção 4.2.1). O simulador de processo Thermoflow não
permite modificar o tipo desta máquina, deixando livre para ser otimizado apenas o
138
ciclo Rankine do ciclo combinado (o ciclo Brayton continua o mesmo durante todo o
processo iterativo).
Assim, antes de iniciar o processo de otimização foi necessário definir o limite
de busca (domínio de busca) de cada variável de otimização escolhida.
Segundo a referência [81], os pares pressão e temperatura usuais na entrada de
uma turbina a vapor podem ser visualizados abaixo:
Pressão do
Vapor [bar]
Acima de 17,24
1 a 200
201 a 315
Temperatua do Vapor [°C]
316 a 399 400 a 450 451 a 482
483 a 510 511 e Maior
17,25 a 41,37
41,38 a 62,01
62,02 a 103,43
103,44 e Maior
Figura 7.3 - Par Pressão e Temperatura na Entrada da Turbina a Vapor [81]
Para a pressão do vapor de alta, de acordo com a Figura 7.3, foram estabelecidos
como limites máximo e mínimo os valores de 131 e 100 bar, respectivamente. De
acordo com a mesma figura, para a temperatura do vapor de alta pressão foram
estabelecidos como limites máximo e mínimo os valores de 560 e 450 °C,
respectivamente.
Para a temperatura do vapor de média pressão foram adotados os mesmo limites
do vapor de alta pressão. Entretanto, para determinar as limites de busca da pressão de
média e das vazões mássicas de vapor, foram realizadas várias simulações no
Thermoflow em torno dos valores verificados na UTE-Base, sendo de 27,6 bar para a
pressão de média, 362,4 t/h para vazão mássica de alta pressão e 433,3 t/h para vazão
mássica de reaquecimento.
Os valores da pressão de média e das vazões de vapor foram variados até que o
simulador Thermoflow indicasse os limites físicos dos materiais (utilizados nos
equipamentos e nas tubulações em geral) ou indicasse o limite da geração de vapor, pois
139
esta grandeza está diretamente ligada a qualidade e a quantidade dos exaustos da turbina
a gás utilizada.
Desta forma, os limites máximos e mínimos para pressão de média foram de 35
e 20 bar, respectivamente. Foi definido para a vazão mássica de alta pressão o limite de
busca máximo de 400 t/h e mínimo de 300 t/h. Os limites de busca para a vazão mássica
referente ao reaquecimento foram de 470 t/h e 370 t/h, respectivamente.
Na Tabela 7.3 a seguir, e de acordo com Figura 7.3, estão disponíveis os limites
de busca de cada variável de otimização adotada neste trabalho.
Tabela 7.3 - Limites de Busca das Variáveis de Otimização
Variáveis de Otimização
Pressão de Vapor de Alta [bar]
Temperatura de Vapor de Alta [°C]
Pressão de Vapor de Reaquecimento [bar]
Temperatura de Vapor de Reaquecimento [°C]
Vazão Mássica de Vapor de Alta Pressão [t/h]
Vazão Mássica de Vapor para Reaquecimento [t/h]
7.3.1
Limites de Busca
Máximo Mínimo
131
100
560
450
35
20
560
450
400
300
470
370
Otimizando Com o Algoritmo da ED_mod
Ao otimizar a UTE-Base segundo o algoritmo da ED_mod (apresentado no
capítulo 6, seção 6.1.2), obteve-se o resultado ótimo de R$ 23,449x106.
O valor do VPL5% da UTE-Base (ultima linha da Tabela 7.2) foi de R$ (28,562x106). Diante disso, é notório que ao otimizar tecnicamente o projeto térmico
houve um ganho de R$ 52,011x106, pois além de tornar o valor da função objetivo
positivo ainda o cresceu em R$ 23,449x106.
Na Tabela 7.4 e Tabela 7.5 são apresentados os valores das variáveis de
otimização e os valores das variáveis de saída do Thermoflow TF1 , TF2 , TF3 e TF4 ,
140
respectivamente. Alguns dos valores apresentados também podem ser visualizados na
Figura 7.4, que mostra o resultado da UTE-Ótima simulado no Thermoflow.
Tabela 7.4 - Variáveis de Otimização: UTE-Base vrs UTE-Ótima
Variáveis de Otimização
Pressão de Vapor de Alta [bar]
Temperatura de Vapor de Alta [°C]
Pressão de Vapor de Reaquecimento [bar]
Temperatura de Vapor de Reaquecimento [°C]
Vazão Mássica de Vapor de Alta Pressão [t/h]
Vazão Mássica de Vapor para Reaquecimento [t/h]
UTE
Base
124,0003
566,0001
27,6
566,0001
362,4166
433,3301
Ótima
100,0037
483,5489
24,9089
482,9103
398,2997
390,1601
141
Figura 7.4 - Fluxograma da UTE-Ótima
142
Tabela 7.5 - Variáveis de Saída do Thermoflow: UTE-Base vrs UTE-Ótima
Variáveis de Sáida do Thermoflow
UTE
Base
505,8266
Ótima
476,0897
6,1991
6,5863
TF 3 - Custo do Investimento [10 US$]
355,651
292,372
TF 4 - Potência Bruta [MW]
Índice de Custo de Instalação [US$/kW]
6
Função Objetivo - VPL 5% [10 R$]
521,8339
681,540
-28,562
490,9094
595,572
23,449
TF 1 - Potência Líquida [MW]
6
TF 2 - Heat Rate [10 Btu/MWh]
6
Na terceira linha da Tabela 7.5 constam os valores dos custos (na referência
Thermoflow) da UTE-Base e UTE-Ótima, respectivamente. Para que fosse possível
comparar os custos de cada projeto, foi necessário dividir estes custos pelos respectivos
valores de potência bruta instalada em kW.
Na quinta linha da Tabela 7.5 constam os índices do custo de instalação de cada
projeto. Percebe-se que para cada kW instalado o UTE-Ótima economiza
aproximadamente US$ 86, apesar de ser menos eficiente (segundo linha da Tabela 7.5)
e proporcionar uma potência líquida menor do que a UTE-Base (primeira linha da
Tabela 7.5).
Estas observações foram corroboradas ao analisar a Tabela 7.4, a qual apresenta
os valores das variáveis de otimização de cada projeto.
Apesar do projeto da UTE-Base fazer uso de uma quantidade menor de vazão
mássica de vapor de alta pressão (362,4166 t/h), a quantidade total de vapor que a
caldeira recuperadora de calor deste projeto evapora e superaquece é maior do que a do
projeto da UTE-Ótima. A UTE-Base evapora e superaquece 433,3301 t/h de vapor,
enquanto a UTE_Ótima evapora e superaquece 398,2997 t/h da vapor. Este é o primeiro
indício do projeto da UTE-Ótima ser mais barato, pois a cadeira deste projeto produz
uma quantidade menor de vapor, sendo refletido no custo do equipamento.
143
O segundo indício pode ser notado ao reparar na quantidade de níveis de pressão
de vapor. A UTE-Base foi projetada com três níveis de pressão, entretanto, o nível de
média pressão do projeto da UTE-Ótima foi eliminado. Isso pode ser notado ao verificar
que a vazão mássica de reaquecimento no projeto ótimo (390,1601 t/h) é menor do que
a vazão total de vapor (398,2997 t/h), não havendo nenhuma vazão mássica adicional
sendo evaporada e superaquecida no nível de média pressão. Esta observação também
pode ser vista no fluxograma do processo de geração da UTE-Ótima (Figura 7.4), onde
mostra apenas dois níveis de pressão na caldeira, sendo um referente a baixa pressão e o
outra a alta pressão de vapor. Certamente o custo da caldeira recuperadora de calor da
UTE-Ótima é menor.
Na Tabela 7.6 a seguir estão listados os custos dos principais equipamentos (na
referência Thermoflow). Na terceira linha desta tabela estão referenciados os custos das
cadeiras recuperadoras de calor de cada projeto. A diferença em relação a este item
corresponde a aproximadamente US$ 26x106.
Outra diferença nos custos pode ser notada ao observar a segunda linha da
Tabela 7.6. Como as condições do vapor (pressão e temperatura) de alta pressão (Tabela
7.4, primeira e segunda linha) são menos severas no projeto da UTE-Ótima, isso é
refletido no custo da turbina a vapor.
O custo total de todos os equipamentos do projeto da UTE-Ótima é
aproximadamente US$ 168x106. A última linha da Tabela 7.6 mostra que a diferença do
custo total de todos os equipamentos entre os dois projetos é de aproximadamente US$
32x106.
Na terceira linha da Tabela 7.5 estão listados os custos dos investimentos (na
referência Thermoflow) da UTE-Base e da UTE-Ótima, respectivamente. A diferença
entre os dois investimentos corresponde a US$ 63,279x106. Assim, é possível observar
144
que a economia obtida somente entre os principais equipamentos (Tabela 7.6)
corresponde a 50,1% da total economia obtida no custo do investimento. Esta
observação mostra que os parâmetros do processo de geração (listados na Tabela 7.3,
pág. 139) escolhidos como variáveis de otimização impactaram diretamente no
dimensionamento dos principais equipamentos, o que certamente proporcionou a
geração de distintos projetos durante o processo iterativo na busca pelo melhor projeto
de geração segundo a função objetivo escolhida.
Tabela 7.6 - Custo dos Principais Equipamentos - Referência Thermoflow
Custo de Referência dos Principais Equipamentos
Principais Equipamentos
Base
Turbina
TurbinaaaGás
Gás[US$]
(2x)
98,154
Turbina
TurbinaaaVapor
Vapor[US$]
28,630
Caldeira
Caldeira[US$]
Recuperadora de Calor (2x)
47,113
Condensador
Condensador[US$]
2,516
Outros
23,278
Outros
Total [US$]
199,691
6
UTE [10 US$]
Ótima
Diferença
98,154
0
24,069
4,561
20,939
26,174
2,336
0,180
22,527
0,751
168,025
31,666
Além do projeto ótimo se mostrar mais barato e tecnicamente mais simples, pois
eliminou um nível de pressão, ele também se adequou melhor aos moldes impostos pelo
Setor Elétrico Brasileiro, pois o risco ao investimento foi mitigado, conforme será visto
a seguir.
A Figura 7.5 mostra o perfil de probabilidade dos possíveis VPL da UTE-Ótima
construído a partir dos 10.000 cenários de VPL deste projeto ótimo. Neste gráfico estão
identificados os pontos dos cenários mais otimista, mais pessimista, de VPL igual à
zero, de VPL médio e o valor da função objetivo do projeto ótimo, o VPL5%.
145
Figura 7.5 - VPL Ótimo pelo Método da ED_mod
Na Tabela 7.7, estão identificados os valores de cada uma dos pontos presentes
na Figura 7.5.
Tabela 7.7 - VPL da UTE-Ótima pelo Método da ED_mod
Valor Presente Líquido da
UTE-ÓTIMA - ED _mod
VPL Otimista
VPL Pessimista
VPL igual a zero
VPL médio
VPL5% - Função Objetivo
6
VPL [10 R$]
70,272
-7,134
0
38,649
23,449
Probabilidade de Ocorrer
Valores Menores
100%
0%
0,02%
49,12%
5%
Repare que todos os valores são maiores do que os apresentados para a UTEBase, Tabela 7.2 (pág. 135). O mais interessante foi que durante o processo de
otimização a estratégia de investimento escolhida se manteve constante durante todo o
processo iterativo. Entretanto, o risco foi mitigado. Esta conclusão pode ser observada
ao analisar a probabilidade de ocorrer VPL menores do que zero nos dois projetos.
Enquanto o projeto da UTE-Base promove uma probabilidade de 90,38% (Tabela 7.2),
146
o projeto ótimo garante que a probabilidade de ocorrer VPL menor do zero é de 0,02%
(Tabela 7.7).
Na Tabela 7.8 a seguir é feito um paralelo entre os pontos evidenciados nas
curvas de probabilidade acumulada dos dois projetos. É notório que, além do projeto
ótimo mitigar o risco, garante VPL maiores do que o projeto da UTE-Base.
Tabela 7.8 - Pontos da Curva - UTE-Base vrs UTE-Ótima
Valor Presente Líquido
VPL Otimista
VPL Pessimista
VPL médio
VPL5% - Função Objetivo
6
6
UTE-Base [10 R$]
UTE-Ótima [10 R$]
21,851
-57,721
-12,497
-28,562
70,272
-7,134
38,649
23,449
Em relação ao desempenho do algoritmo da ED_mod durante a otimização, de
acordo com a Figura 7.6, este algoritmo necessitou de aproximadamente 30.000
cálculos da função objetivo para que a convergência fosse atingida, ressaltando que
cada função objetivo corresponde a uma rodada do software Thermoflow, que consumiu
de 4 a 5 segundos.
147
Figura 7.6 - Curva de Convergência: Método da ED_mod
O algoritmo da ED_mod necessitou de uma quantidade excessiva de cálculos de
função objetivo. Por este motivo, nesta dissertação foi modelado e proposto outro
algoritmo que fizesse uso do método da Superfície de Resposta durante os processos
iterativos (o ED2SR), com intuito de diminuir o número de simulações do Thermoflow.
7.3.2
Otimizando Com o Algoritmo da ED2SR
Ao otimizar a UTE-Base segundo o algoritmo da ED2SR (apresentado no
capítulo 6, seção 6.3.2), obteve-se o resultado ótimo de R$ 23,438x106.
De acordo com o que foi apresentado na seção anterior, o valor do VPL5% da
UTE-Ótima otimizada pelo algoritmo da ED_mod foi de R$ 23,449x106. Desta forma,
pode ser afirmado que ambos os métodos de otimização proporcionaram praticamente o
mesmo resultado, pois a diferença relativa entre os dois resultados foi menor do que
0,05%.
Esta observação pode ser corroborada ao comparar os valores das variáveis de
148
otimização e os valores das variáveis de saída do Thermoflow TF1 , TF2 , TF3 e TF4
obtidos pelos dois métodos de otimização, disponibilizados nas Tabela 7.9 e Tabela
7.10.
Tabela 7.9 - Variáveis de Saída do Thermoflow: ED_mod vrs ED2SR
Variáveis de Otimização
Pressão de Vapor de Alta [bar]
Temperatura de Vapor de Alta [°C]
Pressão de Vapor de Reaquecimento [bar]
Temperatura de Vapor de Reaquecimento [°C]
Vazão Mássica de Vapor de Alta Pressão [t/h]
Vazão Mássica de Vapor para Reaquecimento [t/h]
UTE-Ótima
ED _mod ED2SR
100,0037 100,0129
483,5489 483,5814
24,9089 24,9207
482,9103 482,9051
398,2997 398,2929
390,1601 390,1607
Tabela 7.10 - Variáveis de Saída do Thermoflow: ED_mod vrs ED2SR
Variáveis de Sáida do Thermoflow
TF 1 - Potência Líquida [MW]
UTE-Ótima
ED _mod
ED2SR
476,0897
476,0853
6,5863
6,5864
TF 3 - Custo do Investimento [10 US$]
292,372
292,373
TF 4 - Potência Bruta [MW]
Índice de Custo de Instalação [US$/kW]
6
Função Objetivo - VPL 5% [10 R$]
490,9094
595,572
23,449
490,9051
595,580
23,438
6
TF 2 - Heat Rate [10 Btu/MWh]
6
De acordo com os resultados divulgados nas Tabela 7.9 e Tabela 7.10, é notório
que ambos os métodos de otimização proporcionaram praticamente os mesmos
resultados. Desta forma, a análise física e econômica entre a UTE-Base e UTE-Ótima
(otimizada pelo algoritmo da ED_mod) realizada na seção 7.3.1 é replicante ao projeto
ótimo obtido pelo algoritmo da ED2SR.
Outra forma de observar que os resultados podem ser considerados idênticos é
analisar a curva de probabilidade acumulada, construída a partir dos 10.000 cenários de
VPL deste projeto ótimo.
149
Figura 7.7 - VPL Ótimo pelo Método da ED2SR
Na Figura 7.7 estão identificados os pontos dos cenários mais otimista, mais
pessimista, de VPL igual à zero, de VPL médio e o valor da função objetivo do projeto
ótimo, o VPL5%. Os valores de cada um dos pontos presentes pondem ser vistos na
Tabela 7.11. Note que o projeto otimizado pelo método da ED2SR garante uma
probabilidade de ocorrer VPL menor do zero de 0,02%, mesma observação feita ao
projeto otimizado pelo método da ED_mod (terceira linha da Tabela 7.7).
Tabela 7.11 - VPL da UTE-Base pelo Método da ED2SR
Valor Presente Líquido da
UTE-ÓTIMA - ED2SR
VPL Otimista
VPL Pessimista
VPL igual a zero
VPL médio
VPL5% - Função Objetivo
6
VPL [10 R$]
70,261
-7,145
0
38,638
23,438
Probabilidade de Ocorrer
Valores Menores
100%
0%
0,02%
49,12%
5%
Na Tabela 7.12 a seguir é feito um paralelo entre os pontos evidenciados nas
150
curvas de probabilidade acumulada nos dois projetos ótimos (Figura 7.5 e Figura 7.7,
respectivamente).
É notório que os dois projetos geraram resultados com diferenças relativas muito
pequenas (menores do que 0,2%), levando a acreditar que os resultados podem ser
considerados idênticos sem nenhuma perda de fidelidade.
Tabela 7.12 - Pontos da Curva da UTE-Ótima - ED_mod vrs ED2SR
6
Valor Presente Líquido
UTE-Ótima [10 R$]
ED_mod
ED2SR
Diferença
Relativa
VPL Otimista
70,272
70,261
0,016%
VPL Pessimista
-7,134
-7,145
0,153%
VPL médio
38,649
38,638
0,029%
VPL5% - Função Objetivo
23,449
23,438
0,047%
Em relação à performance do algoritmo da ED2SR durante a otimização, de
acordo com a Figura 7.8 (linha tracejada), este algoritmo necessitou de
aproximadamente 5.000 cálculos de função objetivo para que a convergência fosse
atingida.
151
Figura 7.8 - Comparação das Convergências
Apesar da diferença entre as respostas ótimas obtidas nos dois métodos ser
praticamente imperceptível, o tempo desprendido para convergência do método da
ED2SR foi extremamente inferior, quando comparado com o método da ED_mod.
Esta redução do tempo computacional ocorreu, pois o método da ED2SR
precisou calcular apenas 5.000 vezes a função objetivo, enquanto o algoritmo da
ED_mod precisou calcular 30.000 funções objetivo.
Desta forma, o uso de uma função aproximada (modelada pelo método da
Superfície de Resposta) promoveu uma redução em uma ordem de grandeza no número
de simulações do Thermoflow e garantiu a fidelidade da resposta, quando acoplada ao
método da Evolução Diferenciada.
152
8
CONCLUSÃO
Foi visto no decorrer desta dissertação que investimentos em usinas
termelétricas no Brasil são classificados como de alto risco devido ao alto custo da
geração térmica, incerteza da taxa de câmbio, incerteza do preço do gás natural e
principalmente a baixa probabilidade de ocorrer uma hidrologia desfavorável.
Neste sentido foram encontradas na literatura várias propostas de investimento
em empreendimento termelétrico, segundo as regras impostas pelo Setor Elétrico
Brasileiro [1-3, 47-54] de forma a contornar os riscos inerentes. Entretanto, em todos
estes estudos as características técnicas das térmicas analisadas eram sempre
conhecidas.
Desta forma, a partir de um tipo de investimento flexível (ou seja, poder atuar
tanto no mercado cativo, quanto no livre) definido, o objetivo deste trabalho foi otimizar
os parâmetros do processo de geração de uma termelétrica em busca de um projeto
tecnicamente mais vantajoso, segundo as regras impostas pelo Setor Elétrico Brasileiro.
O estudo de caso apresentado no capítulo 7 mostrou que um empreendedor não
pode apenas levar em consideração o tipo de investimento. Além de escolher o melhor
investimento, segundo as regras impostas pelo Setor Elétrico Brasileiro, é necessário
sempre averiguar se o projeto que se deseja investir é realmente o melhor tecnicamente.
Este trabalho mostrou a grande vantagem em otimizar os parâmetros de uma
térmica, pois encontrou um projeto mais barato e tecnicamente mais simples, além de
praticamente eliminar o risco do seu Valor Presente Líquido ser negativo (menor do que
zero).
Para o cálculo de cada função objetivo durante o processo de otimização não
foram desenvolvidos modelos físicos e termodinâmicos de forma explícita para o ciclo
térmico. Neste trabalho foi utilizado o simulador de processos Thermoflow (software
153
comercial) com o propósito de resolver os balanços de massa e de energia, para que o
Valor Presente Líquido fosse então calculado. Com isso, por utilizar um simulador de
processo o problema de otimização demostrou ser computacionalmente demorado, pois
cada simulação do Thermoflow requereu de 4 a 5 segundos.
Por este motivo, este trabalho propôs com sucesso o desenvolvimento e a
aplicação de um método de Superfície de Resposta em conjunto com o algoritmo de
otimização para auxiliar a redução do tempo computacional citado, a partir da redução
do número total de funções objetivos necessárias durante o processo de otimização.
Sobre a técnica de otimização empregada, conforme apresentado nos capítulos 2
e 6, foi utilizado um método de busca global e heurístico, a Evolução Diferenciada. Na
seção 6.4 do capítulo 6 foi demonstrada a eficácia deste método quando aplicada a
funções teste. Nesta mesma seção foi empregado o modelo de Superfície de Resposta
em conjunto com a Evolução Diferenciada, quando da otimização das mesmas funções
teste. Quando o método de otimização fez uso da Superfície de Resposta o número de
cálculo de função objetivo necessários para encontrar o ponto ótimo diminuiu
significativamente.
Este dois métodos de otimização (um sem Superfície de Resposta e o outro com
Superfície de Resposta) foram empregados na função objetivo proposta neste trabalho, o
Valor Presente Líquido de uma térmica segundo as regras impostas pelo Setor Elétrico
Brasileiro, associada a um risco de apenas 5%.
Ambos os métodos convergiram para o mesmo resultado. Entretanto o algoritmo
proposto, que utilizou o modelo de Superfície de Resposta, reduziu em uma ordem de
grandeza a quantidade de funções objetivo necessárias para achar o melhor projeto.
Como propostas, novos trabalhos poderiam testar o emprego da Superfície de
Resposta em problemas reais de engenharia em diferentes áreas, uma vez que este
154
método mostrou excelentes resultados quando testado em funções de teste e no
problema real proposto.
Os algoritmos empregados neste trabalho não foram testados para problemas
com variáveis de otimização maiores do que 10. Assim, trabalhos futuros poderiam
testar o emprego de algoritmos de otimização com Superfície de Resposta para otimizar
problemas desta magnitude e verificar a sua performance.
155
REFERÊNCIAS BIBLIOGRÁFICAS
[1]
CAVALCANTE, A. F., JUNIOR, J. L. T., 2007, “Avaliação de Investimentos de
Capital na Geração Termoelétrica Usando a Teoria das Opções Reais: Um
Estudo de Caso Utilizando a Equação de Bellman”. VII Encontro Brasileiro de
Finanças, IBMEC, São Paulo, SP, Brasil, 26-27, jul. 2007.
[2]
CASTRO, A. L., 2000, Avaliação de Investimento de Capital em Projetos de
Geração Termelétrica no Setor Elétrico Brasileiro Usando a Teoria de Opções
Reais. Dissertação de M.Sc., PUC, Rio de Janeiro, RJ, Brasil.
[3]
OLIVEIRA, B. N., 2008, Modelo de Comercialização de Energia pela Opção de
Disponibilidade na Geração Termelétrica. Dissertação de M.Sc., UFRJ, Rio de
Janeiro, RJ, Brasil.
[4]
STORN, R., PRICE, K., 1995, Differential Evolution – a simple and efficient
adaptive scheme for global optimization over continuous spaces. In: Technical
Report, International Computer Science Institute, Berkley.
[5]
PIRES, T. S., 2010, Método de Superfície de Resposta Aplicado à Otimização
Termoeconômica de Sistemas de Cogeração Modelados em um Simulador de
Processos. Dissertação de M.Sc., COPPE/UFRJ, Rio de Janeiro, RJ, Brasil.
[6]
CHENG, AH-D., 2012, “Multiquadric and its shape parameter - A numerical
investigation of error estimate, condition number, and round-off error by
arbitrary precision computation”, Engineering Analysis with Boundary Elements,
v.36, pp. 220–239.
[7]
STORN, R., PRICE, K., 1997, “Differential evolution – A simple and efficient
heuristic for global optimization over continuous spaces”. Journal of Global
Optimization., v. 11, pp. 341-359.
[8]
CCEE – CÂMARA DE COMERCIALIZAÇÃO DE ENERGIA ELÉTRICA,
2012. Disponível em: <http://www.ccee.org.br>. Acesso em 20 de fev. 2012.
[9]
TOLMASQUIM, M. T., SZKLO, A. S., SOARES, J. B., 2002, Análise da
Viabilidade de Introdução de Gás Natural em Setores Selecionados. Relatório
Técnico. Convênio FINEP/CT-Petro Rio de Janeiro.
[10]
COMERCIALIZADORA ENERGISA, 2012. Disponível em:
<http://www.grupoenergisa.com.br/Comercializadora/Mercado%20Livre/Dicion
Dicionariod.aspx>. Acesso: 27 de mar. 2013.
[11]
BREST, J., GREINER, S., BOSKOVIC B., MERNIK, M., ZUMER, V., 2006,
“Self-adapting control parameters in differential evolution: a comparative study
on numerical benchmark problems”. IEEE Transactions on Evolutionary
Computation, v. 10, n. 6, pp. 646-657.
[12]
PANT, M., THANGARAJ, R., SINGH, V. P., 2009, “Optimization of
156
mechanical design problems using improved differential evolution algorithm”.
International Journal of Recent Trends in Engineering, v. 1, pp. 21-25.
[13]
VESTERSTROM, J., THOMSEN, R., 2004, “A Comparative Study of
Differential Evolution, Particle Swarm Optimization, and Evolutionary
Algorithms on Numerical Benchmark Problems”. Evolutionary Computation, v.
2, pp. 1980-1987.
[14]
HOLLAND, J. H., 1975, “Adaptation in Natural and Artificial System: An
Introductory Analysis with Applications to Biology, Control, and Artificial
Intelligence”. University of Michigan Press.
[15]
EIBEN, A. E., SMITH, J. E., 2003, Introduction to evolutionary computation.
Springer, Berlin.
[16]
KENNEDY, J., EBERHART, R. C., 1995, “Particle Swarm Optimization”.
Proceedings of the IEEE International Conference on Neural Networks., v. 4,
pp. 1942-1948.
[17]
DORIGO, M., STÜTZLE, T., 2004, “Ant Colony Optimization”. MIT Press,
Cambridge, MA.
[18]
KARABOGA, D., OKDEM, S., 2004, “A simple and Global Optimization
Algorithm for Engineering Problems: Differential Evolution Algorithm”. Turk J.
Elec. Engin., v. 12, n. 1, pp. 53-60.
[19]
STORN, R., PRICE, K., 1996, “Minimizing the real functions of the ICEC 1996
contest by differential evolution”. In: IEEE Int. Conf. Evol. Comput., pp. 842844.
[20]
PRICE, K. V., 1997, “Differential evolution vs. the functions of the 2nd ICEO,”
in Proc. IEEE Int. Conf. Evol. Comput, (Abril), pp. 153–157.
[21]
STORN, R., 1999, “System design by constraint adaptation and differential
evolution”, IEEE Transactions on Evolutionary Computation, v.3, n. 1, pp. 2234.
[22]
LAMPINEN J, ZELINKA I, 2000, “On stagnation of the differential evolution
algorithm”. In: Osmera P (ed) Proceedings of 6th international mendel
conference on soft computing, pp. 76-83.
[23]
QIN, A. K., SUGANTHAN, P. N., 2005, “Self-adaptive Differential Evolution
Algorithm for Numerical Optimization”. Proc. of the 2005 Congress on
Evolutionary Computation, v. 2, pp. 1785-1791.
[24]
HUANG, V. L., QIN, A. K., SUGANTHAN, P. N., 2006, “Self-adaptive
differential evolution algorithm for constrained real-parameter optimization”.
Proc. IEEE Congress on Evolutionary Computation, pp. 17-24.
[25]
NERI, F., TIRRONEN, V., 2008, “On memetic differential evolution
frameworks: a study of advantages and limitations in hybridization”. In:
Proceedings of the IEEE world congress on computational intelligence, pp.
157
2135-2142.
[26]
LIU J., LAMPINEN J., 2002a, “On setting the control parameter of the
differential evolution algorithm”. In: Proceedings of the 8th international
Mendel conference on soft computing, pp. 11-18
[27]
LIU J., LAMPINEN J., 2002b, “A fuzzy adaptive differential evolution
algorithm”. In: Proceedings of the 17th IEEE region 10th international
conference on computer, communications, control and power engineering, vol.
1, pp. 606–611.
[28]
ZIELINSKI K.,WEITKEMPER P., LAUR R., KAMMEYER K-D., 2006,
“Parameter study for differential evolution using a power allocation problem
including interference cancellation”. In: Proceedings of the IEEE congress on
evolutionary computation, pp. 1857-1864.
[29]
ZHANG, J., SANDERSON, A. C., 2009, “Adaptive differential evolution-A
robust approach to multimodel problem optimization”. Series of Adaptation,
Learning, and Optimization, New York: Springer-Verlag, (ago.).
[30]
NORMAN, N., IBA, H., 2008, “Accelerating Differential Evolution Using an
Adaptive Local Search”. IEEE Transaction on Evolutionary Computation, v.12,
n. 1 (feb.), pp. 107-125.
[31]
WANG, Y., CAI, Z., ZHANG, Q., 2011, “Differential evolution with composite
trial vector generation strategies and control parameters”. IEEE Trans. on Evol.
Comput, v. 15, n. 1, pp. 55–66.
[32]
QIN, A. K., HUANG, V. L., SUGANTHAN, P. N., 2009, “Differential
evolution algorithm with strategy adaptation for global numerical optimization”,
IEEE Trans. On Evolutionary Computation, v.13, n. 2 (abril), pp. 398-417.
[33]
MALLIPEDDI, R., SUGANTHAN, P. N., PAN, Q. K., TASGETIREN, M. F.,
2011, “Differential Evolution Algorithm with ensemble of parameters and
mutation strategies,” Applied Soft Computing, v. 11, n. 2 (mar.), pp. 1679-1696.
[34]
ISLAM, SK. M., DAS, S., GHOSH, S., ROY, S., SUGANTHAN, P. N., 2012,
“An Adaptive Differential Evolution Algorithm with Novel Mutation and
Crossover Strategies for Global Numerical Optimization”, IEEE Trans. on SMCB, v. 42, n. 2, pp. 482-500.
[35]
ZHANG, J., SANDERSON, A. C., 2009, “JADE: Adaptive Differential
Evolution with Optional External Archive”, IEEE Trans. Evolut. Comput., v. 13,
n. 5, pp. 945-958.
[36]
YAN, Y., GUO, C., GONG, W., 2012, “Improved Adaptive Differential
Evolution based on Weighted Mean”. Journal of Comput. Information Systems,
v. 8, n. 9, pp. 3723-3737.
[37]
GONG, W., FIALHO, A., HUI LI, Z. C., 2011, “Adaptive strategy selection in
differential evolution for numerical optimization: An empirical study”. Inf. Sci.
v. 181, n. 24, pp. 5364-5386.
158
[38]
DAS, S., SUGANTHAN, P. N., 2011, “Differential Evolution: A Survey of the
State-of-the-art”. IEEE Transactions on Evolutionary Computation, v. 15, n. 1
(Fev.), pp. 4-31.
[39]
STORN, R., 1999, “System Design by Constraint Adaptation and Differential
Evolution”. IEEE Transactions on Evolutionary Computation, v. 3, n. 1, pp. 2234.
[40]
CHENG, S. L.; HWANG, C., 2001, “Optimal Approximation of Linear Systems
by a Differential Evolution Algorithm”. IEEE Transactions on Systems, Man,
and Cybernetics - Part A: Systems and Humans, v. 31, n. 6, pp. 698–707.
[41]
BABU, B. V.; MUNAWAR, S. A., 2001, “Optimal Design of Shell and Tube
Heat Exchanger by Different strategies of Differential Evolution”.
PreJournal.com - The Faculty Lounge.
[42]
BERGAMASCHI, P.R.; SARAMAGO, S.F.P., COELHO, L.S, 2005, “Evolução
Diferencial aplicada à otimização do volume do espaço de trabalho de um
manipulador”. VII SBAI / II IEEE LARDS. Sao Luis.
[43]
MARIANI, V. C., COELHO, L. S., SAHOO, P. K., 2011, “Modified differential
evolution approaches applied in exergoeconomic analysis and optimization of a
cogeneration system”. Expert Systems with Applications, v. 38, n. 11 (Out.), pp.
13886-13893.
[44]
SAYYAADI H., 2009, “Multi-objective approach in thermoenvironomic
optimization of a benchmark cogeneration system”, Applied Energy, v. 86, n. 6,
(jun.), pp. 867-879.
[45]
TSATSARONIS, G., 2007, “Optimization of combined cycle power plants using
evolutionary algorithms”. Chemical Engineering and Processing, v. 46, pp.
1151-1159.
[46]
BURATTI, R. M., 2008, Estratégia de contratação de Energia Elétrica para
uma concessionária de distribuição. Dissertação de M.Sc., PUC-PR, Curitiba,
PR, Brasil.
[47]
SOARES, B. L., 2008, Seleção de Projetos de Investimento em Geração de
Energia Elétrica. Dissertação de M.Sc., PUC/RIO, Rio de Janeiro, RJ, Brasil.
[48]
FODRA, M., 2012, Viabilidade Econômica da Venda de Energia Elétrica em
Co-Geração Sob Condições de Risco: Um Estudo de Caso. Tese de D.Sc.,
UNESP, Botucatu, SP, Brasil.
[49]
SZCZERBACKI, C. F., 2007, Formação de Preços de Energia Elétrica para o
Mercado Brasileiro. Dissertação de M.Sc., PUC, Rio de Janeiro, RJ, Brasil.
[50]
FONTOURA, C. F. V. T., 2011, Avaliação de Projeto de Investimento em Usina
Termelétrica à Capim-Elefante: Uma Abordagem Pela Teoria de Opções Reais.
Dissertação de M.Sc., PUC, Rio de Janeiro, RJ, Brasil.
159
[51]
NASCIMENTO, W. J. D., 2008, Conversão de Termelétricas para BiCombustível em Ambiente de Incerteza: Uma Abordagem por Opções Reais.
Dissertação de M.Sc., PUC, Rio de Janeiro, RJ, Brasil.
[52]
TATONI, W. M, 2012, Avaliação de Projetos de Investimento em Cogeração de
Energia Utilizando Bagaço de Cana-de-Açúcar em Biorrefinarias a Partir do
Uso da Teoria das Opções Reais. Dissertação de M.Sc., FGV-EESP, São Paulo,
SP, Brasil.
[53]
JUNIOR, B. S, 2012, Avaliação da Atratividade de Negócios em Geração
Distribuída e Economia de Energia Elétrica Piloto Aplicado Dentro dos Estudos
do PIR da RAA. Dissertação de M.Sc., USP, São Paulo, SP, Brasil.
[54]
MARTINS, D. M. R., 2008, Setor Elétrico Brasileiro: Análise do Investimento
de Capital em Usinas Termelétricas. Dissertação de M.Sc., PUC, Rio de Janeiro,
RJ, Brasil.
[55]
COLAÇO M. J. SILVA DULIKRAVICH G. S. SAHOO D., 2008a, “A response
surface method-based hybrid optimizer”. Inverse Problems in Science and
Engineering, v. 16, n. 6 (Set.), pp. 717-741.
[56]
COLAÇO M. J. ORLANDE H. R. B. DULIKRAVICH G. S., 2006, “Inverse
and Optimization Problems in Heat Transfer”. Journal of the Brazilian Society of
Mechanical Sciences and Engineering, v. 28, n. 1 (Mar.), pp. 1-24.
[57]
COLAÇO, M. J., SILVA, W. B., MAGALHÃES, A. C., DULIKRAVICH, G.
S., 2008b, “Response Surface Methods Applied to Scarce and Small Sets of
Training Points – A Comparative Study”. In: Proceedings of EngOpt 2008 –
International Conference on Engineering Optimization, Rio de Janeiro, Jun.
[58]
LEITÃO V. M. A., 2001, “A Meshless Method for Kirchhoff Plate Bending
Problems”. International Journal of Numerical Methods in Engineering, v. 52,
n. 10 (Dez.), pp. 1107-1130.
[59]
LEITÃO V. M. A., 2004, “RBF-based meshless methods for D elastostatic
problems”. Engineering Analysis with Boundary Elements, v. 28, n. 10 (Out.),
pp. 1271-1281.
[60]
ICB - ÍNDICE DE CUSTO BENEFÍCIO DE EMPREENDIMENTOS DE
GERAÇÃO TERMELÉTRICA - METODOLOGIA DE CÁLCULO, 2011.
Disponível em: <http://www.aneel.org.br>. Acesso em 20 de fev. 2012.
[61]
CMO – CUSTO MARGINAL DE OPERAÇÃO, 2012. Disponível em: <
http://www.epe.gov.br>. Acesso em 02 de fev. de 2013.
[62]
EPE - EMPRESA DE PESQUISA ENERGÉTICA, 2013. Disponível em:
<www.epe.gov.br>. Acesso em 23 de jan. 2013.
[63]
NOTA TÉCNICA - DEE/DPG - RE- 001/2009-r1. “Projeção dos Preços dos
Combustíveis para Determinação do CVU das Termelétricas para Cálculo da
Garantia Física e dos Custos Variáveis da Geração Termelétrica (COP e CEC)”.
Disponível em: <www.epe.gov.br>. Acesso em 23 de jan. 2013.
160
[64]
EIA - ENERGY INFORMATION ADMINISTRATION, 2013. Disponível em:
<www.eia.gov>. Acesso em 19 de mar. de 2013.
[65]
ODDONE, D. C., 2001, Cogeração: Uma Alternativa para Produção de
Eletrecidade. Dissertação de M.Sc., USP, São Paulo, SP, Brasil.
[66]
BRANCO, F. P., 2005, Analise Termoeconomica de uma UTE a GN - Ciclo
Aberto e Combinado. Dissertação de M.Sc., UNESP, Ilha Solteira, SP, Brasil.
[67]
PORTARIA Nº 108, DE 8 DE MARÇO DE 2012. Disponível em:
<www.ccee.gov.br>. Acesso em 12 de ago. 2013.
[68]
PORTARIA Nº 169, DE 22 DE MARÇO DE 2012. Disponível em:
<www.ccee.gov.br>. Acesso em 12 de ago. 2013.
[69]
ELETROBRAS, 2013. Disponível em: <http://www.eletrobras.org.br>. Acesso
em 13 de ago. de 2013.
[70]
MPX Energia, 2012. Disponível em: <http://www.mpx.com.br>. Acesso em 07
de ago. de 2013.
[71]
COSTA, A. N., 2008, Otimização da Lucratividade de Plantas de Cogeração:
Modelagem do Problema PCLM. Dissertação de M.Sc., COPPE/UFRJ, Rio de
Janeiro, RJ, Brasil.
[72]
ASSIS, R., 2010, Apoio à Decisão em Manutenção na Gestão de Activos
Físicos. 7 ed. Coimbra, Lidel-Zamboni.
[73]
SOUZA, G. C. U. I., 2006, Avaliação de Títulos Conversíveis com Opções de
Compra e Venda Implícitas em Contrato. Tese de D.Sc., PUC, Rio de Janeiro,
RJ, Brasil.
[74]
VIOLA, M. L. L., 2006, Teoria De Valores Extremos E Cópulas: Distribuição E
Valor Extremo E Cópulas Arquimedianas E Generalizadas Trivariadas.
Dissertação de M.Sc., IMECC-UNICAMP, Campinas, SP, Brasil.
[75]
STORN, R., PRICE, K., 2005, “Differential Evolution: A Practical Approach to
Global Optimization”. Springer-Verlag, Berlin.
[76]
ONG, Y.-S., KEANE, A. J., 2004, “Meta-Lamarckian learning in memetic
algorithms”. IEEE Trans. on Evol. Comput, v. 8, n. 2 (Abril), pp. 99-110.
[77]
NOTA TÉCNICA - EPE-DEE-IT-56/2013 “Preços Médios de Referência dos
Combustíveis Vinculados ao CVU de Usinas Termelétricas”. Disponível em:
<http://www.epe.gov.br>. Acesso em 08 ago. de 2013.
[78]
BANCO
CENTRAL
DO
BRASIL,
2013.
<http://www.bcb.gov.br>. Acesso em 20 ago. de 2013.
[79]
CMO – CUSTO MARGINAL DE OPERAÇÃO, 2013. Disponível em:
<http://www.epe.gov.br>. Acesso em 08 ago. de 2013.
Disponível
em:
161
[80]
Receita
Federal
do
Brasil,
2013.
Disponível
<http://www.receita.fazenda.gov.br>. Acesso em 27 jul. de 2013.
em:
[81]
NEMA Standards Publication N° 8M 23-1991, R1997, R2002, Steam Turbines
For Mechanical Drive Service. Rosslyn, IHS sob licença da NEMA - National
Electrical Manufacturers Association.