Manual Referencia Dessem GTRD - simsee.org · Ressalta-se que cada um destes modelos necessita, para sua utilização, de uma série de modelos auxiliares para fornecer os dados necessários

CEPEL

Manual de Referência doModelo DESSEM

Versão 8.2a (GTRD)

Julho de 2003

CEPEL Modelo DESSEM ii

Manual de referência – Versão 8.2a (GTRD) Julho/2003

APRESENTAÇÃO

Este relatório tem o objetivo de servir como referência para o conhecimento da estratégia desolução e funcionalidades contempladas atualmente pelo modelo DESSEM. Estas serãoobjeto de análise pelo Grupo de Trabalho para Reformulação do DESSEM (GTRD), para aelaboração de um termo de referência com as modificações necessárias para futura utilizaçãodo modelo como ferramenta de apoio à realização da programação diária pelo ONS.

Procurou-se dar enfoque, de forma paralela, tanto do aspecto metodológico dasfuncionalidades do modelo como de sua formulação matemática, tal como implementada emsua versão atual (8.2a), disponibilizada pelo CEPEL ao ONS em Junho de 2003. A fim defacilitar a análise, foram omitidos alguns desenvolvimentos que, embora já implementados econstantes dessa versão, não serão considerados pela GTRD, por estarem além do escopopreviamente definido pelo Submódulo 8.1 para a realização da programação diária. Destacam-se neste aspecto o tratamento das restrições de unit commitment das unidades térmicas erestrições que considerem individualmente as unidades hidroelétricas, como as zonasproibidas de geração.

O documento está organizado da seguinte forma: o item introdutório apresenta o eixoprincipal da cadeia de modelos desenvolvida pelo CEPEL para o planejamento e operação dosistema elétrico brasileiro. No item 2, descreve-se de forma geral o problema de planejamentoda operação energética e no item 3 aborda-se especificamente o problema da programaçãodiária, onde descreve-se sucintamente suas características e formulação matemática eapresenta-se a estratégia de solução adotada atualmente no modelo DESSEM. O item 4 listaos tipos de funcionalidades atualmente existentes no modelo, a serem analisadas com o intuitode atender às recomendações do submódulo 8.1. Nos itens 5, 6 e 7 são detalhadas asfuncionalidades presentes em cada tipo. Finalmente, no item 8, mostra a formulação doproblema de fluxo em redes DC considerado pelo DESSEM nos casos com rede elétrica e aestratégia integrada de solução do problema na otimização do despacho como um todo.

Informações pertinentes à utilização em si do modelo, tais como arquivos e formatos dosdados de entrada e saída de resultados, podem ser obtidas por consulta ao Manual do Usuáriodo modelo DESSEM.

dino

Note

Restrições de unit commitment de térmicas não são necessárias num sistema hidrotérmico como o SIN, onde o acompanhamento da carga pode ser feito exclusivamente com geração hidrelétrica.

CEPEL Modelo DESSEM iii


ÍNDICE

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

2 O PROBLEMA DE PLANEJAMENTO DA OPERAÇÃO ................................ 32.1 SISTEMAS PURAMENTE TÉRMICOS .............................................................................. 32.2 SISTEMAS HIDROTÉRMICOS......................................................................................... 32.3 SISTEMAS INTERLIGADOS............................................................................................ 5

3 PROGRAMAÇÃO DIÁRIA DA OPERAÇÃO - MODELO DESSEM ............. 63.1 REPRESENTAÇÃO TEMPORAL ...................................................................................... 63.2 REPRESENTAÇÃO DO SISTEMA..................................................................................... 73.3 FORMULAÇÃO DO PROBLEMA ..................................................................................... 73.4 ESTRATÉGIA DE SOLUÇÃO ......................................................................................... 10

4 TIPOS DE FUNCIONALIDADES DO MODELO DESSEM........................... 12

5 FUNCIONALIDADES ENVOLVENDO CARACTERÍSTICAS /RESTRIÇÕES DO SISTEMA .............................................................................................. 145.1 DESCRIÇÃO ............................................................................................................... 14

5.1.1 Carga por subsistema (estudos sem rede elétrica)............................................... 145.1.2 Intercâmbio entre subsistemas ............................................................................. 145.1.3 Custo de déficit por profundidade de corte de carga........................................... 145.1.4 Restrição de transporte entre Itaipu 50 e 60 Hz e os subsistemas SU/SE ........... 155.1.5 Contratos de importação/exportação de energia................................................. 155.1.6 Consideração de gerações fixas (gerações não otimizadas ou simuladas) ......... 155.1.7 Consumo nas unidades elevatórias ...................................................................... 155.1.8 Restrições elétricas .............................................................................................. 15

5.2 MODELAGEM............................................................................................................. 165.2.1 Termos na função objetivo ................................................................................... 165.2.2 Equação de Atendimento a demanda ................................................................... 165.2.3 Limites de intercâmbio ......................................................................................... 165.2.4 Restrição de transporte entre Itaipu 50 e 60 Hz e os subsistemas SU/SE ........... 165.2.5 Limites de contrato de importação/exportação de energia.................................. 175.2.6 Restrições elétricas .............................................................................................. 17

6 FUNCIONALIDADES ENVOLVENDO CARACTERÍSTICAS /RESTRIÇÕES DAS USINAS HIDROELÉTRICAS.......................................................... 186.1 DESCRIÇÃO DAS FUNCIONALIDADES REFERENTES AO BALANÇO HÍDRICO NO SISTEMA18

6.1.1 Afluências naturais determinísticas aos aproveitamentos ................................... 186.1.2 Evaporação .......................................................................................................... 186.1.3 Retiradas de água para outros usos..................................................................... 186.1.4 Limites físicos de armazenamento nos reservatórios........................................... 186.1.5 Enchimento de volume morto............................................................................... 186.1.6 Tempo de viagem da água entre aproveitamentos ............................................... 196.1.7 Representação de unidades elevatórias ............................................................... 196.1.8 Vertimento máximo nos reservatórios .................................................................. 20

6.2 DESCRIÇÃO DAS FUNCIONALIDADES REFERENTES ÀS RESTRIÇÕES OPERATIVAS ...... 20

CEPEL Modelo DESSEM iv


6.2.1 Limites de armazenamento................................................................................... 206.2.2 Volume de espera para controle de cheias ........................................................... 206.2.3 Limites de turbinamento, vertimento, desvio e defluência total........................... 206.2.4 Limites de bombeamento nas unidades elevatórias ............................................. 20

6.3 DESCRIÇÃO DAS FUNCIONALIDADES REFERENTES À GERAÇÃO DAS USINAS ............ 206.3.1 Função de Produção Energética.......................................................................... 206.3.2 Reserva de potência ............................................................................................. 216.3.3 Manutenção programada ..................................................................................... 216.3.4 Variação da vazão turbinada entre estágios consecutivos................................... 226.3.5 Variação da geração entre estágios consecutivos................................................ 22

6.4 DESCRIÇÃO DAS FUNCIONALIDADES REFERENTES ÀS PARTICULARIDADES DOSISTEMA ................................................................................................................... 22

6.4.1 Variação diária e horária na cota da Régua 11................................................... 226.4.2 Operação do canal Pereira Barreto..................................................................... 226.4.3 Operação das bacias do Alto Tietê e Paraíba do Sul .......................................... 23

6.5 MODELAGEM............................................................................................................. 236.5.1 Termos na função objetivo ................................................................................... 236.5.2 Equação de balanço hídrico ................................................................................ 236.5.3 Restrição de vertimento máximo .......................................................................... 246.5.4 Volume evaporado ................................................................................................ 246.5.5 Limites de armazenamento físico e operacional .................................................. 246.5.6 Restrição de volume de espera............................................................................. 246.5.7 Limites de turbinamento, vertimento, desvio de água por canal e defluência total246.5.8 Limites de bombeamento...................................................................................... 246.5.9 Limites de geração ............................................................................................... 246.5.10 Restrição de reserva de potência ......................................................................... 246.5.11 Função de produção energética ........................................................................... 256.5.12 Restrições de variação da vazão turbinada ......................................................... 256.5.13 Restrições de variação horária e diária no nível da Régua 11............................ 256.5.14 Limite da vazão no canal Pereira Barreto ........................................................... 25

7 FUNCIONALIDADES DAS USINAS TÉRMICAS........................................... 267.1 DESCRIÇÃO ............................................................................................................... 26

7.1.1 Rampa para tomada/alívio de carga.................................................................... 267.1.2 Reserva de potência ............................................................................................. 267.1.3 Manutenção programada ..................................................................................... 267.1.4 Custos e limites físicos de geração por unidade .................................................. 267.1.5 Limites operativos de geração para usinas e unidades ....................................... 26

7.2 MODELAGEM............................................................................................................. 267.2.1 Termos na função objetivo ................................................................................... 267.2.2 Limites de geração das usinas e unidades térmicas ............................................ 277.2.3 Restrição de reserva de potência ......................................................................... 27

8 CONSIDERAÇÃO DA REDE ELÉTRICA NO CÁLCULO DO DESPACHOPELO MODELO DESSEM .................................................................................................. 28

8.1 FORMULAÇÃO DC DA REDE ELÉTRICA..................................................................... 288.2 DESCRIÇÃO DA METODOLOGIA ADOTADA NO DESSEM .......................................... 288.3 OBSERVAÇÕES ........................................................................................................... 30

9 REFERÊNCIAS..................................................................................................... 32

CEPEL Modelo DESSEM v


ANEXO I: NOTAÇÕES UTILIZADAS............................................................................... 33

ANEXO II - PROGRAMAÇÃO DINÂMICA DUAL DETERMINÍSTICA (PDDD)..... 38

II.1 Algoritmo de solução............................................................................................... 40II.2 Consideração do tempo de viagem.......................................................................... 43II.2.1 Exemplo com tempo de viagem .................................................................................. 46

ANEXO III: FORMULAÇÃO MATEMÁTICA COMPLETA DO PROBLEMA DEDESPACHO RESOLVIDO PELO MODELO DESSEM................................................... 48

CEPEL Modelo DESSEM 1


1 INTRODUÇÃO

O planejamento da operação de sistemas hidrotérmicos de geração de energia elétrica abrangeum largo espectro de atividades, desde o planejamento plurianual até a programação diária daoperação dos reservatórios. Devido ao porte e complexidade deste problema, é necessária asua divisão em diversas etapas. Em cada etapa, são utilizados modelos com diferentes grausde detalhe para representação do sistema, abrangendo períodos de estudos com horizontesdistintos, de médio prazo, curto prazo e programação diária. Para a determinação da operaçãodo sistema integrado nacional (SIN) brasileiro, o Operador Nacional do Sistema (ONS)dispõe de uma cadeia de modelos desenvolvidos pelo CEPEL [1].

Na etapa de médio prazo, o ONS utiliza o modelo NEWAVE [2], realizando estudos mensaiscom horizontes de até 5 anos, objetivando definir parcelas de geração hidro e termoelétricaque minimizem o valor esperado do custo de operação, representado o sistema através de ummodelo equivalente. No curto prazo, determinam-se, a cada semana, metas individuais degeração para as usinas térmicas e hidroelétricas, bem como os intercâmbios de energia entresubsistemas, a partir das informações da etapa de médio prazo. Para esta etapa, utiliza-se omodelo DECOMP [3]. Finalmente, na programação diária, define-se, a partir das informaçõesdo curto prazo, uma programação de geração sujeita às condições operativas da rede elétrica.

Para esta etapa de programação diária, foi concebido o modelo DESSEM [4]. Seu objetivo édeterminar metas de geração hidro ou termoelétrica, por usina1, que minimizem o custo totalde operação do sistema, composto por uma parcela de custo presente, ao longo do horizontede estudo (geração térmica, contratos de energia com sistemas externos, déficit noatendimento) e outra de custo futuro, cuja função é fornecida pelo modelo DECOMP erepresenta o valor esperado dos custos de geração térmica e déficit de energia no futuro.

A seguir apresenta-se um diagrama esquemático contendo o horizonte e o intervalo dediscretização de cada um dos modelos disponibilizados pelo CEPEL para o planejamento eprogramação da operação do SIN.

Horizonte de estudo Intervalo de Discretização

Figura 1.1 – Cadeia de modelos desenvolvidos pelo CEPELpara o planejamento e programação da operação do SIN

1 Alternativamente, o modelo permite uma modelagem detalhada por unidade geradora hidro ou termoelétrica

NEWAVE

DECOMP

DESSEM

Função deCusto Futuro

Função deCusto Futuro

MÉDIO PRAZO(até 5 anos)

CURTO PRAZO(até 1 ano)

PROGRAMAÇÃO DIÁRIA(até 14 dias)

mensal

Semanal / mensal

½ hora, 1 hora ou superior(patamares cronológicos)



Ressalta-se que cada um destes modelos necessita, para sua utilização, de uma série demodelos auxiliares para fornecer os dados necessários para a determinação do planejamento,tais como previsão de carga, geração de cenários ou previsão de afluências, cronograma demanutenção das usinas, entre outros.



2 O PROBLEMA DE PLANEJAMENTO DA OPERAÇÃO

O objetivo básico do planejamento da operação é obter, para cada etapa, as metas de geraçãode cada usina (hidro e termoelétrica) do sistema de forma a atender a demanda e minimizar ovalor esperado do custo de operação ao longo do período de planejamento. Este custo écomposto pelo custo variável de combustível das usinas termoelétricas e pelo custo atribuídoàs interrupções de fornecimento de energia, representado por uma função de penalização dosdéficits de energia (custo do déficit).

2.1 Sistemas Puramente Térmicos

Em sistemas de geração compostos somente de unidades termoelétricas, o custo de cada usinadepende basicamente do custo de combustível. Portanto, o problema de operação consiste emdeterminar a combinação de usinas que minimize o custo total de combustível (carvão, óleo,nucleares, etc.) necessário para atender a demanda, satisfazendo os limites de geração emcada usina térmica.

Em sua versão mais simples, este problema se resolve colocando as usinas em ordemcrescente do custo de produzir um MWh adicional (o custo incremental) e ajustando aoperação às flutuações da demanda. O último gerador alocado é chamado de geradormarginal, pois ele atenderá a qualquer acréscimo extra de carga que ocorra. Todos osgeradores com custo inferior ao marginal operam em sua capacidade máxima e os geradorescom custo de operação superior ao do marginal não entram em operação. Ainda que existamfatores adicionais que tornam este problema mais complexo (perdas de energia, limitações naslinhas de transmissão, custos de partida, limites na taxa de variação da produção energética,etc.), o problema de operação termoelétrico tem características básicas, resumidas a seguir:

a. é desacoplado no tempo, quer dizer, uma decisão de operação hoje não tem efeito nocusto de operação do próximo estágio;

b. as unidades têm um custo direto de operação, quer dizer, o custo de operação de umaunidade não depende do nível de geração de outras unidades. Além disso, a operaçãode uma unidade não afeta a capacidade de geração ou a disponibilidade de outraunidade;

c. a confiabilidade do fornecimento de energia depende somente da capacidade total degeração disponível e não da estratégia de operação das unidades do sistema.

2.2 Sistemas Hidrotérmicos

Ao contrário dos sistemas puramente termoelétricos, sistemas com uma porcentagemsubstancial de geração hidroelétrica podem utilizar a energia armazenada "grátis" nosreservatórios do sistema para atender a demanda, substituindo assim a geração dispendiosadas unidades termoelétricas.

Entretanto, o volume de água afluente aos reservatórios é desconhecido, pois dependebasicamente das chuvas que irão ocorrer no futuro. Além disso, a disponibilidade de energiahidroelétrica é limitada pela capacidade de armazenamento nos reservatórios. Isto introduzuma relação entre uma decisão de operação em uma determinada etapa e as conseqüênciasfuturas desta decisão. Por exemplo, se decidirmos utilizar as provisões de energiahidroelétrica para atender o mercado hoje e no futuro ocorrer uma seca, pode vir a sernecessário utilizar geração termoelétrica de custo elevado ou, mesmo, interromper ofornecimento de energia. Por outro lado, se optarmos por fazer uso mais intensivo de geraçãotermoelétrica, conservando elevados os níveis dos reservatórios e ocorrerem vazões altas no



futuro, poderá haver vertimento no sistema, o que representa um desperdício de energia e, emconseqüência, um aumento desnecessário no custo de operação. Esta situação está ilustrada naFigura 2.1.

minimizar o custo de

combustível esvaziando os reservatórios

manter os reservatórios cheios e usar

geração termoelétrica

baixas

baixas

altas

altas

déficit

vertimento

Conseqüências operativas

Afluências Futuras Decisão

operação econômica

operação econômica

Figura 2.1 - Processo de decisão para sistemas hidrotérmicos

O problema de planejamento da operação de sistemas hidrotérmicos possui características quepodem ser assim resumidas como:

• é acoplado no tempo, ou seja, é necessário avaliar as conseqüências futuras de umadecisão no presente. A solução ótima é um equilíbrio entre o benefício presente do usodo água e o beneficio futuro de seu armazenamento, medido em termos da economiaesperada dos combustíveis das unidades térmicas.

• a este problema dinâmico se agrega o problema da irregularidade das vazões afluentesaos reservatórios, que variam sazonalmente e regionalmente. Ademais, as previsõesdas afluências futuras são, em geral, pouco precisas. A incerteza com respeito asvazões, aliada a incerteza com respeito à demanda de energia, faz do planejamento daoperação de sistemas hidrotérmicos um problema essencialmente estocástico.

• é acoplado no espaço, ou seja, há interdependência na operação de usinashidroelétricas, pois a quantidade de água liberada em uma usina afeta a operação deoutra situada a jusante.

• o valor da energia gerada por uma hidroelétrica somente pode ser medido em termosda economia resultante nos custos de geração térmica ou déficits evitados e nãodiretamente como uma função apenas do estado da usina.

• os objetivos de economia de operação e confiabilidade de atendimento são claramenteantagônicos: a máxima utilização da energia hidroelétrica disponível a cada etapa é apolítica mais econômica, pois minimiza os custos de combustível. Entretanto, estapolítica é a menos confiável, pois resulta em maiores riscos de déficits futuros. Por sua



vez, a máxima confiabilidade de fornecimento é obtida conservando o nível dosreservatórios o mais elevado possível. Entretanto, isto significa utilizar mais geraçãotermoelétrica e, portanto, aumentar os custos de operação. O equilíbrio entre oscustos de operação e a confiabilidade é obtido através do custo do déficit, querepresenta o impacto econômico associado à interrupção do fornecimento. Adeterminação do custo do déficit é um problema muito complexo, porém fundamentalpara a determinação da política de operação mais adequada para o sistema. Se o custodo déficit é muito baixo, resulta em uma utilização excessiva dos reservatórios e,portanto, em maiores riscos de racionamento no futuro. Se o custo de déficit é muitoalto, resulta em uma utilização excessiva dos recursos termoelétricos do sistema e,portanto, em custos de operação elevados.

2.3 Sistemas Interligados

A existência de interligações com os sistemas vizinhos permite uma redução dos custos deoperação, através do intercâmbio de energia e um aumento da confiabilidade de fornecimento,através da repartição das reservas.No caso de sistemas puramente térmicos, os custos de combustível proporcionam ummecanismo natural de coordenação para a troca de energia entre os sistemas. Se o custo deoperação da térmica mais cara operando no sistema A (a térmica marginal) é US$ 45/MWh eo custo correspondente no sistema B é US$ 40/MWh, é intuitivo que a operação maiseconômica para o sistema A é receber energia de B.

Observa-se que os intercâmbios de energia entre os sistemas térmicos resultam na otimizaçãoglobal dos custos de operação do sistema interligado. Em outras palavras, os resultados que seobtêm com sistemas independentes fazendo intercâmbios com base nos custos das térmicasmarginais são os mesmos que se obteriam se os sistemas fossem operados de forma conjunta,como um único sistema. Também é importante observar que a coordenação entre os sistemasse faz unicamente através de custos, isto é, o sistema B não tem informações sobre o sistemaA (tipo e disponibilidade dos equipamentos, demanda, etc.) ou vice-versa. Isto facilita oscontratos de intercâmbio, tendo em vista que não é necessário compartilhar informaçõescomerciais.

No caso de sistemas hidrotérmicos, é necessário inicialmente determinar o valor da geraçãohidroelétrica. Como discutido nas seções anteriores, o valor da energia hidroelétrica é o valorda geração térmica que se poderia substituir hoje ou no futuro. Este valor se calcula comouma etapa do processo de determinação da política ótima.

Com este conceito, pode-se representar uma hidroelétrica como uma “térmica” cujo “customarginal de operação” (CMO) é o valor da água. Entretanto, é importante observar que estevalor não se mede de maneira isolada em cada usina, pois depende da operação conjunta dosistema. Em outras palavras, se a política ótima hidrotérmica de cada empresa for calculadade maneira isolada, os intercâmbios de energia posteriores, mesmo baseados nos valores daágua de cada empresa, não resultam na operação mais econômica possível.

Em resumo, para obter os ganhos operativos máximos de um sistema hidrotérmicointerligado, é necessário operar o sistema de maneira integrada, isto é, otimizando de formaconjunta a operação de todos os subsistemas, com o objetivo de minimizar o custo deoperação total.



3 PROGRAMAÇÃO DIÁRIA DA OPERAÇÃO - MODELODESSEM

3.1 Representação Temporal

O modelo DESSEM determina o despacho ótimo para a programação diária de sistemashidrotérmicos interligados, que minimiza o custo total de operação, composto das parcelas decusto presente, ao longo do período de estudo (geração térmica, contratos com sistemasexternos e déficit de energia) e custo futuro (expectativa de geração térmica e déficit deenergia no futuro), cuja função é fornecida pelo modelo DECOMP e acoplada ao final dohorizonte de estudo.

O horizonte de estudo é de até 14 dias, com possibilidades de discretização em intervalos demeia-hora, 1 hora ou em patamares cronológicos de duração variável. Cada intervalo(instante) de tempo na discretização será tratado neste documento pelo termo estágio.

Pode-se dividir o período de estudo no DESSEM em duas partes:

• Período de interesse: pode incluir até 5 dias de estudo, para os quais se utilizaria umadiscretização de meia hora;

• Período estendido: envolveria os demais dias de estudo, nos quais o nível dediscretização seria de 1 hora.

Caso seja de interesse, o intervalo de discretização em cada um destes períodos pode serampliado, de forma a considerar intervalos de tempo maiores, conforma mostra a figura 3.1:

Figura 3.1 – Discretização temporal do modelo DESSEM

A representação do estudo em intervalos maiores de tempo é feita por meio de patamarescronológicos. Para cada um destes patamares, os dados do sistema seriam constituídos poruma média dos valores verificados ao longo das horas que constituem o patamar. Para o casoda carga, por exemplo, seria feita uma agregação da curva de carga horária em intervaloscronológicos superiores a uma hora, conforme mostra a figura 3.2:

D

Período de interesse: até 5 dias Período estendido

Intervalo de meia hora ou superior Intervalo horário ou superior

Horizonte Variável de Estudo: até 14 dias

Início do estudo

FCF do DECOMP



Figura 3.2 – Representação do primeiro dia de estudo do DESSEM (t = 1 até 24h)em cinco patamares cronológicos, t’=1 até 5.

Ao final do período de estudo do modelo DESSEM, acopla-se a função de custo futuro (FCF)produzida pelo modelo de planejamento de curto prazo (DECOMP).

3.2 Representação do Sistema

O modelo considera de forma individual as usinas hidroelétricas e termoelétricas e, emsituações em que se tornar adequado, poderá fazer a representação por unidade hidro outermoelétrica2. Devido à necessidade de se contemplar a regularização diária dosreservatórios, praticamente todas as usinas a fio d’água no DECOMP serão tratadas comoreservatórios no modelo DESSEM. Outra diferença grande do DESSEM em relação aoDECOMP é que, no modelo DESSEM, praticamente em todos os aproveitamentos seráconsiderado o tempo de viagem da água para a usina de jusante.

A transmissão pode ser modelada de maneira simplificada, com a demanda por subsistema -SE, S, NE, N - e os intercâmbios entre eles (estudos sem rede), ou de forma detalhada(estudos com rede), com a representação de cada circuito e a carga por barra, através de umarepresentação DC da rede elétrica.

As afluências as usinas são consideradas determinísticas ao longo de todo o período deestudo.

As diversas características/restrições do problema da programação diária podem ser divididosem 3 tipos:

• características/restrições do sistema, que envolvem conjuntamente as usinashidroelétricas e termoelétricas;

• características/restrições para as usinas hidroelétricas;

• características/restrições para as usinas térmicas.

A descrição e detalhamento das características/restrições de cada tipo são feitas nos itens 5 a7.

3.3 Formulação do Problema

O problema a ser resolvido pode ser resumido matematicamente da seguinte forma:

2 Por exemplo, em estudos com rede elétrica, se unidades de uma mesma usina estiverem conectadas em barrasdiferentes do sistema

t

MW

t’=1 t’= 2 t’= 3 t’= 4 t’= 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

COSE

Note

O DESSEM é praticamente o DECOMP com as adaptações decorrentes da redução da discretização.

COSE

Note

Para uma programação em base horária ou de meia-hora, a representação da rede elétrica tem que ser a mais detalhada possível.



Tttestágiocadapararestriçõesticascaracterísasujeito

VEcontrDefctGTCMinimize TT

T

ttttt

,,1,/:

)( ),,( 1

K=

+∑

=α

(3.1)

onde o termo ),,( tttt EcontrDefctGTC refere-se ao custo imediato (ou presente), que é funçãodas gerações térmicas (GT), déficits de energia (Defct) e energia importada /exportada(Econtr) com sistemas externos, para cada estágio t do horizonte de estudo (consideram-se Testágios no estudo). O termo )( TT Vα representa o custo futuro em função do vetor dearmazenamentos VT para as usinas do sistema ao final do estudo. Esta função é fornecida pelomodelo DECOMP e sua formulação será detalhada mais adiante. A soma dos custos presentesde todos os estágios mais o custo futuro é denominada de Função Objetivo do problema.

Pode-se entender o problema como sendo composto de T subproblemas (um para cadainstante de tempo do estudo), que no entanto encontram-se acoplados pelas restrições queenvolvem a interdependência entre as operações dos diversos estágios3, conforme mostra afigura 3.3.

Figura 3.3 – Representação do problema da programação diária como uma seqüência desubproblemas para cada estágio, com interdependência temporal entre eles, e

acoplamento ao final com a função de custo futuro do modelo DECOMP

A FCF do modelo DECOMP é resultado da resolução do problema de curto prazo. Pode-seutilizar a função da primeira ou segunda semana de estudo do DECOMP, dependendo dohorizonte em que se está fazendo o estudo no DESSEM.

Esta função, denotada por αT(VT), é multidimensional, pois VT é o vetor de armazenamentosem todos os reservatórios do sistema ao final do estudo do DESSEM (final do estagio T).Assim, sua representação no plano (Figura 3.4 a seguir) é apenas ilustrativa.

3 pode-se citar, como exemplos de tais restrições, o balanço hídrico das usinas (volume final de um estágio igualao volume inicial do estágio seguinte), o tempo de viagem da água entre as usinas, e as restrições de variação degeração/defluência das usinas entre estágios consecutivos.

...1

V

$

2 T-1 T

FCF do DECOMPSubproblemas de cada estágio

Problema como um todo



Figura 3.4 – Gráfico ilustrativo da função de custo futuro do modelo DECOMP

A FCF é portanto uma função linear por partes, ou seja, para um dado vetor dearmazenamentos V*, o valor da função corresponde ao menor valor que seja igual ou superiorao valor que todos os segmentos4 que compõem a função assumem ao cruzar o ponto V*.Matematicamente, esta função assume a expressão:

αα minimize)( =TT V

sujeito a

∑=

+≥NR

iT iVi

111 )()( δπα

∑=

+≥NR

iT iVi

122 )()( δπα

...

∑=

+≥NR

ipTp iVi

1)()( δπα

onde:

NR: número de reservatórios no estudo;p: número de segmentos que compõem a função de custo futuro;α: variável escalar auxiliar para denotar αT;

VT(i) Volume, ao final do estudo do DESSEM (estágio T), para o i-ésimo reservatório nosistema (é o i-ésimocomponente do vetor VT);

πk(i) “inclinação” do segmento k, em relação ao eixo em que está representado o volume doreservatório i. O módulo deste valor é conhecido como “valor da água para a usina”, erepresenta o decréscimo no custo futuro considerando-se isoladamente um acréscimounitário de água no reservatório i;

δk termo independente para a equação do k-ésimo segmento que compõe a função. 4 Para o desenho mostrado, cada segmento é uma reta. No entanto, na formulação geral da função para NRreservatórios, cada segmento é um hiperplano definido no espaço ℜ(NR+1).

VT

αT(VT)

V*

α (V*)



A FCF indica então como varia o valor esperado do custo futuro de operação (após ohorizonte de estudo do DESSEM) em relação ao estado em que se atingir os reservatórios dosistema no final do estudo do DESSEM. Assim, a menor ou maior utilização da água em cadareservatório pelo modelo DESSEM estará relacionada à avaliação, além das condiçõesoperativas de cada estágio, de como varia a FCF de acordo com o esvaziamento / enchimentodos reservatórios do sistema.

3.4 Estratégia de Solução

A estratégia de solução adotada envolve inicialmente a decomposição do problema em Tsubproblemas, um para cada estágio t, onde cada subproblema contém apenas as restrições doestágio a que se referem, modeladas de forma linear. Os subproblemas são então resolvidosem sequência por um pacote de otimização OSL [5], que se baseia no algoritmo Simplex deProgramação Linear [6]. A otimização do problema como um todo, ou seja, considerando ainterdependência entre as operações de todos os estágios, é feita através da técnica deProgramação Dinâmica Dual Determinística (PDDD) [7], descrita com detalhes no anexo Ideste manual.

A técnica de PDDD consiste em realizar sucessivas iterações, cada uma composta de umasimulação no sentido direto (simulação forward, resolvendo os subproblemas dos estágios 1ao T-1), seguida de uma simulação no sentido inverso (recursão backward), resolvendo ossubproblemas dos estágios T ao 2. Durante esta recursão backward, constrói-se umaaproximação linear para a FCF de cada estágio t de T-1 a 1, em função dos resultados daresolução do subproblema para o estágio seguinte na própria recursão. Como resultado desucessivas iterações, obtém-se uma FCF para cada estágio, linear por partes, que indica, paracada vetor de armazenamentos nos reservatórios do sistema ao final do estágio, o custo futurode operação, incluindo não só o período até o final do horizonte de estudo do DESSEM, mastambém o período após este horizonte.

Ao final de cada iteração, obtém-se um limite inferior (custo presente + custo futuro doestágio 1) e um limite superior (somatório dos custos presentes dos estágios 1 a T-1 + custototal do estágio T) para o valor da função objetivo do problema como um todo. Aconvergência é atingida quando a diferença entre estes dois valores torna-se inferior a umacerta tolerância de otimalidade, definida pelo usuário.

Depois de atingida a convergência da PDDD, o modelo realiza uma simulação final, na qualresolve em seqüência os subproblemas dos estágios 1 a T, obtendo os resultados da operaçãopara publicação.



O fluxograma da Figura 3.5 esquematiza o algoritmo de resolução do problema deprogramação diaria, resolvido pelo modelo DESSEM com a aplicação da PDDD.

Figura 3.5 - Fluxograma representativo da resolução do problema de programação diária pelo modelo DESSEM

O Anexo II deste relatório descreve a metodologia da PDDD aplicada a um problema deprogramação linear multi-estágio genérico.

1ª iteração (k = 1)

k = k + 1 Simulação forward:Resolve em sequência estágios t = 2 até Τ (naprimeira iteração, resolve a partir do estágio 1).

Obtem-se LSUP: custo imediato dos estágios 1a Τ + custo futuro do estágio Τ

Recursão backward:Resolve em sequência estágios t = T-1 até 1,

obtendo novos cortes para a FCF de cada estágio

Obtem-se LINF: custo imediato + custo futurodo estágio 1

INÍCIO

N

S

(LSUP-LINF) < TOL ? LINF

Teste deConvergência

Simulação FinalResolve em sequência estágios t = 1 até Τ.Obtem-resultados da operação para publicação

FIM



4 TIPOS DE FUNCIONALIDADES DO MODELO DESSEMAs funcionalidades existentes no modelo DESSEM podem ser divididas em três tipos, deacordo com a natureza dos elementos considerados:

Funcionalidades envolvendo características/restrições do sistema:

• Carga por subsistema (estudos sem rede elétrica)• Intercâmbio entre subsistemas• Custo de déficit por profundidade de corte de carga• Restrição de transporte entre Itaipu 50 e 60 Hz e os subsistemas SU/SE• Contratos de importação/exportação de energia• Consideração de gerações fixas (gerações não otimizadas ou simuladas)• Consumo nas unidades elevatórias• Restrições elétricas

Funcionalidades envolvendo características/restrições para as usinas hidroelétricas:

Referentes ao balanço hídrico das usinas:

• Afluências naturais determinísticas aos aproveitamentos• Evaporação• Retiradas de água para outro usos• Limites físicos de armazenamento nos reservatórios• Enchimento de volume morto• Tempo de viagem da água entre aproveitamentos• Representação das unidades elevatórias• Vertimento máximo nos reservatórios

Referentes às restrições operativas para as usinas:

• Limites de armazenamento dos reservatórios• Volume de espera para controle de cheias• Limites de turbinamento, vertimento, desvio e defluência total em cada usina hidroelétrica• Limites de bombeamento nas unidades elevatórias

Referentes às gerações das usinas:

• Função de produção energética para usinas hidroelétricas• Reserva de potência por usina• Manutenção das unidades geradoras• Variação da vazão turbinada entre estágios consecutivos• Variação da geração entre estágios consecutivos (“bang-bang”)

Referentes às particularidades do sistema:

• Representação de Unidades Elevatórias• Variação diária e horária na cota da Régua 11 (a jusante de Itaipu)• Operação do canal Pereira Barreto (restrições de desnível máximo entre Ilha Solteira e

Três Irmãos)• Operação das bacias do Alto Tietê e Paraíba do Sul na otimização ou por um módulo de

simulação

COSE

Note

Esse tipo de restrição é um indício da baixa qualidade da solução obtida pelo modelo, provavelmente culpa de sua simplificação linear



Funcionalidades envolvendo características/restrições para as usinas térmicas:

• Rampa para tomada/alívio de carga

• Reserva de potência por usina• Manutenção das unidades geradoras

• Custos e limites físicos de geração por unidade

• Limites operativos de geração para usinas e unidades

As funcionalidades desses três tipos são detalhadas respectivamente nas seções 5, 6 e 7 queseguem. Cada uma destas seções é dividida em duas partes: na primeira é feita a descrição dascaracterísticas / restrições e na segunda é apresentada a formulação matemática tal comoimplementada no modelo.

Descreve-se à parte, no item 8, o tratamento das características/restrições referentes à redeelétrica, que envolvem dois pontos básicos:• Carga por barra na configuração da rede elétrica

• Representação da rede elétrica por um fluxo DC, considerando restrições de violação doslimites individuais dos circuitos ou de somatório de fluxos nos circuitos

A descrição das notações utilizadas na formulação matemática é apresentada no Anexo I.



5 FUNCIONALIDADES ENVOLVENDO CARACTERÍSTICAS /RESTRIÇÕES DO SISTEMA

Foram incluídas neste tipo as restrições que promovem um acoplamento entre as usinashidroelétricas e térmicas do sistema, ou seja, quando ambos os tipos de usinas estãoenvolvidos na mesma restrição.

5.1 Descrição

5.1.1 Carga por subsistema (estudos sem rede elétrica)Considera-se a demanda por subsistema, para cada estágio do estudo. Assim, é possívelrepresentar a curva de carga horária ou, no caso de representação em patamares cronológicos,a agregação desta curva em intervalos de tempo maiores.

5.1.2 Intercâmbio entre subsistemasRepresentam-se, nos estudos sem rede, os intercâmbios entre os diversos subsistemas, para osquais pode-se definir limites mínimos e máximos de fluxo de energia, para cada estágio doestudo. Nos estudos com rede, embora esta representação seja desnecessária (pois os próprioscircuitos que compõem a interligação estarão sendo considerados na rede elétrica), em algunscasos é interessante se manter um limite para os intercâmbios entre os subsistemas5, paraevitar que, durante a busca do despacho ótimo, o modelo realize despachos preliminaresmuito acima dos limites das interligações, encontrando muitas violações de circuitos ao longodo processo, conforme mostram os resultados de [8].

5.1.3 Custo de déficit por profundidade de corte de cargaO modelo permite que o custo do não atendimento à carga (déficit de energia), em cadasubsistema, por estágio, seja fornecido através de uma curva linear por partes, que relacionacustos distintos às diferentes parcelas de carga não atendida.

Profundidade dedéficit (%)

5 15 25 100

Cdefct (Edefct ) ($)

1200 R$/MWh

1800 R$/MWh

2900 R$/MWh

600 R$/MWh

Figura 5.1 – Exemplo de representação de uma curva de custo de déficit por patamares

Observe que, à medida em que o percentual de corte de carga cresce, o custo unitário dodéfciti aumenta também, pois cargas mais importantes vão deixando de ser atendidas.

5 Ressalta-se, no entanto, que esses limites devem estar coerentes com os limites dos circuitos que perfazem asinterligações



5.1.4 Restrição de transporte entre Itaipu 50 e 60 Hz e os subsistemas SU/SEA representação da capacidade de transporte de energia entre Itaipu e os subsistemas Sul eSudeste é realizada pela incorporação opcional de um conjunto de restrições específico paraeste fim, que incluem a incorporação do nó de Ivaiporã ao sistema. A Figura 5.2 ilustra arepresentação da interligação de Itaipu com a inclusão desse nó. Ressalta-se que a usina estásituada na região Sudeste, mas a sua geração total é repartida, por meio de restrições, entre aslinhas de transmissão de 50Hz e 60Hz. A carga da ANDE está incorporada à geração 50Hz.

60Hz

Itaipu (IT)

ANDE

50Hz

Ivaiporã (IV)

SE

S

Figura 5.2 – Representação da interligação de Itaipu

5.1.5 Contratos de importação/exportação de energiaA incorporação de contratos de importação e/ou exportação de energia permite arepresentação de contratos de compra e/ou venda de energia entre os subsistemas em estudo eoutros subsistemas cuja operação não esteja representada no estudo (ex. Argentina).

Os contratos de importação / exportação são representados no problema de otimizaçãoconsiderando a convenção de sinal positivo / negativo, respectivamente, para importação /exportação, na função objetivo e restrição de atendimento à demanda do subsistemacontratante.

É possível definir, para cada estágio, custos e limites mínimo e máximo de energia contratadapor contrato.

5.1.6 Consideração de gerações fixas (gerações não otimizadas ou simuladas)Podem existir no sistema usinas de pequeno porte e bacias especiais, cujas inclusões dentrodo processo de otimização centralizada do sistema se tornem inadequadas ou desnecessárias.As gerações fixas, por estágio, correspondentes a essas usinas e bacias podem ser fornecidasao modelo DESSEM para serem incluídas nas equações de atendimento a demanda, comoabatimento de mercado.

5.1.7 Consumo nas unidades elevatóriasÉ o consumo específico de energia por unidade de vazão bombeada ((MW/(m3/s)) para umaunidade elevatória (vide item 6.1.7). O consumo total de energia em cada estágio é computadocomo acréscimo de demanda no subsistema em que se situa a unidade.

5.1.8 Restrições elétricasExistem pontos no sistema elétrico brasileiro que merecem especial atenção, por problemas detransmissão ou por sobrecarga. Para considerar estes pontos, nos estudos sem rede, foramincluídas restrições elétricas especiais, que traduzem limitações de geração em conjuntos deusinas hidroelétricas e/ou usinas térmicas. Cada restrição elétrica considera limites inferior e



superior de geração para um conjunto de usinas com fatores de participação diferentes. Essasrestrições são consideradas por estágio.

Mesmo nos estudos com rede pode ser necessário incluir algumas deste tipo, para representarrestrições relativas a controle de tensão ou suporte de reativo, que não são vistas pelamodelagem DC implementada para a rede elétrica.

5.2 Modelagem

Descreve-se a seguir as expressões matemáticas que modelam e relacionam asfuncionalidades descritas no item anterior.

5.2.1 Termos na função objetivo

• Custos de contrato

∑ ∑= =

T

t

NCIE

itt iEcontr(i)Ccontr

1 1)( .

• Custo de déficit

)( 1 1

∑ ∑= =

T

t

NSIST

ktt )k(EdefctCdefc

5.2.2 Equação de Atendimento a demandaA demanda deve ser atendida (ou complementada com déficit) para cada subsistema, a cadaestágio de tempo no estudo. A equação de atendimento a demanda relaciona vários dos termosdescritos anteriormente:

( )

TtNSIST,k

kDkEdefctjEcontrjGPQjGBE

rkfkrfjBombjCons -iGTiGH

tt

kNCIE

jt

kNPQ

jt

kNBE

jt

rtt

kNUNE

jt

kNUT

it

kNUH

it

k

1,..., , ...,1

, )( )( )( )( )(

),(),( )().( )( )(

)(

1

)(

1

)(

1

)(

1

)(

1

)(

1

==

=+±++

+−++

∑∑∑

∑∑∑∑

===

Ω∈===

5.2.3 Limites de intercâmbio

TtrNSISTkrkfrkf ktt ,...,1 , ,,...,1, ),(),(0 =Ω∈=≤≤

5.2.4 Restrição de transporte entre Itaipu 50 e 60 Hz e os subsistemas SU/SE

),(),(),(),()(60_ )(60_ )(50_ )(

IVSfIVSEfSIVfSEIVfItGHItGHItGHItGH

ttttt

ttt

−−+=+=

IV

SE S

It

GH_50t (It)(excluindo parceladestinada à ANDE)

ft (IV,SE)

ft (SE,IV)ft (IV,S)

ft (S,IV)

GH_60t (It)



5.2.5 Limites de contrato de importação/exportação de energia

TtNCIEjjEcontrjEcontrjEcontr ttt 1,..., ,,...,1 , )( )( )( ==≤≤

5.2.6 Restrições elétricas

TtNRE k

kLSREiGTkiiGHkikLIREkNUTR

itt

TtRE

kNUHR

it

HtREt

,...,1 ,,...,1

, )( )( ),( )( ),( )()(

1

)(

1

==

≤+≤ ∑∑==

κκ



6 FUNCIONALIDADES ENVOLVENDO CARACTERÍSTICAS /RESTRIÇÕES DAS USINAS HIDROELÉTRICAS

As funcionalidades das usinas hidroelétricas podem ser referentes a:

• balanço hídrico no sistema;

• restrições operativas;

• geração das usinas;

• particularidades do sistema.

Cada um desses grupos é descrito a seguir.

6.1 Descrição das Funcionalidades Referentes ao Balanço Hídrico noSistema

6.1.1 Afluências naturais determinísticas aos aproveitamentosO modelo DESSEM não considera cenários de afluências ao longo do horizonte de estudo. Asinformações determinísticas sobre as vazões naturais aos aproveitamentos são utilizadas paracalcular as respectivas vazões incrementais, e posterior inclusão nas equações de balançohídrico do sistema.

6.1.2 EvaporaçãoA representação da evaporação da água nos reservatórios é opcional. O usuário escolhe osreservatórios nos quais deverá ser considerada a evaporação, e o modelo calcula, a cadaestágio, o volume evaporado referente ao volume inicial do reservatório e à duração doestágio, descontando este valor no balanço hídrico da usina.

6.1.3 Retiradas de água para outros usosO modelo DESSEM permite representar desvios de água a montante dos aproveitamentos(para irrigação, abastecimento...) a cada estágio, conforme uma vazão definida pelo usuário.

6.1.4 Limites físicos de armazenamento nos reservatóriosAs informações cadastrais sobre os volumes mínimo e máximo dos reservatórios constituemos limites físicos de armazenamento. A diferença entre esses dois valores é denominadavolume útil. Para que as usinas tradicionalmente representadas no modelo DECOMP comofio d’água possam ser modeladas como reservatórios no DESSEM, é preciso atualizar estasinformações cadastrais (normalmente para estas usinas o volume mínimo é igual ao máximo).Além disso, deve-se fornecer os coeficientes do polinômio cota-volume para essas usinas.

6.1.5 Enchimento de volume mortoO enchimento de volume morto é representado através de uma vazão mínima paraarmazenamento em cada estágio. A partir desta vazão, o modelo calcula volumes metamínimos de modo que o volume armazenado no reservatório da usina em cada estágio seja, nomínimo, igual a estes valores. Com o objetivo de se manter uma vazão mínima defluentedurante o período de enchimento, também é possível informar os valores mínimos dedefluência para cada estágio. O prazo de enchimento de volume morto é um dado de entradado modelo.



6.1.6 Tempo de viagem da água entre aproveitamentosPara diversos aproveitamentos hidroelétricos não se pode considerar que todo o volumedefluente das usinas de montante, durante um estágio, esteja disponível a jusante no mesmoestágio. Neste caso, o volume afluente à usina de jusante pode ser obtido através dapropagação das defluências de montante. Uma maneira simplificada de realizar estapropagação consiste em considerar o tempo de viagem, isto é, intervalo de tempo médio queas defluências consomem até alcançar a usina a jusante.

O modelo DESSEM é capaz de representar o tempo de viagem da água entre aproveitamentosadjacentes em uma bacia hidroelétrica. Estes tempos de viagem são considerados atualmentenas restrições de balanço hídrico das usinas hidroelétricas, durante a etapa de simulação final.Encontra-se em etapa de conclusão a consideração, na FCF de cada estágio, dos termosreferentes às defluências das usinas com tempo de viagem para jusante (vide anexo II).

6.1.7 Representação de unidades elevatóriasUma unidade elevatória consiste em um equipamento utilizado para retirar água de umreservatório a uma cota inferior (reservatório de jusante) e envia-la para um reservatório emuma cota superior (reservatório de montante), isto é, desviar o caminho natural do cursod’água (vide Figura 6.1). Tal operação é vantajosa quando se obtém um ganho deprodutividade, ao se desviar a água em relação ao caminho natural desta.

Figura 6.1 - Exemplo de uma unidade elevatória, quebombeia água da usina A para a usina B.

O modelo utilizado para representar uma unidade elevatória consiste somente de uma variávelpara representar o volume bombeado, Bombt(j), j = 1, …, NUNE, onde NUNE é o númerototal de unidades elevatórias existentes.

Cada unidade elevatória é considerada nas restrições de balanço hídrico das usinas demontante e jusante à ela. Nestas restrições, o volume bombeado em uma determinada unidadeelevatória é acrescido ao volume afluente no estágio para a usina de montante, e deduzido dovolume final para a usina de jusante:

VAt+1 = VA

t + VaflAt – QA

t – SAt – Bombt(j)

VBt+1 = VB

t + VaflBt – QB

t – SBt + Bombt(j)

onde Vafl indica o volume afluente às usinas, sejam ela incrementais ou resultantes daoperação de usinas a montante.

O consumo de energia das unidades elevatórias é considerado na equaçã de demanda dosubsistema a que pertencem (item 5.2.2).

A

B

UnidadeElevatória j

QA + SA

QB + SB

Bomb(j)

VaflA

VaflB



6.1.8 Vertimento máximo nos reservatóriosUma usina somente pode verter se o volume armazenado em um determinado estágio estiveracima do volume referente à oleira do vertedouro, com execção das usinas que apresentamdescarga de fundo. As restrições de vertimento máximo são impostas para limitar aquantidade de água a ser vertida por cada usina hidroelétrica em cada estágio. Leva-se emconsideração as afluências e defluências ocorridas no estágio e a quantidade de água noreservatório disponível a partir da cota da crista do vertedor. A formulação detalhada destarestrição encontra-se em [9].

6.2 Descrição das Funcionalidades Referentes às RestriçõesOperativas

6.2.1 Limites de armazenamentoO modelo DESSEM possibilita, em cada estágio, a operação dos reservatórios através delimites mínimo e máximo para controle do volume armazenado nos mesmos.

6.2.2 Volume de espera para controle de cheiasDe modo a considerar as restrições operativas para controle de cheias, o modelo DESSEMrepresenta os volumes de espera para amortecimento de cheias definidos a priori. Estesvolumes são informados ao modelo através do fornecimento do máximo volume armazenávelpara cada usina em cada estágio.

6.2.3 Limites de turbinamento, vertimento, desvio e defluência totalO modelo DESSEM possibilita, em cada estágio, o controle das vazões turbinada, vertida,desviada e defluente total das usinas hidroelétricas, através de limites mínimo e máximo paraessas vazões.

6.2.4 Limites de bombeamento nas unidades elevatóriasO modelo possibilita a operação das unidades elevatórias através de limites mínimo e máximopara controle do volume bombeado nas mesmas, em cada estágio.

6.3 Descrição das Funcionalidades Referentes à Geração das Usinas

6.3.1 Função de Produção EnergéticaA Função de Produção Energética (FPE) das usinas hidroelétricas faz a correspondência entrea quantidade de água que passa pelas turbinas e o montante de energia gerado.

A potência gerada por uma usina depende da vazão turbinada e da altura de queda, e esta porsua vez é função da diferença entre as cotas de montante e jusante (canal de fuga) no instantede tempo considerado. A cota de montante é obtida a partir do volume no reservatório, pelopolinômio volume-cota, enquanto a cota de jusante depende da vazão turbinada e, para asusinas em que o vertimento influencia na cota do canal de fuga, também da vazão vertida.Desta forma, a FPE é uma restrição não linear (e em muitos casos não convexa) envolvendo oo volume armazenado, a vazão turbinada e em alguns casos a vazão vertida.

Devido à não linearidade das FPE’s das usinas, sua representação no modelo DESSEM é feitaatravés de uma Função de Produção Energética Aproximada (FPEA), que é uma função linearpor partes, construída por meio de aproximações de Taylor de primeira ordem sobre a FPE, apartir de uma grade de pontos definida à priori para as variáveis armazenamento eturbinamento. Modela-se a FPEA considerando como região viável para o problema

COSE

Note

Uma das maiores limitações dos modelos da cadeia do CEPEL é aproximar a função de produção das usinas hidrelétricas por funções lineares por partes. Essa limitação é ainda pior na otimização de curtíssimo prazo, onde os rendimentos dos conjuntos turbina-gerador, altamente não lineares, são o principal objetivo da otimização.

COSE

Note

Ao definir a grade de pontos se está definindo os pontos candidatos a solução.



linear a região situada abaixo das diversas inequações lineares (cortes) que aproximam aFPE. Como a função objetivo do problema é de minimização de custos, espera-se que aoperação decidida pelo modelo respeite estritamente os pontos situados exatamente sobre oscortes.

A figura 6.2 exemplifica a construção de uma FPEA a partir da FPE original, considerandoapenas a variável turbinamento:

Figura 6.2 - Exemplo da construção de uma FPEA a partir da curva FPE

Se a FPE original não for côncava, a região abaixo desta curva torna-se não convexa, e destaforma as aproximações lineares por partes não conseguem representar adequadamente ospontos abaixo da curva, havendo perda de informação em relação à FPE. Assim, torna-senecessário determinar a envoltória convexa da região definida pela FPEA, fazendo entãomodificações na FPE original. Só depois deste procedimento é que se realizam asaproximações lineares. É importante ainda minimizar uma função dos desvios entre a FPE e aFPEA até então obtida, para corrigir eventuais distorções entre as duas curvas que tenhamsido introduzidas ao longo do processo.

Ressalta-se que, como a FPE pode ser influenciada por mais de uma variável, o procedimentoé realizado em um espaço multidimensional, o que torna os cálculos mais complexos. Umdetalhamento maior da modelagem da função de produção será apresentado e documentadoem breve pelo CEPEL.

6.3.2 Reserva de potênciaO modelo considera, em cada estágio, a reserva de potência alocada por usina hidroelétrica,limitando sua geração máxima. O usuário escolhe as usinas para as quais se alocará reserva.Caso seja de interesse, a extensão da modelagem para a alocação de reserva por área éimediata

6.3.3 Manutenção programadaConsideram-se atualmente os cronogramas de manutenção programada de cada unidadegeradora das usinas. Caso se deseje que a representação seja feita por usina, seria necessáriauma modificação na modelagem existente, adotando-se um fator de capacidade por estágiopara a usina como um todo, substituindo a representação das manutenções isoladas de cadaunidade geradora ao longo do período de estudo.

Q(m3/s)

FPE

FPEAGH(MW)

Região viável

COSE

Note

O modelo não aloca a reserva de potência no sistema. Ela já vem definida por usina.

COSE

Note

A função de produção, que é motivo da linearização por partes, não considera a variação do rendimento da turbina, do gerador e das perdas hidráulicas com o aumento da vazão turbinada.



6.3.4 Variação da vazão turbinada entre estágios consecutivosO modelo possibilita o controle da variação da vazão turbinada por uma usina em doisestágios consecutivos. Este recurso tem o objetivo de garantir que a vazão turbinada não variebruscamente de um intervalo de tempo para outro.

6.3.5 Variação da geração entre estágios consecutivosO modelo possibilita o controle da variação da geração de uma usina em dois estágiosconsecutivos, para evitar oscilações bruscas indesejáveis para a operação (“bang-bang”). Umestudo minucioso desta questão é apresentado em [11].

6.4 Descrição das Funcionalidades Referentes às Particularidades doSistema

6.4.1 Variação diária e horária na cota da Régua 11A área de confluência dos rios Paraná, Iguaçu, Acaray e Monday, junto à seção R-11 (Régua11), por representar a área fornteiriça do Brasil, Paraguai e Argentina, está sujeita à restriçõesimpostas por acordos internacionais, como o Acordo Tripartite. Existem também restriçõesdevido a ocupação desordenada nas margens do rio Paraná e da existência de uma usinahidroelétrica imediatamente a jusante de Itaipu. A figura 6.2 apresenta o esquema desta regiãode confluência.

Figura 6.2 - Localização da Régua 11.

No modelo DESSEM, as restrições associadas à Régua 11 são consideradas através dainclusão das restrições de variação horária e/ou diária no nível da Régua 11. Por exemplo, nocaso de uma restrição de variação horária, o nível na Régua 11 pode variar no máximo meiometro ao longo do período de uma hora. Já no caso de uma restrição de variação diária, onível na Régua 11 pode variar no máximo dois metros ao longo de qualquer intervalo igual ouinferior a 24 horas compreendido no período de estudo.

6.4.2 Operação do canal Pereira BarretoA operação do canal Pereira Barreto, que localiza-se entre as usinas hidroelétricas de IlhaSolteira e Três Irmãos, é realizada através da consulta a uma tabela de dados fornecida pelaCESP para controlar o valor da vazão e o seu sentido no canal. Esta tabela fornece valores devazão para o canal em função do desnível entre as cotas dos dois reservatórios e do valor damaior cota entre eles. A tabela apresenta uma região cujos níveis conduzem a velocidadesproibitivas de vazão para a navegação.

RÉGUA 11

Salto Caxias

Itaipu

Porto Primavera

Rio Iguaçu

Rio Paraná

COSE

Note

Tipo de restrição que indica a baixa qualidade da solução devido a linearidade do modelo.

COSE

Note

Idem.



A cada estágio, o modelo DESSEM determina a vazão no canal e o seu o sentido (entre IlhaSolteira e Três Irmãos), realizando um controle de modo a não atingir as regiões proibidas.

6.4.3 Operação das bacias do Alto Tietê e Paraíba do SulAlgumas bacias hidrográficas do sistema interligado brasileiro, tais como as bacias do Paraíbado Sul e Alto Tietê, têm sua operação hidroenergética determinada através de procedimentosde simulação específicos que não possuem o objetivo de minimizar o custo de operação dosistema. O Modelo DESSEM permite que estas duas bacias sejam tratadas dentro do processode otimização (juntamente com todas as outras usinas do sistema) ou por algoritmos desimulação.

No caso da simulação, a geração destas bacias é obtida através de simulação da operaçãoconsiderando as regras específicas de operação das mesmas. Para a bacia do Paraíba do Sul,foram incorporadas as regras do modelo SHP “Simulador Hidroenergético do Paraíba doSul”, desenvolvido na divisão de Operação Energética da Eletrobrás. No início de cadaestágio é calculada a operação do Paraíba do Sul através destas regras, retornando-se aoproblema as perdas energéticas nas unidades elevatórias e as gerações nas usinashidroelétricas. Tanto as perdas energéticas quanto a geração são então incluídas na equação deatendimento a demanda do estágio (item 5.2.2). Mais detalhes das regras de operação da baciado Paraíba do Sul podem ser obtidos em [10].

Para a bacia do Alto Tietê, foi também incorporada a simulação das regras específicas para aoperação desta bacia. No início de cada estágio é calculada a operação do Alto Tietê atravésdestas regras, retornando-se ao problema as perdas energéticas nas unidades elevatórias e asgerações nas usinas hidroelétricas. Tanto as perdas quanto a geração são então incluídas naequação de atendimento a demanda do estágio.

Uma forma alternativa de representar aproximadamente as regras de operação dessas baciasno processo de otimização da operação do sistema interligado é através da inclusão derestrições lineares e um novo elemento, as unidades elevatórias. As restrições especificamlimites mínimo e máximo para vazão defluente total, vertida, desviada e volume armazenadono reservatório. Além disso, as informações das unidades elevatórias, tais como,bombeamentos máximo e mínimo, e taxa de consumo são incluídas no problema.

6.5 Modelagem


6.5.1 Termos na função objetivoNão há termos na função objetivo envolvendo a operação das usinas hidroelétricas, apenas afunção de custo futuro, já detalhada na seção 3.3.

6.5.2 Equação de balanço hídricoÉ uma equação para cada usina no estudo (seja ela reservatório ou fio d’água), em cadainstante de tempo, englobando vários dos termos descritos anteriormente.



( ) ( )

( ) ( )

,...,1 ,...,1

, )( )( )( )()()(

)()( )()()( )( )( )(

)()()(

)(

)(1

T tNUH,i

iViViVjSjQjDv

jSjQiDviSiQiAiViV

tttiMtvj

jtvtjtvtiMj

t

iMjtttttttt

enchimevapalterntDV

t

==

−−−++

+++++−+=

∑∑

∑

∈−−

∈

∈+

6.5.3 Restrição de vertimento máximoÉ uma equação para cada reservatório no estudo, em cada instante de tempo, englobandovários dos termos descritos anteriormente.

( ) ( )

( )

,...1,,1),( )(

)()( )()()( )( )( )(

)(

)()()()(

TtNRiiVcristjDv

jSjQjSjQiQiAiViS

iMjt

iMjtt

iMtvjjtvtjtvtttt

DV

tt

==−+

+++++−+≤

∑

∑∑

∈

∈∈−−

6.5.4 Volume evaporado

,...,1 ,...,1 , ))(( )( 1 T tNUH, i iVfiV tevaptevap === −

6.5.5 Limites de armazenamento físico e operacional

TtNUHiiViVt ,...,1 ,,...,1 , )( )( 0 ==≤≤

TtNUHiiViViVttt ,...,1 ,,...,1 , )( )( )( ==≤≤ supinf

6.5.6 Restrição de volume de esperaTtNUHiiiV

tt espera ,...,1 ,,...,1 , )(V )( 0 ==≤≤

6.5.7 Limites de turbinamento, vertimento, desvio de água por canal e defluênciatotal

TtNUHiiQiQiQ ttt ,...,1 ,,...,1 , )( )( )( ==≤≤ maxmin

TtNUHiiSiSiS ttt ,...,1 ,,...,1 , )( )( )( ==≤≤ maxmin

TtNUHiiDviDviDv Dvttt ,...,1 ,,...,1 , )( )( )( ==≤≤ maxmin

TtNUHiiDefiDviSiQiDef ttttt ,...,1 ,,...,1 , )( )()()( )( ==≤++≤ maxmin

6.5.8 Limites de bombeamentoTtNUNEjjBombjBombjBomb ttt ,...,1 ,,...,1 , )( )( )( ==≤≤

6.5.9 Limites de geração

TtNUHiiGHiGH tt ,...,1 ,,...,1 , )( )(0 ==≤≤

6.5.10 Restrição de reserva de potência

TtNUHiiRHiGHiGH ttt 1,..., ,,...,1 , )()( )( ==≥−

COSE

Note

Imposição de reserva de potência por usina!



6.5.11 Função de produção energética

( ) Tt,NUH,iiSiQiViVfiGH ttttiprodt ,...,1 ...,1 , )(),(),(),( )( 1

=== +

6.5.12 Restrições de variação da vazão turbinadaTtNUHiiQiQiQiQ

tttt ,...,1 ,,...,1 , )( )()()( sup1inf ==∆≤−≤∆− −

6.5.13 Restrições de variação horária e diária no nível da Régua 11

( ) ( ) T tItSItQkiusSiusQkQ ttIt

RtRtIg

tR RR,...,1 , )()( )()(

1111 111111 =+++=

TtQQQQ hRtRtR

hR ,...,1 , 111111111 =∆≤−≤∆− −

TtkQQQQ dRktRtR

dR ,...,1,24,...,1 , - 11111111 ==∆≤≤∆− −

)( 111111h

RRh

R HfQ ∆=∆

)( 111111d

RRd

R HfQ ∆=∆

6.5.14 Limite da vazão no canal Pereira BarretoT t VVfQ TRIILS ttTRIILStTRIILS ,...,1 ,),( )()( == −−

T t VVfQ TRIILS ttTRIILStTRIILS ,...,1 ,),( )()( =≤ −−



7 FUNCIONALIDADES DAS USINAS TÉRMICAS

7.1 Descrição

Algumas das restrições descritas e formuladas nesta seção representam as unidades térmicasindividualmente. Entretanto, a representação por usina térmica pode ser feita de formaanáloga, bastando agregar as restrições referentes às unidades de uma mesma usina.

7.1.1 Rampa para tomada/alívio de cargaExiste atualmente no modelo DESSEM a possibilidade de representação detalhada daoperação das unidades térmicas, considerando questões como custo de partidas das unidades,tempos de resfriamento e aquecimento das unidades e rampas máximas para alívio e tomadade carga por unidade. Esta modelagem não foi descrita neste relatório, pois este nível dedetalhamento para a operação das usinas não está sendo requerido pelo Submódulo 8.1,referência de estudo para a GTRD.

No entanto, é possível a implementação de uma modelagem mais simples, considerandoapenas rampas para variação máxima positiva e negativa na geração térmica entre doispatamares de tempo consecutivos.

7.1.2 Reserva de potênciaO modelo considera, em cada estágio, a reserva de potência alocada por usina térmica, emcada estágio de tempo. O usuário escolhe as usinas para as quais se alocará reserva.Caso sejade interesse, a extensão da modelagem para a alocação de reserva por área é imediata.

7.1.3 Manutenção programadaConsideram-se atualmente os cronogramas de manutenção programada de cada unidadegeradora das usinas. Caso se deseje que a representação seja feita por usina, seria necessáriauma modificação na modelagem existente, adotando-se um fator de capacidade por estágiopara a usina como um todo, substituindo a representação das manutenções isoladas de cadaunidade geradora ao longo do período de estudo.

7.1.4 Custos e limites físicos de geração por unidadeOs limites físicos de geração são informações cadastrais sobre a capacidade máxima degeração de cada unidade térmica de uma usina. Já os custos lineares de geração das unidadestérmicas (R$/MWh) são informações operativas que podem alterar de valor ao longo doperíodo de estudo, isto é, nos estágios.

7.1.5 Limites operativos de geração para usinas e unidadesO modelo DESSEM possibilita, em cada estágio, o controle dos limites mínimo e máximo degeração para as usinas térmicas e/ou suas unidades geradoras.

7.2 Modelagem


7.2.1 Termos na função objetivoCustos lineares das unidades térmicas:

COSE

Note

É evidente que num sistema predominantemente hidrelétrico como o brasileiro, onde o acompanhamento da carga pode ser feito exclusivamento com geração hidrelétrica, não há necessidade de modelagem elaborada da operação das termelétricas.



( )∑ ∑= =

T

t

NUT

it

it iGTCterm

1 1 )(

7.2.2 Limites de geração das usinas e unidades térmicas

TtNUTiiGTiGTiGT ttt 1,..., ,,...,1 )()()( ==≤≤

TtNUNTiigtigtigt ttt 1,..., ,,...,1 )()()( ==≤≤

7.2.3 Restrição de reserva de potência

TtNUTiiRTiGTiGT ttt 1,..., ,,...,1 , )()( )( ==≥−

ou

( ) TtNUTiiRTjgtjgtiNUNT

jttt 1,..., ,,...,1 , )( )()(

)(

1

==≥−∑=



8 CONSIDERAÇÃO DA REDE ELÉTRICA NO CÁLCULO DODESPACHO PELO MODELO DESSEM

Este item tem por objetivo explicar a metodologia implementada atualmente no modeloDESSEM para tratamento da rede elétrica segundo uma modelagem DC, considerando-se, nadeterminação do despacho, os limites de fluxos individuais nos circuitos e eventuais limitesde somatório de fluxos em vários circuitos. Maiores detalhes e um estudo de caso mostrandoalguns resultados obtidos com esta técnica podem ser obtidos em [9].

8.1 Formulação DC Da Rede Elétrica

O modelo linearizado em potência ativa, ou fluxo DC, fornece uma aproximação dadistribuição dos fluxos de potência ativa no sistema, no qual se despreza o efeito datensão/potência reativa. O fluxo DC é baseado no forte acoplamento entre a potência ativa e oângulo das tensões e apresenta resultados tanto melhores quanto maior for o nível de tensãodo sistema representado. Este modelo aproximado permite o cálculo dos fluxos ativos comrazoável precisão e com baixo custo computacional.

O modelo DC é obtido através da linearização das equações de fluxo de potência ativa narede. Desprezando-se as perdas, a equação de fluxo de potência ativa F em um circuito l entreduas barras k e m é dada por :

lF = -Vk Vm bl senθkm (8.1)

onde Vk e Vm representam as tensões das barras k e m, bl é a susceptância do circuito l e θkm adiferença angular entre k e m.

Considerando as aproximações Vk ≅ Vm ≅ 1 p.u., senθkm ≅θkm e bl ≅-1/xl, obtém-se:

l

mkkmll x

bFθθ

θ −=−= (8.2)

Como a injeção de potência ativa em cada barra é igual à soma dos fluxos que saem da barra,chega-se a uma formulação matricial do tipo:

p = Bθ (8.3)

onde p é o vetor de injeções de potência ativa (geração menos carga) por barra, θ é o vetor deângulos das tensões nodais e B a matriz de susceptâncias da rede.

Conhecendo p e B pode-se calcular, definida uma barra de referência para o sistema, osângulos de tensão θ em cada barra. Posteriormente, avaliam-se os fluxos nos circuitos pelaequação (8.2).

8.2 Descrição da Metodologia Adotada no DESSEM

Para se obter de forma direta, no subproblema de programação linear (PPL) de cada estágiodo DESSEM, um despacho ótimo que atenda aos limites de fluxo ativo nos circuitos da redeelétrica, seria preciso adicionar duas restrições por circuito, uma para cada sentido de fluxo.Em um sistema de grande porte, isto pode representar milhares de restrições.

Para contornar este inconveniente, adotou-se um esquema iterativo para a resolução dosubproblema para cada estágio, a cada iteração, a fim de se obter um fluxo de potênciaótimo DC. Este processo é composto dos seguintes passos:

COSE

Note

A parte elétrica é representada através modelo linear CC. Através de uma técnica de relaxação, as restrições de transmissão violadas vão sendo consideradas até que uma solução viável seja produzida.



Passo 1: Inicialmente, calcula-se um despacho considerando apenas limites de intercâmbioentre subsistemas. Embora estas restrições possam ser desconsideradas, sua inclusão éimportante para encontrar-se um despacho inicial que não esteja muito distante da soma doslimites nos circuitos que compõem os intercâmbios.

Passo 2: Conhecendo as barras onde se conectam as cargas e unidades geradoras hidráulicas etérmicas (ou qualquer outra injeção de potência, como um contrato de importação ouexportação de energia), têm-se o vetor P de cargas e gerações nas barras da rede.

Passo 3: Uma vez que B é construída a partir da entrada de dados do modelo, calculam-setodos os fluxos na rede resolvendo (8.2) e (8.3). Devido à restrição de atendimento à demandapor subsistema, incluída no PPL, o vetor P já está equilibrado (∑ cargas = ∑ gerações) e aresolução de (8.3) não leva a uma alteração na geração da barra de referência do sistema,mantendo-se então todas as gerações obtidas no PPL.

Passo 4: Comparando cada fluxo com um valor máximo nominal, dada uma tolerânciadesejada, identificam-se quais circuitos foram violados e, para cada um destes, adiciona-se aseguinte restrição no PPL:

wl • p ≤ lFmax , se a violação é no sentido k–m ,

ou (8.4)

wl • p ≥ - lFmax , se a violação é no sentido m–k

onde • denota um produto escalar6 e:

wl = bl (A B-1)l é o vetor de coeficientes (derivadas) da restrição de fluxo no circuito l emrelação às gerações nodais, possuindo dimensão NB (número de barras do sistema). O índice lno parêntesis indica a l-ésima linha na matriz AB-1, que corresponde ao circuito de mesmoíndice.

A = matriz incidência nó-ramo. Cada linha da matriz corresponde a um circuito e possuiapenas dois elementos não nulos: -1 na coluna correspondente à barra “de” do circuito e +1 nacoluna correspondente à barra “para” do circuito. Possui dimensão NC x NB, onde NC é onúmero de circuitos na rede

B-1 = inversa da matriz de susceptâncias (dimensão NB x NB)

Já que se supõem conhecidas as cargas, assim como algumas gerações “não otimizadas”(ambas são abatidas na demanda no PPL e entram no cálculo de p), vários dos termos àesquerda em (8.4) podem ser deslocados para o termo independente. Introduz-se esta equaçãono PPL aplicando a cada unidade geradora hidráulica ou térmica (ou contrato de importação)otimizado, o coeficiente da barra em que está conectado.

Passo 5. Resolve-se o novo PPL com a inclusão de novas restrições. Três situações poderãoocorrer:

a) O PPL se torna inviável, indicando que é impossível se realizar um despacho atendendo atodas as cargas e sem violar os circuitos dentro da tolerância especificada. => Nesta situação,é necessária a definição de prioridades de corte de carga no sistema, a fim de tornar odespacho viável.

6 Um produto escalar entre dois vetores a e b, definido quando ambos os vetores têm a mesma dimensão n, tem

sua expressão dada por: a • b= ∑=

n

iii ba

1. onde ai e bi são as coordenadas de ambos os vetores.



b) Um novo despacho é obtido, mas que pode levar a violações em outros circuitos, devido àredistribuição que teve que ser feita nas gerações anteriores => Neste caso volta-se ao passo 2,verificam-se os novos fluxos, inserem-se no PPL restrições para os novos circuitos violados etorna-se a resolver o PPL.

c) Encontra-se um despacho viável e que não viola nenhum outro circuito => Passa-se então àresolução do próximo estágio na PDDD.

A Figura 8.1 mostra o fluxograma do processo iterativo que se insere na resolução de cadaestágio, quando consideramos os limites da rede elétrica no estudo:

Figura 8.1 – Algoritmo de resolução do subproblema paracada estágio, considerando limites na rede

8.3 Observações

• Como a nova aproximação para a FCF do estágio t-1 só é feita quando se obtém umdespacho viável para o estágio, não estamos considerando soluções inviáveis no cálculodesta função.

• As restrições de violação nos fluxos adicionadas em determinado estágio e iteraçãopermanecem no PPL deste estágio para as iterações seguintes, sendo vistas previamentecomo restrições “potencialmente violáveis”.

Garante-se que o despacho final obtido, se o processo iterativo terminar com uma soluçãoviável no PPL de cada estágio, é uma solução de mínimo custo para o problema de despacho

despacho: gerações térmicas e hidroelétricas,energia importada / exportada, consumos nasusinas elevatórias

N

S

Adiciona ao PPLdo estágio

restrições paranovas violações

Resolve o PPL para estágio t, iteração k(PPL original + restrições de fluxo já adicionadas

nesta iteração ou em iterações anteriores)

Despachorespeita fluxosnos circuitos ?

- constrói cortes para (t-1), se estiver nabackward ou no estágio Τ da PDDD)

- passa para próximo estágio da PDDD

Despachoviável?

Definir critérios de corte de carga, afim de tornar o despacho viável

(não implementado)

S

N



energético com atendimento aos limites de fluxo DC nos circuitos. Note que quando seinserem novas equações de limite de fluxo no PPL de certo estágio ao longo do processoiterativo da PDDD, não estamos modificando o problema, pois na verdade estas restrições jáexistiam implicitamente, apenas não tinham ainda se tornado ativas.

O procedimento adotado é geral para os casos onde se tem várias gerações reunidas em umamesma barra, bastando contabilizar, na geração da barra para resolução do fluxo, o despachode todos os geradores, além de aplicar o fator referente a esta barra, nas inequações de limitesde fluxos nos circuitos construídas, igualmente a todos esses geradores.

Esta metodologia é adotada também para considerar limites de somatório de fluxos em umconjunto de circuitos. Para tal, compõe-se uma inequação semelhante à (8.4), onde o termo àesquerda compreende a soma das expressões de fluxo de todos os circuitos envolvidos narestrição e o termo à direita consiste no limite desejado, podendo-se eventualmente limitartambém o valor do fluxo no sentido oposto ao indicado.

Ressalta-se que a metodologia apresentada não considera as perdas elétricas nos circuitos. Ocálculo do despacho no modelo DESSEM com representação das perdas está em fase deimplementação, conforme mencionado no item 10 deste relatório.



9 REFERÊNCIAS

[1] Maria Elvira P. Maceira, Leslie A. Terry, Jorge M. Damazio, Fernanda S. Costa, Albert C.G. de Melo, Chain of Optimization Models for Setting the Energy Dispatch and Spot Pricein the Brazilian System, Power System Computation Conference - PSCC’02, Sevilla,Spain, June 24-28, 2002.

[2] CEPEL, Especificação Funcional do modelo NEWAVE, Relatório Técnico, 1999.

[3] CEPEL, Manual de Referência – Modelo DECOMP, versão 10.0 - Relatório Técnico,2002

[4] M.E.P.Maceira, L.A.Terry, A.S.L.Diniz, L.C.F. de Souza, F.S.Costa, S.P.Romero,S.Binato, S.M.Amado, C.E.Vilasboas, R.Vilanova, Despacho de Geração Horário comRepresentação Detalhada de Restrições Hidráulicas, VII SEPOPE, 2000.

[5] IBM, OSL-Optimization Subroutine Library - Guide and Reference Release 2.1 – IBM,1995

[6] Dantzig, Linear Programming and Extensions, Princeton University Press, New Jersey,1963

[7] CEPEL, Programação Dinâmica Dual Estocástica Aplicada ao Planejamento daOperação Energética de Sistemas Hidrotérmicos com Representação do ProcessoEstocástico de Afluências por Modelos Auto-regressivos Periódicos, Relatório Técnico,1993

[8] Diniz, A.S.L., Sousa, L.C.F., Maceira, M.E. et al, , Estratégia de Representação DC daRede Elétrica no Modelo de Despacho da Operação Energética - DESSEM , VIIISEPOPE, 2002

[9] CEPEL, Restrições de Vertimento Máximo em Usinas Hidrelétricas – Modelo DESSEMRelatório Técnico, 2000

[10] CEPEL, Descrição das Regras de Operação do Simulador Hidroenergético do Paraíbado Sul (SHP), Relatório Técnico, 2000.

[11] CEPEL, Análise das Variações nas Gerações das Usinas Hidroelétricas, Relatório Técnico para o ONS, 2002.

[12] Benders, J. F., “Partitionimg procedures for solving mixed variables programming problems”, Numer. Math., 4, 238-252, 1962.

[13] Lasdon, L. S., “Optimization theory for large systems”, McMillan Company, New York, 1970.

[14] Pereira, M.V.F., Pinto, L.M.V.G., “Stochastic optimization of a multireservoirhydroelectric system: a decomposition approach”, Water Res. Res., vol 21, n 6, 779-792,June, 1985.



Anexo I: NOTAÇÕES UTILIZADAS

Valor esperado do custo total de operação após o horizonte deplanejamento.

Afluência incremental à i-ésima usina hidroelétrica, ao longo do estágio t.

Bombeamento da j-ésima usina elevatória, no estágio t.

)( jBombt Bombeamento mínimo para a j-ésima usina elevatória, no estágio t.

)( jBombt Bombeamento máximo para a j-ésima usina elevatória, no estágio t.

Custo da energia contratada no j-ésimo contrato de importação/exportaçãode energia, no estágio t (custo negativo para energia exportada).

Custo de déficit de energia no k-ésimo subsistema, no estágio t.

Termo independente da k-ésima equação de função de produção da usinahidroelétrica i.

Fator referente ao turbinamento na k-ésima equação de função de produçãoda usina hidroelétrica i.

Fator referente ao vertimento na k-ésima equação de função de produçãoda usina hidroelétrica i.

Fator referente ao volume armazenado na k-ésima equação de função deprodução da usina hidroelétrica i.

Taxa de consumo de energia da j-ésima usina elevatória, em relação àvazão bombeada.

Função de custo de geração para a i-ésima unidade térmica, no estágio t,modelada como uma função linear de )(igtt .

Variação horária máxima negativa para a vazão turbinada da usina i, noestágio t.

Variação horária máxima positiva para a vazão turbinada da usina i, noestágio t.

Limite superior de variação horária na vazão da Régua 11 (positiva ounegativa).

Limite superior de variação horária na cota da Régua 11 (positiva ounegativa).

Limite superior de variação diária na vazão da Régua 11 (positiva ounegativa).

Limite superior de variação diária na cota da Régua 11 (positiva ounegativa).

Demanda de energia no k-ésimo subsistema, no estágio t.

Defluência máxima para a i-ésima usina hidroelétrica, no estágio t.

Defluência mínima para a i-ésima usina hidroelétrica, no estágio t.

Demanda de energia na k-ésima ilha da rede elétrica, no estágio t.

Tα

)(iAt

)(•itCterm

)( jCcontrt

)(kCdefct

)( jCons

)( jBombt

)(kDt

)(min iDeft

)(max iDeft

)(kDilht

hQR11∆

hHR11∆

dQR11∆

dQR11∆

)(sup iQt

∆

)(inf iQ t∆

)(iC kI

)(iC kQ

)(iC kS

)(iC kV



Desvio de água da i-ésima usina hidroelétrica, no estágio t.

Máximo desvio de água da i-ésima usina hidroelétrica, no estágio t.

Mínimo desvio de água da i-ésima usina hidroelétrica, no estágio t.

Energia contratada no j-ésimo contrato de importação/exportação deenergia, no estágio t.

Limite inferior de energia contratada no j-ésimo contrato de importação /exportação de energia, no estágio t.

Limite superior de energia contratada no j-ésimo contrato de importação /exportação de energia, no estágio t.

Déficit de energia no k-ésimo subsistema, no estágio t.

),( rkft Intercâmbio de energia do subsistema k para o subsistema r, no estágio t.

),( rkft Limite superior de intercâmbio de energia do subsistema k para osubsistema r, no estágio t.

Função que calcula a evaporação no reservatório da i-ésima usinahidroelétrica, em função de seu volume armazenado.

Função que determina a vazão no canal Pereira Barreto, entre as usinas deIlha Solteira e Três Irmãos, em função dos níveis em ambos reservatórios.

Função que determina a vazão máxima no canal Pereira Barreto, entre asusinas de Ilha Solteira e Três Irmãos, em função dos níveis em ambosreservatórios.

Função de produção energética da i-ésima usina hidroelétrica, expressa emtermos de seu volume médio, vazão turbinada e, em alguns casos, da vazãovertida.

Função que determina a variação de vazão da régua 11, em relação àvariação da cota no mesmo local.

Geração da j-ésima bacia especial, no estágio t.

Geração da i-ésima usina hidroelétrica, no estágio t.

Geração máxima da i-ésima usina hidroelétrica, no estágio t.

Geração em 50 Hz da usina hidroelétrica de Itaipu, no estágio t.

Geração em 60 Hz da usina hidroelétrica de Itaipu, no estágio t.

Geração da j-ésima pequena usina, no estágio t.

Geração da i-ésima unidade térmica, no estágio t.

Limites inferior e superior de geração para a i-ésima unidade térmica, noestágio t.

Geração da i-ésima usina térmica, no estágio t.

Limites inferior e superior de geração para a i-ésima usina térmica, noestágio t.

Ifrom(l) Índice da barra de origem do circuito l.

Ito (l) Índice da barra de destino do circuito l.

)( jEcontrt

)(iDvt

)(iGTt

)(kEdefct

)(igtt

)(max iDvt

)(min iDvt

)( jGBEt

)(50_ ItGH t

)(60_ ItGH t

)( jGPQt

)(•evapf

)(•−TRIILSf

)( jEcontrt

)( jEcontrt

)(11 •Rf

)(•iprodf

)(•−TRIILSf

)(iGH t

)( e )( igtigt tt

)( e )( iGTiGT tt

)(iGH t



ILS Índice da usina de Ilha Solteira.

It Índice da usina de Itaipu.

iusR11 Índice da usina mais a jusante no Rio Iguaçu, na configuração em estudo.

IV, SE e S Índices correspondentes aos nós de Ivaiporã, Sudeste e Sul, narepresentação dos subsistemas no modelo.

Fator de participação da energia consumida na j-ésima unidade elevatória,para o fluxo referente à k-ésima violação de limite de fluxo ocorrida emdeterminado estágio.

Fator de participação da energia contratada no j-ésima contrato, para ofluxo referente à k-ésima violação de limite de fluxo ocorrida emdeterminado estágio.

Fator de participação da geração da i-ésima unidade hidroelétrica, para ofluxo referente à k-ésima violação de limite de fluxo ocorrida emdeterminado estágio.

Fator de participação da geração da i-ésima unidade térmica, para o fluxoreferente à k-ésima violação de limite de fluxo ocorrida em determinadoestágio.

Fator de participação da i-ésima usina hidroelétrica, na k-ésima restriçãoelétrica, no estágio t.

Fator de participação da i-ésima usina térmica, na k-ésima restriçãoelétrica, no estágio t.

Fator de participação da defluência da usina mais à jusante no Rio Iguaçu,na composição da vazão na régua 11.

Fator de participação da defluência de Itaipu, na composição da vazão narégua 11.

Limite inferior para a k-ésima restrição elétrica, no estágio t.

Limite inferior para a k-ésima restrição de somatório de fluxos noscircuitos, no estágio t.

Limite superior para a k-ésima restrição de violação de limite de fluxo, noestágio t.

Limite superior para a k-ésima restrição elétrica, no estágio t.

Limite superior para a k-ésima restrição de somatório de fluxos noscircuitos, no estágio t.

Conjunto das usinas que se encontram imediatamente à montante da usinahidroelétrica i, com tempo de viagem da água entre as usinas MENOR doque a duração do estágio t.

Conjunto de usinas com canal de desvio a jusante para a usinahidroelétrica i.

Conjunto das usinas que se encontram imediatamente à montante da usinahidroelétrica i, com tempo de viagem da água entre as usinas MAIOR doque a duração do estágio t.

)(kLIREt

)(iM t

)(iMtvt

),( kiHtREκ

),( kiTtREκ

IgR

k11

ItR

k11

),( kiK Hf

),( kiK Tf

),( kjK Bf

),( kjK Cf

)(kLSFLt

)(iM Dv

)(kLStSMFLX

)(kLSREt

)(kLItSMFLX



)(iQt

NBE(k) Número de bacias especiais no k-ésimo subsistema, barra ou ilha da redeelétrica.

NB Número de barras na rede elétrica

NC Número de circuitos na rede elétrica

NCIE[(k)] Número total de contratos de importação / exportação de energia. Quandoseguido por um índice k, representa o número de contratos referentes ao k-ésimo subsistema, barra ou ilha da rede elétrica.

NCUTDC Número de inequações (cortes de Benders) que compõem a função decusto futuro do estágio T.

NILH Número total de ilhas da rede elétrica no problema.

NLviolt Número de circuitos violados durante a resolução do subproblema doestágio t.

NPQ(k) Número de pequenas usinas no k-ésimo subsistema ou ilha da rede elétrica.

Número de restrições elétricas no problema.

Número de restrições de somatório de fluxos no problema, para o estágio t.

NSIST Número total de subsistemas.

NUNE[(k)] Número total de unidades elevatórias no problema. Quando seguido porum índice k, representa o número de unidades elevatórias referentes ao k-ésimo subsistema, barra ou ilha da rede elétrica.

NUNT[(k)] Número total de unidades térmicas no problema. Quando seguido por umíndice k, representa o número de unidades referentes à k-ésima usinatérmica, subsistema, barra ou ilha da rede elétrica.

NUSHR(k) Número de usinas hidroelétricas na k-ésima restrição elétrica.

NUSIHDv Número total de usinas hidroelétricas com canal de desvio no problema.

NUH Número total de usinas hidroelétricas no problema.

NUT Número total de usinas térmicas no problema.

NUSTR(k) Número de usinas térmicas na k-ésima restrição elétrica.

Conjunto de subsistemas diretamente conectados ao subsistema k.

Fator referente ao armazenamento da i-ésima usina hidroelétrica, na k-ésima inequação da função de custo futuro do estágio T.

Ângulo de tensão da k-ésima barra da rede elétrica, na solução do estágio t.

Turbinamento da i-ésima usina hidroelétrica, no estágio t.

Turbinamento máximo da i-ésima unidade hidroelétrica, no estágio t.

Limite superior para o turbinamento da i-ésima usina hidroelétrica, noestágio t.

Limite inferior para o turbinamento da i-ésima usina hidroelétrica, noestágio t.

Vazão no canal Pereira Barreto entre as usinas de Ilha Solteira e TrêsIrmãos, no estágio t.

NRE

)(iQt

tTRIILSQ −

kΩ

tkθ

)(ikTVπ

tNRSFL

)(min iQt

)(max iQt



Vazão na Régua 11, no estágio t.

Vazão na Régua 11 k horas anteriores ao estágio t.

Reserva de potência para a i-ésima usina hidroelétrica, no estágio t.

Reserva de potência para a i-ésima usina térmica, no estágio t.

Vertimento da i-ésima usina hidroelétrica, no estágio t.

Limite superior para o vertimento da i-ésima usina hidroelétrica, noestágio t.

Limite inferior para o vertimento da i-ésima usina hidroelétrica, no estágiot.

Último estágio do período de estudo do DESSEM.

Termo independente da k-ésima inequação da função de custo futuro doestágio T.

TRI Índice da usina de Três Irmãos.

Tempo de viagem da água da usina hidroelétrica j até a próxima usina ajusante.

Volume útil armazenado no reservatório da i-ésima usina hidroelétrica, aofinal do estágio t.

Volume útil máximo do reservatório da i-ésima usina hidroelétrica.

Captação de água para usos alternativos, no reservatório da i-ésima usinahidroelétrica, no estágio t.

Volume referente à crista do vertedouro, para o reservatório da i-ésimausina hidroelétrica.

Vazão para enchimento de volume morto do reservatório da i-ésima usinahidroelétrica, no estágio t.

Volume útil de espera, para controle de cheias, do reservatório da i-ésimausina hidroelétrica, no estágio t (porção cheia do reservatório).

Volume evaporado do reservatório da i-ésima usina hidroelétrica, noestágio t.

Limite inferior para a operação do reservatório da i-ésima usinahidroelétrica, no estágio t.

Limite superior para a operação do reservatório da i-ésima usinahidroelétrica, no estágio t.

Reatância do l-ésimo circuito na rede elétrica.

)( jtv

)(iVt

)(iV

)(iSt

kTFCFTI

)(iV tespera

)(iVtsup

)(min iSt

)(max iSt

)(iRH t

)(iRTt

ktRQ − 11

lx

tRQ 11

T

)(iV taltern

)(iVcrist

)(iV tenchim

)(iV tevap

)(iV tinf



Anexo II - Programação Dinâmica Dual Determinística (PDDD)

O conceito de programação dinâmica dual será ilustrado considerando o problema deoperação hidrotérmica (1) de dois estágios em que a afluência a cada usina hidroelétrica, emqualquer estágio do período de planejamento, é conhecida. Este problema pode serrepresentado por:

min/

c x c xs a

A x bE x A x b

1 1 2 2

1 1 1

1 1 2 2 2

+

≥+ ≥

(1)

Este problema pode ser interpretado como um processo de decisão em dois estágios: no 1°estágio, escolhe-se uma decisão x1 viável, x1

* , tal que A x b1 1 1* ≥ . Dada a decisão x1

* , resolve-se o problema de otimização do 2° estágio:

min/

*

c xs a

A x b E x

2 2

2 2 2 1 1≥ −(2)

Observe que x1*é um valor conhecido e passa para o lado direito do conjunto de restrições do

problema. Os vetores x1 e x2 representam os volumes finais das usinas hidroelétricas, asvazões turbinadas, as gerações térmicas, etc., no primeiro e segundo estágio respectivamente.O objetivo do processo é minimizar a soma dos custos do primeiro e segundo estágiosc x c x1 1 2 2+ . O custo do segundo estágio, c x2 2

* , onde x2* é a solução ótima do problema de 2°

estágio, pode ser visto como uma função da decisão x1 , do problema de 1° estágio, isto é:α1 1 2 2

2 2 2 1 1

( ) min/

x c xs a

A x b E x

=

≥ −(3)

O problema de dois estágios (1) pode, portanto, ser reescrito como:min ( )

/c x x

s aA x b

1 1 1 1

1 1 1

+

≥

α(4)

onde α1 1( )x é o valor da solução ótima de (3) para cada x1 . A estrutura do problema (4) ésemelhante à recursão de programação dinâmica, onde c x1 1 representa o custo imediato eα1 1( )x fornece informações acerca das consequências futuras (2° estágio) da decisão x1 .

O princípio de decomposição de Benders 9 é uma técnica que permite construir aproximaçõespara a função α1 1( )x , baseadas na solução iterativa dos problemas de 1 e 2° estágios, comodescrito a seguir.

O comportamento da função custo futuro, α1 1( )x , pode ser caracterizado a partir do dual doproblema de 2° estágio, sob a hipótese de linearidade deste problema. Seja o dual doproblema (2):



( )α π

π

1 1 2 1 1

2 2

( ) max/

x b E xs a

A c

= −

≤(5)

onde o vetor π representa as variáveis duais. Note que na formulação dual a decisão x1 estána função objetivo. Consequentemente, o conjunto de restrições π A c2 2≤ define uma regiãoviável para o problema (5) que não depende da decisão de 1° estágio x1 e pode ser

caracterizada pelo conjunto de todos os pontos extremos ou vértices π π π π= 1 2, , ,K P .

Como a solução ótima de um problema de programação linear sempre corresponde a umvértice da região viável, o problema (5) pode, em princípio, ser resolvido por enumeração:

( )

α π

π π π π

1 1 2 1 1

1 2

( ) max/

, , ,

x b E xs a

i

i P

= −

∈ K

(6)

O problema (6) pode ser rescrito da seguinte forma:

( )( )

( )

α α

α πα π

α π

1 1

12 1 1

22 1 1

2 1 1

( ) min/

xs a

b E xb E x

b E xP

=

≥ −≥ −

≥ −M

(7)

sendo α uma variável escalar.

Como α é maior do que ou igual a cada ( )π i b E x2 1 1− para todo i P= 1, ,K , α é, em

particular, maior do que ou igual ao ( ) max , , ,π i b E x i P2 1 1 1− = K . Como a função objetivo

do problema (7) é minimizar α, esta restrição será atendida na igualdade. Portanto, a soluçãodeste problema é igual a solução ótima do problema (6) e, consequentemente, igual a soluçãoótima do problema (5).

Este problema possui uma interpretação geométrica interessante: a função de custo futuroα1 1( )x é uma função linear por partes da decisão x1 , cujos componentes são hiperplanosdefinidos por ( )π i b E x2 1 1− . Isto implica que a função de custo futuro pode ser caracterizada

sem discretização dos valores x1 , sendo suficiente o conhecimento dos coeficientes π i doshiperplanos suporte.

Sabe-se da teoria de programação linear que o valor da função objetivo do problema primal edo problema dual coincidem na solução ótima. Como o problema (7) é equivalente aoproblema (5), pode-se concluir que as restrições α π≥ −i b E x( )2 1 1 do problema (7) definema função de custo futuro α1 1( )x do problema original (4). O problema de dois estágios (1)pode então ser rescrito como:



min/

( )( )

( )

c xs a

A x bb E xb E x

b E xP

1 1

1 1 11

2 1 12

2 1 1

2 1 1

00

0

+

≥− − ≤− − ≤

− − ≤

α

π απ α

π αM M M M M

(8)

Desta forma, o problema original de decisão de dois estágios (1) pode ser escrito somente emfunção das variáveis do problema de 1° estágio mais a variável escalar α.

Apesar do conjunto de restrições ( ) π i b E x i P2 1 1 1− =, , ,K , ser muito grande, somentealgumas delas estarão ativas na solução ótima. Isto sugere o uso de técnicas de relaxação, naqual as restrições são sucessivamente adicionadas ao problema a cada iteração, fornecendouma aproximação mais precisa para a função custo futuro até alcançar uma tolerânciapreestabelecida.

II.1 Algoritmo de soluçãoO algoritmo programação dinâmica dual (PDD) para o problema de dois estágios é descritopelos seguintes passos:a Inicialize o contador de iterações J = 0 ; o limite superior para a solução ótima z = +∞ ;

a aproximação inicial da função custo futuro $ ( )α 1 1 0x = , para todo x1 (isto significaque não está disponível nenhuma informação sobre o conjunto de pontos extremos ouvértices π) e um valor para a tolerância de convergência TOL.

b Resolva o problema relaxadomin $

/

( ) $ , ,

c xs a

A x bb E x j Jj

1 1

1 1 1

2 1 1 0 1

+

≥− − ≤ =

α

π α K

(9)

c Seja ( )x1* *, $α a solução ótima do problema (9). Pode-se ver que:

z c x= +1 1* *$α (10)

é um limite inferior para a solução do problema original (4), pois o problema (9) é umaversão relaxada do problema (8).

d Dada a decisão x∗1, resolva o problema de 2° estágio:

α 1 1 2 2

2 2 2 1 1

( ) min/

*

*

x c xs a

A x b E x

=

≥ −(11)

e Seja x∗2 a solução ótima do problema (11). O par (x∗

1, x∗2) é uma solução viável do

problema (4), mas não necessariamente a solução ótima. Um limite superior pode sercalculado por:

z z c x c x= +min , * *1 1 2 2 (12)



f Se (z z− ≤ TOL) o problema está resolvido e a solução ótima corresponde ao par

( )x x1 2* *, associado ao limite superior

_z. Caso contrário, execute os procedimentos de (g)

a (h) para melhorar a aproximação da função de custo futuro.g Seja π* o vetor de multiplicadores simplex associados às restrições do problema (11).

Sabe-se da teoria de programação linear que este vetor é uma solução básica viável doproblema dual (5), e portanto um vértice da região viável π A c2 2≤ . Este vértice podeser usado para formar uma nova restrição do tipo π α*( )b E x2 1 1 0− − ≤ , denominadacorte de Benders, que será adicionada ao problema relaxado (9).

h Faça J J J= + =1, *π π e vá para (b).

Os passos a) – h) descrevem o algoritmo de decomposição de Benders. Maiores detalhespodem ser encontrados em [12], [13] e [14]. Na prática, é usual reescrever o corte de Bendersπ α*( )b E x2 1 1 0− − ≤ descrito no passo (g) da seguinte forma: seja w* o valor da soluçãoótima do problema (11) e π* o vetor de multiplicadores simplex associado. Da igualdade desoluções ótimas dos problemas primal e dual podemos escrever:

w b E x* * *( )= −π 2 1 1 (13)

Colocando (π* b2) em evidência, obtém-se:

π π* * * *b w E x2 1 1= +

Substituindo na expressão π α*( )b E x2 1 1 0− − ≤ , obtém-se uma expressão alternativa para ocorte de Benders:

w E x x* * *( )+ − − ≤π α1 1 1 0 (14)O algoritmo PDD pode ser estendido para problemas multi-estágios. Seja o problema dedecisões sequenciais com T estágios:

min/

c x c x c xs a

A x bE x A x b

E x A x b

E x A x b

T T

T T T T T

1 1 2 2

1 1 1

1 1 2 2 2

2 2 3 3 3

1 1

+ + +

≥+ ≥

+ ≥

+ ≥− −

L

M

(15)

Este problema pode ser representado por:min ( )

/c x x

s aA x b

1 1 1 1

1 1 1

+

≥

α(16)

onde α1 1( )x representa as consequências da decisão de 1° estágio, x1 , nas decisões dosdemais estágios. Esta função é calculada por:



TTTTT

TT

bxAxE

bxAxExEbxA

asxcxcx

≥+

≥+−≥

++=

−− 11

33322

11222

2211

/min)(

M

Lα

(17)

Repetindo este procedimento (T-2) vezes obtém-se:α αT T T T T T

T T T T T

x c x xs a

A x b E x

− − − − − −

− − − − −

= +

≥ −

2 2 1 1 1 1

1 1 1 2 2

( ) min ( )/ (18)

onde α T Tx− −1 1( ) é a função do T-ésimo estágio:α T T T T

T T T T T

x c xs a

A x b E x

− −

− −

=

≥ −

1 1

1 1

( ) min/ (19)

O algoritmo de solução do problema multi-estágios é:a Inicialize o contador de iterações J = 0 ; o limite superior para a solução ótima z = +∞ ;

a aproximação inicial da função custo futuro $ ( )α t tx = 0 , para todo xt e para todot T= 1, ,K (isto significa que não está disponível nenhuma informação sobre o conjuntode pontos extremos ou vértices π associados a cada estágio) e um valor para a tolerânciade convergência TOL.

b Resolva o problema aproximado de 1° estágio:min $

/

( ) $ , ,

c xs a

A x bb E x j Jj

1 1 1

1 1 1

2 2 1 1 1 0 1

+

≥− − ≤ =

α

π α K

(20)

cuja solução ótima: ( , $ )* *x1 1α

c calcule z c x= +1 1 1* *$α

d Repita para t = 2,…,T (simulação “forward”)Dado *

1−tx , resolva o problema aproximado do t-ésimo estágio:$ ( ) min $

/

( ) $

( ) $

*

*

α α

π απ α

t t t t t

t t t t t

tj

t t t t

tj

tj

t t t t

x c xs a

A x b E xb E x

ou w E x x

− −

− −

+ +

+ +

= +

≥ −− − ≤

+ − − ≤

1 1

1 1

1 1

1 1

00

(21)

As equações A x b E xt t t t t≥ − − −1 1* representam as restrições do estágio t e as equações

π αtj

t t t tb E x j J+ + − − ≤ =1 1 0 1( ) $ , ,K representam a aproximação da função custofuturo $ ( )α t tx exceto para t = T, onde $ ( )α T Tx é sempre igual a zero.A solução ótima deste problema é ( , $ )* *xt tα



e O vetor ( , , )* *x xT1 K é uma solução viável do problema (15), mas não necessariamente asolução ótima. Portanto,

z z c xt tt

T

=

=

∑min , *

1

é um limite superior da solução ótima.f Se (z z− ≤ TOL) o problema está resolvido e a solução ótima corresponde ao vetor

( , , )* *x xT1 K associado ao limite superior _z. Caso contrário, vá para (g).

g Faça J = J + 1.Repita para t = T,…,2 (recursão “backward”)

Resolva o problema de otimização:

0ˆ)(0ˆ)(

/ˆmin

*11

11

*11

≤−−+≤−−

−≥

+

++

++

−−

ttttj

tj

t

ttttj

t

ttttt

ttt

xxEwouxEb

xEbxAas

xc

απαπ

α

(22)

As equações π αtj

t t t tb E x j J+ + − − ≤ =1 1 0 1( ) $ , ,K representam a aproximação dafunção custo futuro $ ( )α t tx exceto para t = T, onde $ ( )α T Tx é sempre igual a zero.Seja j

tπ o vetor de multiplicadores simplex associado ao conjunto de restrições doproblema (22) na solução ótima. j

tπ medem a variação do custo de operação doestágio t até o final do período de planejamento T devido a variações marginais nosníveis de armazenamento dos reservatórios no início do estágio t (ou final do estágio( )t − 1 ), representados por xt−1

* . Estes multiplicadores são usados para formar umanova restrição do tipo π αt

jt t t tb E x( ) $− − ≤− − −1 1 1 0 (corte de Benders) que será

adicionada à função $ ( )α t tx− −1 1 , obtendo-se uma nova aproximação.h Vá para (b).

O algoritmo de programação dinâmica dual acima descrito apresenta diversas vantagens: nãonecessita discretizar o espaço de estados; fornece a cada iteração limites superior e inferiorpara a solução ótima, permitindo um balanço entre o esforço computacional e a precisãodesejada; o problema de otimização é decomposto numa sucessão de subproblemas de umúnico estágio, podendo tomar partido de algoritmos especializados.

II.2 Consideração do tempo de viagemPara alguns aproveitamentos hidrelétricos não se pode considerar que todo o volumedefluente das usinas de montante durante o estágio t esteja disponível no mesmo estágio.Neste caso, o volume afluente à usina pode ser obtido através da propagação das defluênciasde montante. Uma maneira simplificada de realizar esta propagação consiste em considerar otempo de viagem – intervalo de tempo médio que as defluências consomem até alcançar ausina a jusante – empregando um fator de proporcionalidade ao volume defluente. A Figura4.2 abaixo ilustra este conceito para um sistema composto por duas usinas hidrelétricas emcascata.

Pode-se observar neste exemplo que as afluências que alcançam o reservatório de jusante noestágio t foram defluídas do reservatório de montante nos estágios t, t-1, t-2 e t-3, com fatoresde proporcionalidade iguais a k k k kt t t t

01

12

23

3, , e − − − respectivamente.



Considerar o tempo de viagem como um fator fixo de propagação dos volumes permite queeste aspecto seja incorporado ao modelo DECOMP, resultando em uma modelagem maisprecisa para aqueles sistemas hidrotérmicos em que o tempo de viagem não pode serdesprezado, como por exemplo, a bacia do rio São Francisco.Para considerar estes aspectos, a formulação do modelo deve ser estendida, como serámostrado a seguir.

Figura II.1 - Representação do tempo de viagem

Considere-se o problema multi-estágio abaixo, no qual o estado do sistema está definidoapenas em função das variáveis de um estágio:

min/

c x c x c xs a

A x bE x A x b

E x A x b

E x A x b

T T

T T T T T

1 1 2 2

1 1 1

1 1 2 2 2

2 2 3 3 3

1 1

+ + +

≥+ ≥

+ ≥

+ ≥− −

L

M

Como foi visto no item 4.1, a solução deste problema pode ser obtida através de um algoritmode PDDD. Um dos passos fundamentais deste algoritmo consiste no refinamento daaproximação $α t da função de custo futuro α t tx( ) . Esta função $α t é dada por:

$ ( ) , ,α πt tj

t t tb E x j J≥ − =+ +1 1 1K

ou, equivalentemente$ ( ) , ,*α πt t

jtj

t tj

tw E x x j J≥ + − =+ +1 1 1K

onde wtj+1 é a solução do problema do estágio (t+1) na iteração j e πt

j+1 as variáveis duais

associadas. Cabe ressaltar que, para cada iteração j, a aproximação da função de custo futuro

t-2t-3 t-1t-4 t

kt0kt−1

1kt−22kt−3

3



$α t incorpora um novo segmento linear que depende apenas da solução xtj* e que é expressa

apenas em função de xt, vetor que define o estado do sistema no estágio t.

Por outro lado, quando o problema é tal que o vetor de estado no estágio t depende dediversos estágios anteriores, a estrutura do corte de Benders se altera. Por exemplo,considerando-se um problema de 6 estágios no qual o estado depende de três estágios, aformulação fica:

min/

c x c x c x c x c x c xs a

A x bE x A x bE x E x A x bE x E x E x A x b

E x E x E x A x bE x E x E x A x b

1 1 2 2 3 3 4 4 5 5 6 6

1 1 1

11

1 2 2 2

12

1 21

2 3 3 3

13

1 22

2 31

3 4 4 4

23

2 32

3 41

4 5 5 5

33

3 42

4 51

5 6 6 6

+ + + + +

≥+ ≥+ + ≥+ + + ≥

+ + + ≥+ + + ≥

Esta formulação pode ser identificada com o caso em que o tempo de viagem entre osaproveitamentos não pode ser desprezado. Empregando-se um raciocínio análogo ao que foidesenvolvido no item anterior, este problema pode ser decomposto numa sucessão desubproblemas de um único estágio. Para o estágio t = 6, temos:

w c xs a A x b E x E x E x

6 6 6

6 6 6 51

5 42

4 33

3

=≥ − − −

min/

Considerando-se π 6j o vetor de multiplicadores simplex associados às restrições e w j

6 o valorda solução ótima do problema para xt

j , t = 3, 4 e 5, podemos escrever:w b E x E x E xj j j j j

6 6 6 51

5 42

4 33

3= + − −π ( )e, portanto, o corte de Benders para o estágio 5 é:

[ ]α π5 6 6 51

5 5 42

4 4 33

3 3≥ + − + − + −w E x x E x x E x xj j j j j( ) ( ) ( )

O problema de um estágio resultante para o estágio 5 fica, então:

[ ]

w c xs a

A x b E x E x E xw E x x E x x E x xj j j j j

5 5 5 5

5 5 5 41

4 32

3 23

2

5 6 6 51

5 5 42

4 4 33

3 3

= +

≥ − − −

≥ + − + − + −

min/

( ) ( ) ( )

α

α π

onde o vetor x5 é a variável de decisão deste estágio. Podemos escrever, portanto:

[ ]

w c xs a

A x b E x E x E x aE x w E x E x x E x x bj j j j j j

5 5 5 5

5 5 5 41

4 32

3 23

2

6 51

5 5 6 6 51

5 42

4 4 33

3 3

= +

≥ − − −

+ ≥ + + − + −

min/

( )( ) ( ) ( )

α

π α π

Sendo π 5j e λ5

j os vetores de multiplicadores simplex associados às restrições (a) e (b) e w j5 o

valor da solução ótima do problema para xtj , t = 2, 3 e 4, temos:



[ ][ ]

w b E x E x E x

w E x

j j j j j

j j j j5 5 5 4

14 3

23 2

32

5 6 6 51

5

= − + + +

+

π

λ π

Desta forma, o corte de Benders para o estágio 4 fica:

[ ][ ]

α π

λ π

4 5 5 41

4 4 32

3 3 23

2 2

5 6 42

4 4 33

3 3

≥ + − + − + −

+ − + −

w E x x E x x E x x

E x x E x x

j j j j j

j j j j

( ) ( ) ( )

( ) ( )

Assim, para o estágio 4, a formulação do problema é:

[ ][ ]

w c xs a

A x b E x E x E x aE E x w E x E x x E x x b

E x E x x

j j j j j j j j

j j j j

4 4 4 4

4 4 4 31

3 22

2 13

1

5 41

5 6 42

4 4 5 5 41

4 32

3 3 23

2 2

5 6 42

4 33

3 3

= +

≥ − − −

+ + ≥ + + − + −

+ + −

min/

( )( ) ( ) ( ) ( )

( )

α

π λ π α π

λ π

Sendo π 4j e λ4

j os vetores de multiplicadores simplex associados às restrições (a) e (b) e w j4

valor da solução ótima do problema para xtj , t = 1, 2 e 3, pode ser escrito:

[ ][ ]

w b E x E x E xw E x

E x

j j j j j

j j j j

j j j j

4 4 4 31

3 22

2 13

1

4 5 5 41

4

4 5 6 42

4

= − + +

+ ++

π

λ πλ λ π

O corte de Benders para o estágio 3 fica:

[ ][ ]

[ ]

α π

λ π

λ λ π

3 4 4 31

3 3 22

2 2 13

1 1

4 5 32

3 3 23

2 2

4 5 6 33

3 3

≥ + − + − + −

+ − + −

+ −

w E x x E x x E x x

E x x E x x

E x x

j j j j j

j j j j

j j j j

( ) ( ) ( )

( ) ( )

( )

que, como pode ser deduzido, é a formulação geral do corte para este problema em qualquerestágio. Ressalta-se que a atualização da função de custo futuro leva em conta os efeitos dosestados definidos em estágios anteriores ao estágio corrente, implicando em modificação nafórmula geral de atualização da função de custo futuro.Desta forma, a formulação geral para o corte de Benders para qualquer estágio para umproblema cujo estado é definido por L estágios anteriores é dada por:

[ ]

[ ]∑ ∑∏

∑−

=

−

=−+−+

+−+++

=+

=−+−+−+++

−

+

−+≥

1

1 11111

1

111111

)(

)(

L

k

kL

llt

jlt

kllt

jkt

k

i

jit

L

llt

jlt

llt

jt

jtt

xxE

xxEw

πλ

πα

A generalização para problemas estocásticos é feita de maneira análoga àquela apresentada noitem 4.2, considerando-se valores esperados para w, π e λ.

II.2.1 Exemplo com tempo de viagemEste procedimento pode ser ilustrado para o caso em que o sistema hidroelétrico é compostopor apenas duas usinas em cascata, sendo o tempo de viagem entre elas representados pelosfatores k k k kt t t t

01

12

23

3, , e − − − para os lags 0, 1, 2 e 3 respectivamente, como na Figura 4.1.



As variáveis que fazem a ligação entre os estágios são dadas pelos volumes armazenados nofinal do estágio e pelos volumes turbinados e vertidos no estágio. Agrupando estas variáveisno início do vetor de estados xt temos:

x v v q q s st t t t t t t= ( , , , , , , )1 2 1 2 1 2 K

onde:vt

i = volume final armazenado na usina i no final do estágio tqt

i = volume turbinado na usina i no estágio tst

i = volume vertido na usina i no estágio tSem perda de generalidade, podemos escrever:

x v v d dt t t t t= ( , , , , )1 2 1 2 K

onde dti representa o volume total defluente na usina i no estágio t.

Empregando a fórmula geral acima, o corte de Benders para o estágio t é dado por:

[ ][ ]

[ ]

α π

λ π

λ λ π

t tj

tj

t tj

t t tj

t t tj

t

tj

tj

t tj

t t tj

t

tj

tj

tj

t tj

t

w E x x E x x E x x

E x x E x x

E x x

≥ + − + − + −

+ − + −

+ −

+ + − − − − − −

+ + − − −

+ + +

1 11

12

1 1 23

2 2

1 22

13

1 1

1 2 33

( ) ( ) ( )

( ) ( )

( )

Neste caso, as matrizes Etl que representam o acoplamento entre os estágios são escritas

como:

Ekt

t

11

1 0 0 00 1 0

=−

− −

e, para l > 1,

E ktl

tl=

−

0 0 0 00 0 0

onde ktl é o fator no estágio t para o lag l (defasagem).

Fazendo π π πt t= ( , )1 2 temos:

[ ][ ]

[ ]

α ππ

λ π

λ λ π

t tj

tj

tj

t

tj

tj

t t tj

t t tj

t t tj

t

tj

tj

t tj

t t tj

t

tj

tj

tj

t tj

t

w v vv v k d d k d d k d d

k d d k d d

k d d

≥ − −− − + − + − + −

− − + −

− −

+ +

+ − − − − − −

+ + − − −

+ + +

11

11 1

21

2 2 1 1 11

21

11

12

32

12

1

12

22 1 1

13

11

11

1 22

33 1 1

( ) ( )( ) ( ) ( ) ( ) ( )

( ) ( ) ( )

( ) ( )

ou, agrupando os termos:

[ ][ ][ ]

α π ππ λ π λ λ π

π λ π

π

t tj

tj

tj

t tj

tj

t

tj

t tj

tj

t tj

tj

tj

t tj

t

tj

t tj

tj

t tj

t

tj

t tj

t

w v v v vk k k d d

k k d d

k d d

≥ − − − −− + + −

− + −

− −

+ + +

+ + + + + +

+ − + + − − −

+ − − −

11

11 1 2

12 2

21

11

22

21 2

23

3 1 1

21 1

21 2 1

31

11

1

21 2

32

12

1

( ) ( ) ( ) ( )( ) ( ) ( ) ( )

( ) ) ( )

( ) ( )



Anexo III: FORMULAÇÃO MATEMÁTICA COMPLETA DOPROBLEMA DE DESPACHO RESOLVIDO PELOMODELO DESSEM

A formulação do problema com todas as características/restrições de operação do sistema, éapresentada a seguir. Porém, antes descreve-se novamente a natureza de cadacaracterística/restrição identificada por seu índice na formulação do problema.

a) Características/restrições do sistema• Demanda por subsistema (estudos sem rede elétrica), ou

demanda por barra (estudos com rede elétrica) (2)

• Demanda por ilha elétrica (estudo com rede) (3)

• Limites de interligação entre subsistemas (4)

• Restrições elétricas (5)

• Limites de fluxos nos circuitos (34)

• Somatório de fluxos nos circuitos (estudos com rede) (35)

• Restrição de transporte entre Itaipu 50 e 60 Hz e os subsistemas SU/SE(6)

• Custos e limites de energia importada/exportada nos contratos de energiacom fontes externas (1) – (7)

• Consumo nas unidades elevatórias (2) – (3)

• Custos de déficit de energia por subsistema (1)

• Consideração, como abatimento de carga, da geração em bacias especiais, pequenasusinas e outras fontes não incluídas na otimização (2) – (3)

b) Características/restrições específicas das usinas hidroelétricas• Limites físicos de armazenamento nos reservatórios (8)

• Limites de bombeamento das unidades elevatórias (9)

• Restrições de balanço hídrico por usina (10)

• Afluência natural às usinas (10)

• Retirada de água para usos alternativos (irrigação, abastecimento) (10)

• Enchimento de volume morto (10)

• Tempo de viagem da água entre aproveitamentos (10)

• Volume evaporado nos reservatórios (10) – (11)

• Restrições operacionais de limites de armazenamento, turbinamento, vertimento,defluência total e volume desviado nos reservatórios / usinas (12) – (16)

• Função de produção das usinas hidroelétricas (produtividade variável com a queda) (17)

• Potência máxima nas usinas (17) – (18)



• Engolimento máximo das usinas (19)

• Restrições de reserva por usina (20)

• Restrições de vertimento máximo nos reservatórios (21)

• Volume de espera nos reservatórios, para controle de cheias (22)

• Restrições operacionais de limites de variação de vazão turbinadanas usinas (23)

• Limites de variação horária e diária na cota da Régua 11 (24) – (27)

• Limites permitidos de vazão para a operação do canal Pereira Barreto (28) – (29)

• Função de custo futuro ao final do horizonte de estudo, em função do armazenamento nosreservatórios (fornecida pelo modelo DECOMP) (30)

c) Características/restrições específicas das usinas térmicas• Custos por unidade (1)

• Limites de geração mínima e máxima nas usinas e unidades; (31) - (32)

• Restrições de reserva por usina (33)

O problema a ser resolvido toma então a seguinte forma:

( )

T

NSIST

ktt

T

t

NCIE

jtt

NUT

it

it

kEdefct(k)Cdefc

jEcontr(j)CcontriGTCtermMinimize

α )( .

)( . .)(

1

1 11

+

+

+

∑

∑ ∑∑

=

= == (1)

sujeito a :

( )

TtNSIST,k

kDkEdefctjEcontrjGPQjGBE

rkfkrfjBombjCons -iGTiGH

tt

kNCIE

jt

kNPQ

jt

kNBE

jt

rtt

kNUNE

jt

kNUT

it

kNUH

it

k

1,..., , ...,1

, )( )( )( )( )(

),(),( )().( )( )(

)(

1

)(

1

)(

1

)(

1

)(

1

)(

1

==

=+±++

+−++

∑∑∑

∑∑∑∑

===

Ω∈===

(2)



( )

TtNILH,k

kDilhjEcontrjGPQjGBE

rkfkrfjBombjCons-iGTiGH

t

kNCIE

jt

kNPQ

jt

kNBE

jt

rtt

kNUNE

jt

kNUT

it

kNUH

it

k

1,..., , ...,1

, )( )( )( )(

),(),( )().( )( )(

)(

1

)(

1

)(

1

)(

1

)(

1

)(

1

==

=±++

+−++

∑∑∑

∑∑∑∑

===

Ω∈===

(3)

TtrNSISTkrkfrkf ktt ,...,1 , ,,...,1 , ),(),(0 =Ω∈=≤≤ (4)

TtNRE k

kLSREiGTkiiGHkikLIREkNUSTR

itt

TtRE

kNUSHR

it

HtREt

,...,1 ,,...,1

, )( )( ),( )( ),( )()(

1

)(

1

==

≤+≤ ∑∑==

κκ(5)

),(),(),(),()(60_ )(60_ )(50_ )(

IVSfIVSEfSIVfSEIVfItGHItGHItGHItGH

ttttt

ttt

−−+=+=

(6)

TtNCIEjjEcontrjEcontrjEcontr ttt 1,..., ,,...,1 , )( )( )( ==≤≤ (7)

TtNUHiiViVt ,...,1 ,,...,1 , )( )( 0 ==≤≤ (8)

TtNUNEjjBombjBombjBomb ttt ,...,1 ,,...,1 , )( )( )( ==≤≤ (9)

( ) ( )

( ) ( )

,...,1 ,...,1

, )( )( )( )()()(

)()( )()()( )( )( )(

)()()(

)(

)(1

T tNUH, i

iViViVjSjQjDv

jSjQiDviSiQiAiViV

tttiMtvj

jtvtjtvtiMj

t

iMjtttttttt

enchimevapalterntDV

t

==

−−−++

+++++−+=

∑∑

∑

∈−−

∈

∈+

(10)

IV

SE S

It

GH_50t (It)[1] (excluindo

parcela

ft (IV,SE)

ft (SE,IV)ft (IV,S)

ft (S,IV)

GH_60t (It)



,...,1 ,...,1 , ))(( )( 1 T tNUH, i iVfiV tevaptevap === − (11)

TtNUHiiViViVttt ,...,1 ,,...,1 , )( )( )( ==≤≤ supinf (12)

TtNUHiiQiQiQ ttt ,...,1 ,,...,1 , )( )( )( ==≤≤ maxmin (13)

TtNUHiiSiSiS ttt ,...,1 ,,...,1 , )( )( )( ==≤≤ maxmin (14)

TtNUHiiDefiDviSiQiDef ttttt ,...,1 ,,...,1 , )( )()()( )( ==≤++≤ maxmin (15)

TtNUHiiDviDviDv Dvttt ,...,1 ,...,1 , )( )( )( ==≤≤ maxmin (16)

( ) T tNUH,,iiSiQiViVfiGH ttttiprodt ,...,1 ...,1 )(),(),(),( )( 1

=== + ,

onde

iprodf é uma função linear por partes da forma:

)()()()()( )()()( 111

11 iSiCiQiCiViCiCiGH tStQtVIt +++≤ + (17)

)()()()()( )()()( 221

22 iSiCiQiCiViCiCiGH tStQtVIt +++≤ +

M

)()()()()( )()()( 1 iSiCiQiCiViCiCiGH tnSt

nQt

nV

nIt

iiii +++≤ +

TtNUHiiGHiGH tt ,...,1 ,,...,1 , )( )(0 ==≤≤ (18)

TtNUHiiQiQt ,...,1 ,,...,1 , )( )(0 ==≤≤ (19)

TtNUHiiRHiGHiGH ttt 1,..., ,,...,1 , )( )()( ==≥− (20)

( ) ( )

( )

)( )(

)()( )()()( )( )( )(

)(

)()()()(

iVcristjDv

jSjQjSjQiQiAiViS

iMjt

iMjtt

iMtvjjtvtjtvtttt

DV

tt

−+

+++++−+≤

∑

∑∑

∈

∈∈−−

(21)

TtNUHiiiVtt espera ,...,1 ,,...,1 , )(V )( 0 ==≤≤ (22)

TtNUHiiQiQiQiQtttt ,...,1 ,,...,1 )( )()()( sup1inf ==∆≤−≤∆− − (23)

( ) ( ) T tItSItQkiusSiusQkQ ttIt

RtRtIg

tR RR,...,1 , )()( )()(

1111 111111 =+++= (24)



TtQQQQ hRtRtR

hR ,...,1 111111111 =∆≤−≤∆− − (25)

TtkQQQQ dRktRtR

dR ,...,1,24,...,1 - 11111111 ==∆≤≤∆− − (26)

)( 111111h

RRh

R HfQ ∆=∆ (27a)

)( 111111d

RRd

R HfQ ∆=∆ (27b)

T t VVfQ TRIILS ttTRIILStTRIILS ,...,1 ,),( )()( == −− (28)

T t VVfQ TRIILS ttTRIILStTRIILS ,...,1 ,),( )()( =≤ −− (29)

...,1 , )( )( 1

1 NCUTDC, k iViTINDAM

iT

kTV

kTFCFT =+≥ ∑

=+πα (30)

TtNUNTiigtigtigt ttt 1,..., ,,...,1 )()()( ==≤≤ (31)

TtNUTiiGTiGTiGT ttt 1,..., ,,...,1 )()()( ==≤≤ (32)

TtNUTiiRTiGTiGT ttt 1,..., ,,...,1 , )( )()( ==≥− (33)

ou

( ) TtNUTiiRTjgtjgtiNUNT

jttt 1,..., ,,...,1 , )( )()(

)(

1==≥−∑

=

(33)

T,tNLviol,kkLSFLjBombjConskjK

jEcontrkjKigtkiKighkiK

tt

NUNE

jtt

Bf

NCIE

jt

Cf

NUNT

it

Tf

NUNH

it

Hf

...1 ,...,1 , )( )().( ),(-

)( ),( )( ),( )( ),(

1

111

==≤

+++

∑

∑∑∑

=

=== (34)

T,tNRSFL,kkLSjBombjConskjK-

jEcontrkjKigtkiKighkiKkLI

ttSMFLX

kNUNE

jtt

Bf

kNCIE

jt

Cf

kNUNT

it

Tf

kNUNH

it

HftSMFLX

...1 ,...,1 , )( )().( ),(

)( ),()( ),( )( ),( )(

)(

1

)(

1

)(

1

)(

1

==≤

+++≤

∑

∑∑∑

=

=== (35)

Documents

Manual Referencia Dessem GTRD - simsee.org · Ressalta-se que cada um destes modelos necessita, para sua utilização, de uma série de modelos auxiliares para fornecer os dados necessários