Embed Size (px)
Citation preview
Métodos de Análise da Estabilidade Transitória de
Sistemas de Energia Eléctrica
João Pedro de Carvalho Mateus
Dissertação para obtenção do Grau de Mestre em
Engenharia Electrotécnica e de Computadores
Júri
Presidente: Prof. Paulo José da Costa Branco
Orientador: Prof. José Pedro da Silva Sucena Paiva
Co-orientador: Prof. Maria Eduarda de Sampaio Pinto de Almeida Pedro
Vogais: Prof. José Manuel Dias Ferreira de Jesus
Outubro 2010
III
Agradecimentos
Há várias pessoas sem as quais a realização deste trabalho não teria sido possível.
Começo por agradecer ao meu orientador, Prof. José Pedro Sucena Paiva, por me ter
dado a oportunidade de trabalhar neste tema, cuja curiosidade e interesse nasceram nas
cadeiras por ele leccionadas. À minha co-orientadora, Prof. Maria Eduarda de Almeida Pedro, um
muito obrigado por toda a disponibilidade e prontidão sempre demonstradas.
Aos meus colegas e amigos, com quem partilhei momentos de trabalho e divertimento
durante todo o curso, mas também aos que me acompanharam durante todo o meu percurso
escolar, obrigado a todos por terem de alguma forma feito parte do meu caminho.
A toda a minha família, mas principalmente aos meus pais, agradeço por todo o carinho
e apoio, todas as condições que me proporcionaram, sem as quais não teria sido possível. Ao
meu irmão, obrigado pelos momentos de descontracção. Aos meus primos Joana, Hugo e
Pedrocas, obrigado pela hospitalidade e disponibilidade que sempre me ofereceram.
Por fim, e porque os últimos são os primeiros, um obrigado muito especial à minha
namorada. Obrigado por todas as palavras de apoio nos momentos mais difíceis, todos os
momentos partilhados, e pela compreensão sempre demonstrada. Obrigado por fazeres parte da
minha vida.
A todos, um muito obrigado!
João Mateus
V
Resumo
Esta dissertação aborda o tema da estabilidade transitória dos Sistemas de Energia
Eléctrica. Mais concretamente, tem como objectivos apresentar o estado da arte no que toca a
métodos de análise da estabilidade transitória dos SEE e implementar um método híbrido com o
mesmo propósito.
Relativamente à apresentação dos métodos já existentes, são referidas as principais
características dos Métodos de Integração Numérica, Métodos Directos, Métodos Híbridos e
Técnicas de Inteligência Artificial. Segue-se uma comparação entre estes, analisando vantagens
e desvantagens.
Quanto ao método implementado, procurou-se desenvolver um algoritmo capaz de
calcular os tempos críticos de actuação das protecções para diferentes perturbações a ocorrer
numa rede. Este método híbrido conjuga as vantagens dos Métodos Directos com as dos
Métodos de Integração Numérica. Aos primeiros vai buscar a rápida análise enquanto que dos
segundos obtém as possibilidades de modelação. Utiliza índices para avaliar a estabilidade
transitória que permitem interromper o processo de integração das equações antes do tempo
total de simulação ser atingido. Utiliza o critério das áreas iguais para estudar um modelo
reduzido equivalente do sistema, constituído por uma máquina ligada a um barramento de
potência de curto-circuito infinita, para o qual obtém o valor da margem de estabilidade
transitória. Baseado nos tempos de actuação das protecções e no valor da margem de
estabilidade transitória, estima o valor do tempo crítico que se pretende determinar através de
processos de regressão linear. Para além da apresentação do método implementado são
apresentados os resultados dos testes a que este foi sujeito.
Palavras-chave: Estabilidade transitória de Sistemas de Energia Eléctrica, Métodos
híbridos, Métodos de integração numérica, Critério das áreas iguais, Tempo crítico de actuação
das protecções
VII
Abstract
In this thesis it is studied the topic of Electric Power Systems transient stability. More
specifically, it aims to present the state of the art regarding the methods of analysis of transient
stability of EPS and implement a hybrid method with that objective too.
In the presentation of the state of the art, the main characteristics of different methods are
listed. The referred methods are Numerical Integration Methods, Direct Methods, Hybrid Methods
and Artificial Intelligence Techniques. After this, a comparison between them is made, showing
their advantages and disadvantages.
In the development of the hybrid method, the objective was to obtain an algorithm
capable of getting the critical clearing time for different disturbances occurring in an electrical
network. This method combines advantages of direct formulations with those of hybrid
formulations. From the first it gets de fast analysis, from the second it gets the modelling
capabilities. The numerical integration process is interrupted using indexes that allow to know if
the system is stable or not before the total simulation time. An equivalent model of the system,
made of a single machine connected to an infinite bus, is studied using the equal area criterion,
from which it is computed the transient stability margin. Using a linear regression, based on the
values of the transient stability margin and respective clearing time of several simulations, the
critical clearing time is computed. Besides the presentation of the hybrid method, the results of
tests to which the method was subjected are presented.
Keywords: Electric Power System transient stability, Hybrid methods, Equal area
criterion, Numerical integration methods, Critical clearing time
IX
Índice Lista de Figuras ............................................................................................................................ XI
Lista de Tabelas ......................................................................................................................... XIII
Lista de Símbolos e Abreviações ............................................................................................... XV
1. Introdução .................................................................................................................................. 1
1.1. Contextualização ........................................................................................................... 1
1.2. Objectivos e estrutura da dissertação ........................................................................... 2
2. Modelização do Sistema de Energia Eléctrica .......................................................................... 5
2.1. Máquina Síncrona .......................................................................................................... 5
2.2. Linha de Transmissão ................................................................................................... 7
2.3. Transformador ............................................................................................................... 8
2.4. Carga ........................................................................................................................... 10
2.5. Modelo Global da Rede ............................................................................................... 11
2.6. Contingências Simuladas ............................................................................................ 12
2.7. Conclusões .................................................................................................................. 13
3. Métodos de análise da estabilidade transitória do SEE .......................................................... 15
3.1. Métodos de Integração Numérica ............................................................................... 15
3.2. Métodos Directos ......................................................................................................... 17
3.2.1. Método de Lyapunov ........................................................................................... 18
3.2.2. Método da Função de Energia Transitória (TEF) ................................................ 19
3.2.3. Superfície Limite de Energia Potencial (SLEP) ................................................... 21
3.2.4. Método do Ponto de Equilíbrio Instável de Controlo (BCU) ................................ 26
3.2.5. Método das Áreas Iguais ..................................................................................... 27
3.3. Métodos Híbridos ......................................................................................................... 30
3.3.1. Método híbrido baseado na função de energia transitória .................................. 31
3.3.2. Método híbrido baseado no método das áreas iguais ........................................ 31
3.4. Técnicas de Inteligência Artificial ................................................................................ 32
3.4.1. Reconhecimento de Formas ............................................................................... 32
3.4.2. Redes Neuronais ................................................................................................. 33
3.4.3. Árvores de Decisão ............................................................................................. 33
3.5. Comparação dos diferentes métodos .......................................................................... 34
X
3.6. Conclusões .................................................................................................................. 35
4. Método híbrido implementado ................................................................................................ 37
4.1. Introdução dos dados .................................................................................................. 37
4.1.1. Barramentos ........................................................................................................ 38
4.1.2. Geradores ........................................................................................................... 38
4.1.3. Linhas de Transmissão ....................................................................................... 38
4.1.4. Transformadores ................................................................................................. 39
4.1.5. Contingências ...................................................................................................... 39
4.2. Cálculo dos valores pré-defeito ................................................................................... 39
4.3. Ciclo para determinação do tempo crítico .................................................................. 40
4.3.1. Índices de estabilidade e instabilidade ............................................................... 41
4.3.2. Identificação das máquinas críticas .................................................................... 43
4.3.3. Redução do sistema a uma máquina equivalente .............................................. 45
4.3.4. Cálculo da margem de estabilidade transitória ................................................... 47
4.3.5. Estimativa dos valores do tempo crítico ............................................................. 49
4.4. Resultados fornecidos ................................................................................................. 50
4.5. Ficheiros de código ..................................................................................................... 50
4.6. Conclusões .................................................................................................................. 51
5. Resultados computacionais .................................................................................................... 53
5.1. Testes computacionais ............................................................................................... 53
5.1.1. Resultados obtidos .............................................................................................. 53
5.1.2. Análise dos erros................................................................................................. 54
5.1.3. Correcções a implementar .................................................................................. 56
5.2. Exemplo de aplicação ................................................................................................. 57
5.3. Conclusões .................................................................................................................. 61
6. Conclusões e Propostas para trabalhos futuros ..................................................................... 63
6.1. Conclusões .................................................................................................................. 63
6.2. Propostas para trabalhos futuros ................................................................................ 65
Bibliografia .................................................................................................................................. 67
Anexo 1 ....................................................................................................................................... 69
XI
Lista de Figuras
Figura 2.1 - Esquema equivalente da máquina síncrona .............................................................. 5
Figura 2.2 - Esquema equivalente em π de uma linha .................................................................. 8
Figura 2.3 - Esquema equivalente simplificado de um transformador .......................................... 8 Figura 2.4 - Esquema equivalente de um transformador com regulação de tensão .................... 9 Figura 2.5 - Esquema equivalente de um transformador desfasador ......................................... 10 Figura 2.6 - Representação esquemática do SEE ...................................................................... 11
Figura 3.1 - Fluxograma para simulação da estabilidade transitória .......................................... 16 Figura 3.2 - Fluxograma do método da bissecção para cálculo do tempo crítico ....................... 17 Figura 3.3 - Rede utilizada no exemplo do Método da SLEP ..................................................... 22
Figura 3.4 - Gráfico de 𝑉𝑃𝐸(𝜃1) .................................................................................................... 23
Figura 3.5 - Gráfico de 𝑉(𝜃1) ....................................................................................................... 24 Figura 3.6 - Superfície Limite de Energia Potencial .................................................................... 25 Figura 3.7 - Representação gráfica do critério das áreas iguais ................................................. 27 Figura 3.8 - Representação gráfica da situação estudada no exemplo do Método das Áreas
Iguais ........................................................................................................................................... 30 Figura 3.9 - Exemplo de uma árvore de decisão ........................................................................ 33
Figura 4.1 - Estrutura geral do método implementado ................................................................ 37 Figura 4.2 - Diagrama do ciclo para o cálculo do tempo crítico .................................................. 40 Figura 4.3 - Comportamento dos índices IDCS e IDE, situação instável .................................... 43 Figura 4.4 - Comportamento dos índices IDCS e IDE, situação estável .................................... 43 Figura 4.5 - Comportamento do índice IDTO .............................................................................. 44 Figura 4.6 - Escolha do método de aproximação da curva de potência a utilizar ...................... 49
Figura 4.7 - Relação linear entre 𝜂𝜂 e 𝑡𝑡𝑐𝑙 ...................................................................................... 49 Figura 4.8 - Interacção dos diferentes ficheiros de código no método híbrido implementado ... 51
Figura 5.1 - Potência eléctrica gerada (-.), aproximação polinomial da potência eléctrica gerada
(linha contínua) e potência mecânica (--) da máquina equivalente para a perturbação 22 ....... 55 Figura 5.2 - Ângulo rotórico e potência eléctrica gerada da máquina equivalente para a
contingência 22 ............................................................................................................................ 56
Figura 5.3 - Ângulos rotóricos para a perturbação 2, com 𝑡𝑡𝑐𝑙 = 0,6 s ......................................... 58
Figura 5.4 - Índices IDCS e IDE para a perturbação 2, com 𝑡𝑡𝑐𝑙 = 0,6 s ...................................... 58 Figura 5.5 - Curvas temporais do ângulo rotórico e potência eléctrica gerada da máquina
equivalente, para a perturbação 2, com 𝑡𝑡𝑐𝑙 = 0,6 s ..................................................................... 59 Figura 5.6 - Curva de potência da maquina equivalente (a tracejado) e aproximação sinusoidal
(a cheio), para a perturbação 2, com 𝑡𝑡𝑐𝑙 = 0,6 s ......................................................................... 59
XII
Figura A1.1 - Rede de teste da CIGRE ....................................................................................... 69
XIII
Lista de Tabelas
Tabela 3.1 - Comparação entre os diferentes métodos .............................................................. 34
Tabela 4.1 - Estrutura do ficheiro com os dados dos barramentos ............................................ 38 Tabela 4.2 - Estrutura do ficheiro com os dados dos geradores ................................................ 38 Tabela 4.3 - Estrutura do ficheiro com os dados das linhas ....................................................... 39 Tabela 4.4 - Estrutura do ficheiro com os dados dos transformadores ...................................... 39 Tabela 4.5 - Estrutura do ficheiro com os dados dos defeitos .................................................... 39 Tabela 4.6 - Estrutura do ficheiro com os resultados da simulação ........................................... 50
Tabela 5.1 - Tempos críticos de actuação das protecções para a rede de teste da CIGRE ...... 54 Tabela 5.2 - Conteúdo do ficheiro de entrada utilizado no exemplo de aplicação ..................... 57 Tabela 5.3 - Identificação do conjunto das máquinas críticas para a perturbação 2, com
𝑡𝑡𝑐𝑙 = 0,6 s ..................................................................................................................................... 59 Tabela 5.4 - Dados relativos às iterações realizadas com a aproximação sidusoidal ................ 60 Tabela 5.5 - Dados relativos às iterações realizadas com a aproximação polinomial ................ 60
Tabela A1.1 - Características das linhas da rede de teste da CIGRE ........................................ 69 Tabela A1.2 - Características dos geradores da rede de teste da CIGRE ................................. 70 Tabela A1.3 - Resultados do trânsito de energia da rede de teste da CIGRE ........................... 70
XV
Lista de Símbolos e Abreviações
𝑬𝑬𝒊𝒊′ - força electromotriz transitória da máquina 𝑑𝑑
𝐸𝑖 - módulo da força electromotriz transitória da máquina 𝑑𝑑
𝑽𝑽𝒊𝒊 - tensão aos terminais da máquina 𝑑𝑑
𝑋𝑋𝑑𝑖′ - reactância transitória da máquina 𝑑𝑑
𝑰𝑰𝒊𝒊 - intensidade de corrente fornecida pela máquina 𝑑𝑑
𝐻𝑖 - constante de inércia da máquina 𝑑𝑑
𝜔0 - velocidade angular nominal do sistema
𝛿𝛿𝑖 - ângulo da força electromotriz transitória da máquina 𝑑𝑑, ângulo rotórico da máquina 𝑑𝑑
𝐷𝐷𝑖 - coeficiente de amortecimento da máquina 𝑑𝑑
𝑃𝑃𝑚𝑖 - potência mecânica da máquina 𝑑𝑑
𝑃𝑃𝑒𝑖 - potência eléctrica gerada pela máquina 𝑑𝑑
𝑃𝑃𝑎𝑐𝑖 - potência de aceleração do gerador 𝑑𝑑
𝜔𝑖 - velocidade angular da máquina 𝑑𝑑
𝛿𝛿0 - ângulo rotórico do centro de inércia do sistema
𝑓0 - frequência nominal do sistema
𝑀𝑇 - coeficiente de inércia total
𝑀𝑖 - coeficiente de inércia da máquina 𝑑𝑑
𝜃𝑖 - ângulo rotórico da máquina 𝑑𝑑, referido ao centro de inércia do sistema
𝜔�𝑖 - velocidade angular da máquina 𝑑𝑑, referida ao centro de inércia do sistema
𝑓𝑖(𝜃𝑖) - potência de aceleração da máquina 𝑑𝑑, referida ao centro de inércia do sistema
𝑃𝑃𝐶𝑂𝐼 - potência de aceleração do sistema
𝑽𝑽𝒆𝒆 - tensão na emissão
𝑽𝑽𝒓𝒓 - tensão na recepção
𝑰𝑰𝒆𝒆 - intensidade de corrente na emissão
𝑰𝑰𝒓𝒓 - intensidade de corrente na recepção
𝑅𝑅𝐿 - resistência longitudinal
𝑋𝑋𝐿 - reactância longitudinal
𝑌𝑌𝑇 - admitância transversal
𝑽𝑽𝒑𝒑 - tensão no primário
𝑽𝑽𝒔𝒔 - tensão no secundário
𝑰𝑰𝒑𝒑 - corrente no primário
𝑰𝑰𝒔𝒔 - corrente no secundário
𝑅𝑅𝑇 - resistência total
𝑋𝑋𝑇 - reactância total
𝑽𝑽𝒔𝒔′ - tensão no secundário do transformador ideal
𝑰𝑰𝒔𝒔′ - corrente no secundário do transformador ideal
𝒁𝑪𝑪 - impedância de curto-circuito
XVI
𝑚𝑚′ - relação de transformação
𝒀𝒀 - matriz de admitâncias
𝐺𝐺 - parte real da matriz de admitâncias
𝐵 - parte imaginária da matriz de admitâncias
𝒀𝒀𝒄𝒄𝒄𝒄 - admitância de curto-circuito
𝒎𝒎′ - relação de transformação complexa
𝒀𝒀𝒄𝒄𝒄𝒄 - admitância da carga ligada ao barramento 𝑘
𝐺𝐺𝑐𝑘 - parte real de 𝒀𝒀𝒄𝒄𝒄𝒄
𝐵𝑐𝑘 - parte imaginária de 𝒀𝒀𝒄𝒄𝒄𝒄
𝑃𝑃𝑐𝑘 - potência activa de carga do barramento 𝑘
𝑄𝑐𝑘 - potência reactiva de carga do barramento 𝑘
𝑉𝑘 - módulo da tensão no barramento 𝑘
[𝑬𝑬′] - vector das forças electromotrizes
[𝑽𝑽] - vector das tensões nos barramentos
[𝑰𝑰] - vector das correntes injectadas nos barramentos de geração
[𝒀𝒀𝒎𝒎𝒎𝒎], [𝒀𝒀𝒎𝒎𝒄𝒄], [𝒀𝒀𝒄𝒄𝒎𝒎], [𝒀𝒀𝒄𝒄𝒄𝒄] - matrizes de admitâncias que descrevem a rede aumentada
𝑚𝑚 - número de geradores, equivalente ao número de barramentos de geração
𝑛 - número de barramentos
𝒀𝒀𝒓𝒓𝒆𝒆𝒅 - matriz de impedâncias reduzida
𝒀𝒀𝒅𝒆𝒆𝒇 - admitância de defeito
𝑡𝑡𝑚𝑎𝑥 - tempo máximo de simulação
𝑡𝑡𝑐𝑟 - tempo crítico de actuação das protecções
𝑡𝑡𝑐𝑙 - tempo de actuação das protecções
𝑡𝑡𝑐𝑙0 - estimativa inicial do tempo crítico de actuação das protecções, correspondente a uma
situação instável
𝑡𝑡𝑖𝑛𝑠𝑡 - tempo de actuação das protecções relativo a uma situação instável
𝑡𝑡𝑒𝑠𝑡 - tempo de actuação das protecções relativo a uma situação estável
∆𝑡𝑡 - passo do processo de integração
𝜀𝜀 - tolerância
𝐸(𝑚𝑚) - energia total do sistema
𝑚𝑚𝑒 - ponto(s) de equilíbrio
𝑉(𝑚𝑚) - função de Lyapunov
𝑉𝐾𝐸 - energia cinética do sistema
𝑉𝑃𝐸 - energia potencial do sistema
𝑉𝑐𝑟 - valor crítico da função de Lyapunov
𝑉𝑐𝑙 - valor da função de Lyapunov no instante de eliminação do defeito
𝑓𝑖𝐷 - potência de aceleração do sistema no período de defeito
𝑓𝑖𝑃 - potência de aceleração do sistema no período pós-defeito
𝜃𝑢 - ponto de equilíbrio instável de controlo
XVII
𝜃𝑠 - ponto de equilíbrio estável do sistema na configuração pós-defeito
𝑉𝑃𝐸𝑚𝑎𝑥 - valor máximo da energia potencial do sistema
𝜃𝑐𝑙 - ângulo rotórico no instante 𝑡𝑡𝑐𝑙, referido ao centro de inércia
𝜔�𝑐𝑙 - velocidade angular no instante 𝑡𝑡𝑐𝑙, referida ao centro de inércia
𝜃0 - ângulo rotórico inicial, referido ao centro de inércia
𝜃𝑎𝑝𝑝𝑢 - ponto de equilíbrio instável de controlo aproximado
𝑃𝑃𝑒0 𝑚𝑎𝑥 - potência eléctrica máxima no período pré-defeito
𝑃𝑃𝑒𝐷 𝑚𝑎𝑥 - potência eléctrica máxima no período pré-defeito
𝑃𝑃𝑒𝑃 𝑚𝑎𝑥 - potência eléctrica máxima no período pós-defeito
𝑃𝑃𝑒0 - potência eléctrica no período pré-defeito
𝑃𝑃𝑒𝐷 - potência eléctrica no período de defeito
𝑃𝑃𝑒𝑃 - potência eléctrica no período pós-defeito
𝛿𝛿𝑚𝑎𝑥 - ângulo da máquina máximo atingido durante a oscilação
𝜂𝜂 - margem de estabilidade transitória
𝐴𝐴𝑎 - área de aceleração
𝐴𝐴𝑑 - área de desaceleração
∆𝑡𝑡𝑐𝑟 - diferença entre duas estimativas consecutivas de 𝑡𝑡𝑐𝑟
𝛿𝛿𝑒𝑞 - ângulo rotórico da máquina equivalente
𝜔𝑒𝑞 - velocidade angular máquina equivalente
IDCS - índice de detecção instabilidade
IDE - índice de detecção de estabilidade
IDTO - índice auxiliar para determinação do instante óptimo para a identificação do conjunto de
máquinas críticas
𝐼𝐴𝐴𝐶𝑂𝐼𝑖 - índice utilizado na determinação do conjunto de máquinas críticas, nas situações estáveis
𝛿𝛿𝐶 - ângulo rotórico equivalente do conjunto C
𝑀𝑘 - coeficiente de inércia da máquina 𝑘, pertencente ao conjunto C
𝛿𝛿𝑘 - ângulo rotórico da máquina 𝑘, pertencente ao conjunto C
𝑀𝐶 - coeficiente de inércia equivalente do conjunto C
𝛿𝛿𝑅 - ângulo rotórico equivalente do conjunto R
𝑀𝑗 - coeficiente de inércia da máquina 𝑗𝑗, pertencente ao conjunto R
𝛿𝛿𝑗 - ângulo rotórico da máquina 𝑗𝑗, pertencente ao conjunto R
𝑀𝑅 - coeficiente de inércia equivalente do conjunto R
𝜔𝐶 - velocidade angular equivalente do conjunto C
𝜔𝑘 - velocidade angular da máquina 𝑘, pertencente ao conjunto C
𝜔𝑅 - velocidade angular equivalente do conjunto R
𝜔𝑗 - velocidade angular da máquina 𝑗𝑗, pertencente ao conjunto R
𝑀𝑒𝑞 - coeficiente de inércia do sistema
𝑃𝑃𝑚 𝑒𝑞 - potência mecânica da máquina equivalente
𝑃𝑃𝑒 𝑒𝑞 - potência eléctrica da máquina equivalente
XVIII
𝑃𝑃𝑚𝑘 - potência mecânica da máquina 𝑘, pertencente ao conjunto C
𝑃𝑃𝑚𝑗 - potência mecânica da máquina 𝑗𝑗, pertencente ao conjunto R
𝑃𝑃𝑒𝑘 - potência eléctrica da máquina 𝑘, pertencente ao conjunto C
𝑃𝑃𝑒𝑗 - potência eléctrica da máquina 𝑗𝑗, pertencente ao conjunto R
𝛿𝛿𝑒𝑞0 - ângulo rotórico da máquina equivalente no instante pré-defeito
𝛿𝛿𝑒𝑞𝑢 - ângulo rotórico da máquina equivalente, de valor superior a 𝛿𝛿𝑒𝑞0 , para o qual a potência
eléctrica volta a igualar a potência mecânica
𝑃𝑃𝑒 𝑒𝑞𝑚𝑎𝑥 - potência eléctrica máxima fornecida pela máquina equivalente
𝑐𝑐1, 𝑐𝑐2 e 𝑐𝑐3 - constantes da aproximação polinomial
𝛼𝛼 - factor multiplicativo da tolerância
𝑚𝑚𝜂 - declive da recta que relaciona o tempo crítico e a margem de estabilidade transitória
𝑏𝜂 - ordenada na origem da recta que relaciona o tempo crítico e a margem de estabilidade
transitória
∆𝑡𝑡𝑐𝑐𝑐𝑐𝑚𝑚𝑒𝑒𝑡𝑡 - diferença entre os valores do tempo crítico calculados pelo método híbrido e por métodos
de integração numérica
%∆𝑡𝑡𝑐𝑟𝑚𝑒𝑡 - diferença, em percentagem, entre os valores do tempo crítico calculados pelo método
híbrido e por métodos de integração numérica
𝜂𝜂𝑖 - margem de estabilidade transitória da iteração 𝑑𝑑
𝑡𝑡𝑐𝑟𝑖 - estimativa do tempo crítico da iteração 𝑑𝑑
𝑡𝑡𝑐𝑟 ℎ𝑖𝑏 - tempo crítico obtido pelo método híbrido
𝑡𝑡𝑐𝑟 𝑖𝑛 - tempo crítico obtido por métodos de integração numérica
SEE - Sistema Eléctrico de Energia ou Sistema de Energia Eléctrica
COI - Centre Of Inertia
TEF/FET - Função de Energia Transitória
SLEP - Superfície Limite de Energia Potencial
BCU - Boundary Controlling Unstable equilibrium point
SIME - SIngle Machine Equivalent
CMR - Critical Machines Ranking
1
1. Introdução
Este capítulo introdutório inicia-se com uma contextualização da temática abordada
nesta dissertação, apresentando-se uma pequena síntese da evolução que sofreram os
Sistemas de Energia Eléctrica (SEE) e procurando mostrar-se a necessidade de existência de
métodos para o estudo da estabilidade transitória. Apresentam-se também alguns conceitos
referentes à estabilidade dos SEE. De seguida, o ênfase será posto na dissertação em si,
apresentando os seus objectivos e uma pequena introdução a cada um dos capítulos que a
constitui.
1.1. Contextualização
Desde o primeiro SEE a que pode ser dado esse nome, criado por Thomas Edison em
Setembro de 1882 em Nova Iorque, o qual se resumia a um circuito de lâmpadas alimentadas
por uma central eléctrica a funcionar em corrente contínua, os sistemas eléctricos passaram por
profundas mudanças até chegarem aos modernos SEE da actualidade (1). Apesar de no início a
opção recair sobre a corrente contínua, a invenção do transformador aliada à necessidade de
reduzir as perdas em transmissões a grandes distâncias, mais a facilidade de construir motores e
geradores de corrente alternada, fizeram a opção recair definitivamente sobre a corrente
alternada. Seguidamente, devido às economias de escala, o caminho a seguir foi o de construir
centros produtores de elevada capacidade, mas em número reduzido, ligados por uma rede de
alta tensão. Posteriormente, devido à liberalização do mercado e ao crescimento acentuado das
fontes renováveis, assistiu-se a uma descentralização e multiplicação dos centros produtores.
Este crescimento acentuado aliado à imprevisibilidade das fontes, provocou um aumento de
dificuldade na coordenação da exploração do sistema, o qual, associado ao facto de por motivos
económicos se reduzir a reserva girante, o que faz com o que sistema opere próximo dos seus
limites, mostrou a necessidade de estarem disponíveis eficientes métodos para analisar o
funcionamento do sistema na existência de perturbações, pois existe uma redução na garantia
de fiabilidade devido aos factos mencionados.
Como referido, a temática da estabilidade dos SEE veio a revelar-se da maior
importância. Aquando da existência de uma perturbação são dois os principais tipos de
problemas de estabilidade que podem ocorrer: estabilidade de tensão e estabilidade da marcha
síncrona. Enquanto que a primeira está relacionada com desequilíbrios na potência reactiva que
provocam variações nos níveis de tensão, a segunda, e objecto de estudo nesta dissertação,
está relacionada com desequilíbrios na potência activa que provocam variações de frequência,
sendo também designada por estabilidade transitória (2).
Os estudos de estabilidade podem também ser classificados quanto à sua duração,
distinguindo-se: curta duração (até 20 s), média duração (até 5 minutos) e longa duração (até 20
minutos). Para além da duração, e no que diz respeito à estabilidade da marcha síncrona, podem
2
distinguir-se duas situações: as pequenas perturbações, para as quais o sistema pode ser
linearizado, e as grandes perturbações para as quais não é possível linearizar o sistema e que
são estudadas com recurso a diferentes métodos, os quais podem ser divididos em quatro
grupos: Métodos de Integração Numérica, Métodos Directos, Métodos Híbridos e Técnicas de
Inteligência Artificial. Enquanto que os três primeiros são métodos determinísticos, o último,
proposto mais recentemente, é um método probabilístico (3).
1.2. Objectivos e estrutura da dissertação
O trabalho desenvolvido tinha dois objectivos: o estudo dos métodos existentes para a
análise da estabilidade transitória do SEE; a implementação de um método híbrido para ser
usado no estudo da estabilidade transitória do SEE.
O primeiro objectivo consistiu na procura, em bibliografia da especialidade, de métodos
já desenvolvidos e testados.
Para o segundo objectivo procurou-se implementar e testar um método que conjuga as
vantagens do Método das Áreas Iguais com as vantagens dos Métodos de Integração Numérica.
Esta dissertação encontra-se dividida em seis capítulos, os quais são precedidos pelos
Agradecimentos e por um Resumo e um Abstract da dissertação, aos quais se segue o Índice e
as Listas das Figuras, das Tabelas e dos Símbolos e Abreviações presentes no texto. No final
apresentam-se a Bibliografia e os Anexos. Relativamente aos seis capítulos apresenta-se agora
um pequeno resumo de cada um.
No Capítulo 1, a Introdução, expõe-se um breve resumo da história dos SEE e a
necessidade do estudo da estabilidade transitória.
No Capítulo 2, designado Modelização do SEE, apresentam-se os modelos utilizados
para descrever os diferentes componentes da rede e também as contingências a simular.
No Capítulo 3, denominado Métodos de análise da estabilidade transitória do SEE,
apresentam-se os métodos disponíveis actualmente para o estudo da estabilidade transitória.
Iniciando-se pelos Métodos de Integração Numérica, passando pelos Métodos Directos e
Métodos Híbridos e terminando nas Técnicas de Inteligência Artificial.
No Capítulo 4, cujo título é Método híbrido implementado, explica-se em pormenor o
método implementado, dando especial atenção aos conceitos considerados mais importantes.
No Capítulo 5, intitulado Resultados computacionais, são apresentados os resultados
obtidos nos testes computacionais a que o método implementado foi sujeito. Apresenta-se
também um exemplo de aplicação.
3
No último capítulo, o Capítulo 6, designado Conclusões e Propostas para trabalhos
futuros, são apresentadas, para além das conclusões do trabalho desenvolvido, possíveis linhas
de desenvolvimento para trabalhos futuros.
5
2. Modelização do Sistema de Energia Eléctrica
A simulação computacional do comportamento de um SEE exige a modelização de cada
um dos seus componentes. Enquanto que a dinâmica rotacional dos geradores é descrita por
equações diferenciais, os restantes elementos são descritos recorrendo a equações algébricas
(2). Isto acontece pois, na altura em que se iniciaram as simulações dos SEE, a capacidade de
computação era reduzida, o que obrigava ao uso de processos mais leves computacionalmente.
Seguidamente, serão apresentados os modelos que descrevem cada um dos componentes do
SEE (máquinas síncronas, linhas de transmissão, transformadores, cargas), o modelo que
descreve a totalidade da rede e ainda a estratégia para simular as perturbações.
No programa desenvolvido foi utilizado o modelo clássico, o qual, apesar de algumas
limitações, permite uma avaliação válida da estabilidade transitória para a primeira oscilação, o
que compreende um intervalo de tempo não excedendo 2 s (3). À medida que os modelos para
cada um dos componentes forem apresentados, destacar-se-á a diferença introduzida pela
utilização do modelo clássico.
2.1. Máquina Síncrona
Utilizando o modelo clássico, a máquina síncrona é descrita simplesmente por uma força
electromotriz transitória, de módulo constante, em série com a reactância transitória (Figura 2.1)
(2), ou seja,
𝑬𝑬𝒊𝒊′ = 𝑽𝑽𝒊𝒊 + 𝑗𝑗𝑋𝑋𝑑𝑖′ 𝑰𝑰𝒊𝒊 (2.1)
em que 𝑬𝑬𝒊𝒊′ - força electromotriz transitória da máquina 𝑑𝑑 𝑽𝑽𝒊𝒊 - tensão aos terminais da máquina 𝑑𝑑 𝑋𝑋𝑑𝑖′ - reactância transitória da máquina 𝑑𝑑 𝑰𝑰𝒊𝒊 - intensidade de corrente fornecida pela máquina 𝑑𝑑.
A força electromotriz pode ser representada na forma exponencial, ou seja,
𝑬𝑬𝒊𝒊′ = 𝐸𝑖𝑒𝑒𝑗𝛿𝑖 (2.2)
onde 𝐸𝑖 - módulo da força electromotriz transitória da máquina 𝑑𝑑 𝛿𝛿𝑖 - ângulo da força electromotriz transitória da máquina 𝑑𝑑.
𝑰𝑰𝒊𝒊
𝑬𝑬𝒊𝒊′ 𝑽𝑽𝒊𝒊
𝑗𝑗𝑋𝑋𝑑𝑑𝑑𝑑′
Figura 2.1 - Esquema equivalente da máquina síncrona
6
Esta aproximação implica o desprezo quer do sistema de controlo de frequência quer do
de controlo de tensão e, ainda que o desprezo do sistema de controlo de frequência seja válido
para um período até 5 s, o desprezo do sistema de controlo da tensão apenas é válido em
intervalos inferiores a 2 s. Além de desprezar o efeito dos reguladores mencionados, o modelo
usado também não tem em linha de conta aspectos como a saliência dos pólos e a saturação
dos circuitos magnéticos.
Considera-se que o ângulo rotórico máquina 𝑑𝑑 toma o mesmo valor que o ângulo da força
electromotriz transitória da máquina 𝑑𝑑, ou seja, 𝛿𝛿𝑖 refere-se tanto ao ângulo rotórico da máquina 𝑑𝑑
como ao ângulo da força electromotriz transitória da máquina 𝑑𝑑.
O comportamento mecânico da máquina é descrito utilizando a equação de oscilação
2𝐻𝑖𝜔0
𝑑𝑑2𝛿𝛿𝑖𝑑𝑑𝑡𝑡2
+ 𝐷𝐷𝑖𝑑𝑑𝛿𝛿𝑖𝑑𝑑𝑡𝑡
= 𝑃𝑃𝑚𝑖 − 𝑃𝑃𝑒𝑖 (2.3)
onde 𝐻𝑖 - constante de inércia da máquina 𝑑𝑑 𝜔0 - velocidade angular nominal do sistema 𝛿𝛿𝑖 - ângulo rotórico da máquina 𝑑𝑑 𝐷𝐷𝑖 - coeficiente de amortecimento da máquina 𝑑𝑑 𝑃𝑃𝑚𝑖 - potência mecânica da máquina 𝑑𝑑 𝑃𝑃𝑒𝑖 - potência eléctrica gerada pela máquina 𝑑𝑑.
Pode-se também definir a potência de aceleração do gerador 𝑑𝑑
𝑃𝑃𝑎𝑐𝑖 = 𝑃𝑃𝑚𝑖 − 𝑃𝑃𝑒𝑖 (2.4)
Para proceder à integração numérica no domínio do tempo da equação (2.3) é útil
transformá-la em duas equações diferenciais de primeiro grau. Tendo em conta que
𝑑𝑑𝛿𝛿𝑖𝑑𝑑𝑡𝑡
= 𝜔𝑖 − 𝜔0 (2.5)
onde 𝜔𝑖 representa a velocidade angular da máquina 𝑑𝑑, a equação (2.3) transforma-se em
2𝐻𝑖𝜔0
𝑑𝑑𝜔𝑖
𝑑𝑑𝑡𝑡+ 𝐷𝐷𝑖(𝜔𝑖 − 𝜔0) = 𝑃𝑃𝑚𝑖 − 𝑃𝑃𝑒𝑖 (2.6)
O valor do coeficiente de amortecimento e as variações dos valores das constantes de
inércia são usualmente desprezados. O desprezo do coeficiente de amortecimento só é válido se
a análise da estabilidade transitória se restringir à primeira oscilação devendo, caso contrário, ser
considerado.
Habitualmente opta-se por representar as variáveis de estado que descrevem a dinâmica
rotacional dos geradores utilizando como referência o centro de inércia do sistema (COI - Centre
Of Inertia) (3). No desenvolvimento do método híbrido, esta representação é várias vezes
utilizada, sendo por isso importante referir os cálculos necessários para a utilizar.
7
O ângulo do centro de inércia do sistema é, em cada instante, dado por
𝛿𝛿0(𝑡𝑡) =
1𝑀𝑇
�𝑀𝑖𝛿𝛿𝑖
𝑚
𝑖=1
(2.7)
com
𝑀𝑇 = �𝑀𝑖
𝑚
𝑖=1
(2.8)
e 𝑀𝑖 =
2𝐻𝑖𝜔0
=𝐻𝑖𝜋𝑓0
(2.9)
onde 𝑓0 - frequência nominal do sistema 𝑀𝑇 - coeficiente de inércia total 𝑀𝑖 - coeficiente de inércia da máquina 𝑑𝑑.
Desprezando o amortecimento, as equações de oscilação passam então a
�
𝑑𝑑𝜃𝑖𝑑𝑑𝑡𝑡
= 𝜔�𝑖
𝑀𝑖𝑑𝑑𝜔�𝑖𝑑𝑑𝑡𝑡
= 𝑓𝑖(𝜃𝑖)� (2.10)
com
𝜃𝑖 = 𝛿𝛿𝑖 − 𝛿𝛿0 (2.11)
𝜔�𝑖 =
𝑑𝑑𝛿𝛿𝑖𝑑𝑑𝑡𝑡
−𝑑𝑑𝛿𝛿0𝑑𝑑𝑡𝑡
(2.12)
𝑓𝑖(𝜃𝑖) = 𝑃𝑃𝑚𝑖 − 𝑃𝑃𝑒𝑖 −𝑀𝑖
𝑀𝑇𝑃𝑃𝐶𝑂𝐼 (2.13)
𝑃𝑃𝐶𝑂𝐼 = �(𝑃𝑃𝑚𝑖 − 𝑃𝑃𝑒𝑖)
𝑚
𝑖=1
(2.14)
onde 𝜃𝑖 - ângulo rotórico da máquina 𝑑𝑑, referido ao centro de inércia do sistema 𝜔�𝑖 - velocidade angular da máquina 𝑑𝑑, referida ao centro de inércia do sistema 𝑓𝑖(𝜃𝑖) - potência de aceleração da máquina 𝑑𝑑, referida ao centro de inércia do sistema.
2.2. Linha de Transmissão
Em virtude de a frequência se manter aproximadamente constante, a linha eléctrica pode
ser modelada usando parâmetros concentrados (3). Usa-se então o modelo em 𝜋, representado
na Figura 2.2.
8
𝑽𝑽𝒆𝒆 𝑽𝑽𝒓𝒓
𝑰𝑰𝒆𝒆 𝑰𝑰𝒓𝒓
𝑌𝑌𝑇𝑇2� 𝑌𝑌𝑇𝑇
2�
𝑅𝑅𝐿𝐿 𝑗𝑗𝑋𝑋𝐿𝐿
Figura 2.2 - Esquema equivalente em π de uma linha
onde 𝑽𝑽𝒆𝒆 - tensão na emissão 𝑽𝑽𝒓𝒓 - tensão na recepção 𝑰𝑰𝒆𝒆 - intensidade de corrente na emissão 𝑰𝑰𝒓𝒓 - intensidade de corrente na recepção 𝑅𝑅𝐿 - resistência longitudinal 𝑋𝑋𝐿 - reactância longitudinal 𝑌𝑌𝑇 - admitância transversal.
2.3. Transformador
Os transformadores considerados são: transformadores com dois enrolamentos por fase,
transformadores com regulação de tensão e transformadores desfasadores (2).
Desprezando a corrente de magnetização, o transformador com dois enrolamentos por
fase pode ser representado pelo esquema da Figura 2.3.
𝑽𝑽𝒑𝒑 𝑽𝑽𝒔𝒔
𝑰𝑰𝒑𝒑 𝑰𝑰𝒔𝒔 𝑅𝑅𝑇𝑇 𝑗𝑗𝑋𝑋𝑇𝑇
Figura 2.3 - Esquema equivalente simplificado de um transformador
onde 𝑽𝑽𝒑𝒑 - tensão no primário 𝑽𝑽𝒔𝒔 - tensão no secundário 𝑰𝑰𝒑𝒑 - corrente no primário 𝑰𝑰𝒔𝒔 - corrente no secundário 𝑅𝑅𝑇 - resistência total 𝑋𝑋𝑇 - reactância total.
Do esquema retira-se que
𝑽𝑽𝒑𝒑 = 𝑽𝑽𝒔𝒔 + 𝑗𝑗𝒁𝒕𝑰𝑰 (2.15)
tendo em conta que 𝑰𝑰𝒑𝒑 = 𝑰𝑰𝒔𝒔 = 𝑰𝑰 e onde 𝒁𝒕 = 𝑅𝑅𝑡 + 𝑗𝑗𝑋𝑋𝑡 representa a impedância do transformador.
O transformador com regulação de tensão difere do apresentado anteriormente no facto
de ser possível variar a relação de transformação. Esta variação é efectuada usando um
9
comutador de tomadas instalado num dos enrolamentos. Este transformador pode ser modelado
considerando um transformador ideal com relação de transformação 𝑚𝑚′ em série com a
impedância do transformador, também designada impedância de curto-circuito pois é usualmente
medida realizando um ensaio em curto circuito. O esquema do transformador com regulação de
tensão apresenta-se na Figura 2.4.
𝑽𝑽𝒑𝒑 𝑽𝑽𝒔𝒔
𝑰𝑰𝒑𝒑 𝑰𝑰𝒔𝒔 𝑍𝑍𝐶𝐶𝐶𝐶
𝑽𝑽𝒔𝒔′
𝑰𝑰𝒔𝒔′ 𝑚𝑚′ : 1
Figura 2.4 - Esquema equivalente de um transformador com regulação de tensão
onde 𝑽𝑽𝒑𝒑 - tensão no primário 𝑽𝑽𝒔𝒔′ - tensão no secundário do transformador ideal 𝑽𝑽𝒔𝒔 - tensão no secundário 𝑰𝑰𝒑𝒑 - corrente no primário 𝑰𝑰𝒔𝒔′ - corrente no secundário do transformador ideal 𝑰𝑰𝒔𝒔 - corrente no secundário 𝒁𝑪𝑪 - impedância de curto-circuito 𝑚𝑚′ - relação de transformação.
Do esquema retira-se que
𝑽𝑽𝒑𝒑 = 𝑚𝑚′(𝑽𝑽𝒔𝒔 + 𝑗𝑗𝒁𝒄𝒄𝒄𝒄𝑰𝑰𝒔𝒔) (2.16)
A matriz de admitâncias 𝒀𝒀 deste transformador assume a forma de
𝒀𝒀 = �
𝒀𝒀𝒄𝒄𝒄𝒄𝑚𝑚′2 −
𝒀𝒀𝒄𝒄𝒄𝒄𝑚𝑚′
−𝒀𝒀𝒄𝒄𝒄𝒄𝑚𝑚′ 𝒀𝒀𝒄𝒄𝒄𝒄
� (2.17)
sendo 𝒀𝒀𝒄𝒄𝒄𝒄 = 1𝒁𝒄𝒄𝒄𝒄� a admitância de curto-circuito do transformador.
O transformador desfasador apresenta, em relação ao anterior, o facto de a relação de
transformação poder agora ser complexa, ou seja
𝒎𝒎′ = 𝑚𝑚′𝑒𝑒𝑗𝛽
Tal como com o transformador com regulação de tomadas, o esquema do transformador
desfasador é constituído por um transformador ideal com relação de transformação 𝒎𝒎′, agora
complexa, em série com a impedância de curto-circuito, como se observa na Figura 2.5.
10
𝑽𝑽𝒑𝒑 𝑽𝑽𝒔𝒔
𝑰𝑰𝒑𝒑 𝑰𝑰𝒔𝒔 𝑍𝑍𝐶𝐶𝐶𝐶
𝑽𝑽𝒔𝒔′
𝑰𝑰𝒔𝒔′ 𝒎𝒎′ : 1
Figura 2.5 - Esquema equivalente de um transformador desfasador
onde 𝑽𝑽𝒑𝒑 - tensão no primário 𝑽𝑽𝒔𝒔′ - tensão no secundário do transformador ideal 𝑽𝑽𝒔𝒔 - tensão no secundário 𝑰𝑰𝒑𝒑 - corrente no primário 𝑰𝑰𝒔𝒔′ - corrente no secundário do transformador ideal 𝑰𝑰𝒔𝒔 - corrente no secundário 𝒁𝑪𝑪 - impedância de curto-circuito 𝒎𝒎′ - relação de transformação complexa.
Do esquema retira-se que
𝑽𝑽𝒑𝒑 = 𝒎𝒎′(𝑽𝑽𝒔𝒔 + 𝑗𝑗𝒁𝒄𝒄𝒄𝒄𝑰𝑰𝒔𝒔) (2.18)
A matriz de admitâncias 𝒀𝒀 deste transformador assume a forma de
𝒀𝒀 = �
𝒀𝒀𝒄𝒄𝒄𝒄𝑚𝑚′2 −
𝒀𝒀𝒄𝒄𝒄𝒄𝒎𝒎′∗
−𝒀𝒀𝒄𝒄𝒄𝒄𝒎𝒎′ 𝒀𝒀𝒄𝒄𝒄𝒄
� (2.19)
2.4. Carga
Utilizando o modelo clássico todas as cargas do sistema são transformadas em
impedâncias constantes (elasticidade 2), não tendo portanto em conta quer as variações de
tensão quer de frequência (3).
A admitância que descreve a carga ligada ao barramento 𝑘 é obtida por
𝒀𝒀𝒄𝒄𝒄𝒄 = 𝐺𝐺𝑐𝑘 + 𝑗𝑗𝐵𝑐𝑘 =𝑃𝑃𝑐𝑘 − 𝑗𝑗𝑄𝑐𝑘
𝑉𝑘2 (2.20)
com 𝒀𝒀𝒄𝒄𝒄𝒄 - admitância da carga ligada ao barramento 𝑘 𝐺𝐺𝑐𝑘 - parte real de 𝒀𝒀𝒄𝒄𝒄𝒄 𝐵𝑐𝑘 - parte imaginária de 𝒀𝒀𝒄𝒄𝒄𝒄 𝑃𝑃𝑐𝑘 - potência activa de carga do barramento 𝑘 𝑄𝑐𝑘 - potência reactiva de carga do barramento 𝑘 𝑉𝑘 - módulo da tensão no barramento 𝑘.
11
2.5. Modelo Global da Rede
Como referido anteriormente, a rede é representada através de equações algébricas,
sendo caracterizada por uma matriz de admitâncias reduzida. Esta matriz reduzida tem como
base a matriz de admitâncias utilizada no trânsito de energia, a qual é aumentada e
posteriormente reduzida (2). Para efectuar o aumento da matriz de admitâncias utilizada no
trânsito de energia é necessário realizar duas alterações:
• adição da reactância transitória de cada gerador entre os barramentos de
geração e o nó interno do gerador (nó fictício);
• adição das admitâncias equivalentes das cargas.
Em resultado do exposto anteriormente e da utilização do modelo clássico para
descrever os componentes do SEE, este último é representado pelo esquema da Figura 2.6.
𝑗𝑗𝑋𝑋𝑑𝑑1′
𝑬𝑬𝟏𝟏′
𝑰𝑰𝟏𝟏
𝑗𝑗𝑋𝑋𝑑𝑑𝑑𝑑′
𝑬𝑬𝒊𝒊′
𝑰𝑰𝒊𝒊
𝑗𝑗𝑋𝑋𝑑𝑑𝑚𝑚′
𝑬𝑬𝒎𝒎′
𝑰𝑰𝒎𝒎
Rede de Transmissão
𝒀𝒀𝒄𝒄𝟏𝟏
𝒀𝒀𝒄𝒄𝒄𝒄
𝒀𝒀𝒄𝒄𝒄𝒄
Rede Aumentada
Figura 2.6 - Representação esquemática do SEE
Na figura anterior identificam-se: 𝑚𝑚 - número de geradores, equivalente ao número de barramentos de geração 𝑛 - número de barramentos 𝑬𝑬𝒊𝒊′ - força electromotriz transitória da máquina 𝑑𝑑 𝑋𝑋𝑑𝑖′ - reactância transitória da máquina 𝑑𝑑 𝑰𝑰𝒊𝒊 - intensidade de corrente fornecida pela máquina 𝑑𝑑 𝒀𝒀𝒄𝒄𝒄𝒄 - admitância da carga ligada ao barramento 𝑘.
Tendo em conta as alterações referidas anteriormente, as equações nodais da rede
aumentada escrevem-se
�[𝑰𝑰][0]� = �
[𝒀𝒀𝒎𝒎𝒎𝒎] [𝒀𝒀𝒎𝒎𝒄𝒄][𝒀𝒀𝒄𝒄𝒎𝒎] [𝒀𝒀𝒄𝒄𝒄𝒄]� �
[𝑬𝑬′][𝑽𝑽] � (2.21)
onde
12
[𝑬𝑬′] - vector das forças electromotrizes transitórias [𝑽𝑽] - vector das tensões nos barramentos [𝑰𝑰] - vector das correntes injectadas nos barramentos de geração.
Nos barramentos de carga as correntes injectadas são nulas pois as cargas são
representadas por admitâncias constantes.
Partindo do sistema de equações (2.21) obtém-se
[𝑰𝑰] = [𝒀𝒀𝒎𝒎𝒎𝒎][𝑬𝑬′] + [𝒀𝒀𝒎𝒎𝒄𝒄][𝑽𝑽]
[0] = [𝒀𝒀𝒄𝒄𝒎𝒎][𝑬𝑬′] + [𝒀𝒀𝒄𝒄𝒄𝒄][𝑽𝑽] (2.22)
Eliminando [𝑽𝑽] tem-se
[𝑰𝑰] = ([𝒀𝒀𝒎𝒎𝒎𝒎] − [𝒀𝒀𝒎𝒎𝒄𝒄][𝒀𝒀𝒄𝒄𝒄𝒄]−1[𝒀𝒀𝒄𝒄𝒎𝒎])[𝑬𝑬′] (2.23)
A rede pode então ser representada por
[𝑰𝑰] = [𝒀𝒀𝒓𝒓𝒆𝒆𝒅][𝑬𝑬′] (2.24)
sendo 𝒀𝒀𝒓𝒓𝒆𝒆𝒅, a matriz de admitâncias reduzida, dada por
𝒀𝒀𝒓𝒓𝒆𝒆𝒅 = [𝒀𝒀𝒎𝒎𝒎𝒎] − [𝒀𝒀𝒎𝒎𝒄𝒄][𝒀𝒀𝒄𝒄𝒄𝒄]−1[𝒀𝒀𝒄𝒄𝒎𝒎] (2.25)
A matriz de admitâncias 𝒀𝒀 que descreve a rede, neste caso 𝒀𝒀𝒓𝒓𝒆𝒆𝒅, pode ser dividida na
sua parte real, matriz de condutâncias, e parte imaginária, matriz de susceptâncias, ou seja
𝒀𝒀 = 𝐺𝐺 + 𝑗𝑗𝐵 (2.26)
A potência activa gerada pelo gerador 𝑑𝑑 pode ser determinada utilizando a expressão
𝑃𝑃𝑒𝑖 = Re{𝑬𝑬𝒊𝒊′𝑰𝑰𝒊𝒊∗} (2.27)
sendo 𝑬𝑬𝒊𝒊′ a respectiva força electromotriz transitória e 𝑰𝑰𝒊𝒊 a corrente injectada no respectivo
barramento, obtida a partir da primeira equação de (2.22). A equação (2.27) pode ser
desenvolvida, obtendo-se
𝑃𝑃𝑒𝑖 = 𝐸𝑖2𝐺𝐺𝑖𝑖 + �𝐸𝑖𝐸𝑗�𝐺𝐺𝑖𝑗 cos�𝛿𝛿𝑖 − 𝛿𝛿𝑗� + 𝐵𝑖𝑗 sin�𝛿𝛿𝑖 − 𝛿𝛿𝑗��
𝑚
𝑗=1𝑗≠𝑖
(2.28)
onde 𝐺𝐺𝑖𝑗 e 𝐵𝑖𝑗 se referem à posição 𝑑𝑑𝑗𝑗 da matriz 𝒀𝒀, respectivamente, de condutâncias e
susceptâncias.
2.6. Contingências Simuladas
As contingências aqui abordadas resumem-se a curto-circuitos trifásicos simétricos na
extremidades das linhas junto aos barramentos.
13
Em termos temporais, considera-se que a perturbação ocorre sempre 0,1 s após o início
da simulação, ou seja, 𝑡𝑡𝑑𝑒𝑓 = 0,1 s. Quando ocorre o defeito, o seu efeito é traduzido
numericamente pela adição, ao nó respectivo na matriz de admitâncias 𝒀𝒀 que descreve a rede,
de uma admitância de defeito 𝒀𝒀𝒅𝒆𝒆𝒇 de módulo elevado, por exemplo �𝒀𝒀𝒅𝒆𝒆𝒇� = 108. O defeito é
eliminado quando as protecções actuam, no instante de tempo 𝑡𝑡𝑐𝑙 (clearing time), podendo
ocorrer uma de duas situações: o defeito é eliminado espontaneamente ou o defeito é eliminado
pela retirada de uma linha de serviço, sendo necessário reflectir esse facto alterando a matriz 𝒀𝒀.
2.7. Conclusões
Neste capítulo foram introduzidos os modelos que permitem descrever o comportamento
dos diferentes componentes da rede (máquinas síncronas, linhas, transformadores, cargas) e da
própria rede. Foi também referida a noção de centro de inércia do sistema e quais as alterações
provocadas nas equações de oscilação. Por fim apresentou-se a estratégia adoptada para
introduzir o efeito provocado pelas contingências simuladas.
Como referido, a utilização do modelo clássico, apenas é válida em intervalos inferiores a
2 s, o que corresponde à primeira oscilação. É portanto necessário, caso se pretendam simular
maiores intervalos de tempo, a modelação de elementos como os sistemas de controlo de
frequência e de tensão, além de uma modelação mais precisa das cargas que tenha em conta as
variações de frequência e de tensão.
15
3. Métodos de análise da estabilidade transitória do SEE
Na fase de planeamento, e também durante a exploração, é necessário garantir a
capacidade do SEE suportar determinadas perturbações graves, como por exemplo, a ocorrência
de curto-circuitos ou a perca de grupos geradores significativos. A resposta do SEE quando
sujeito a estas perturbações pode ser simulada recorrendo a métodos computacionais, que
descrevem o comportamento dinâmico do sistema e permitem assim tirar conclusões acerca da
sua estabilidade transitória.
Seguidamente serão apresentados alguns métodos desenvolvidos com o propósito de
avaliar a estabilidade transitória do sistema. Serão apresentados inicialmente os Métodos de
Integração Numérica e os Métodos Directos. De seguida serão apresentados métodos mais
recentes como os Métodos Híbridos e as Técnicas de Inteligência Artificial. Por fim, apresenta-se
uma comparação entre os diferentes métodos apresentados.
3.1. Métodos de Integração Numérica
Um SEE é descrito por dois sistemas de equações (3). Genericamente, tem-se um
sistema de equações diferenciais
�̇� = 𝑓(𝑚𝑚,𝑦, 𝑡𝑡) (3.1)
e um sistema de equações algébricas não-lineares
𝑔(𝑚𝑚,𝑦, 𝑡𝑡) = 0 (3.2)
O primeiro sistema inclui a dinâmica rotacional dos geradores, bem como as equações
diferenciais que caracterizam os sistemas de regulação de tensão e de velocidade, caso estes
sejam considerados. Quanto maior for o detalhe da modelização maior será a ordem deste
sistema. O segundo sistema descreve o comportamento da rede.
O objectivo desta abordagem é resolver, no domínio do tempo, os sistemas de equações
acima referidos. Os algoritmos utilizados podem ser classificados quanto ao método utilizado na
integração numérica das equações diferenciais e quanto ao processo para resolver o conjunto
dos sistemas de equações. A primeira classificação divide os métodos em métodos implícitos ou
explícitos caso a solução obtida num instante de tempo dependa ou não, respectivamente, das
grandezas nesse mesmo instante de tempo (4). A segunda classificação divide os algoritmos em
dois grupos: alternado caso os sistemas sejam resolvidos separada e alternadamente e
simultâneo caso estes sejam resolvidos em conjunto.
Na Figura 3.1 é apresentado o fluxograma genérico usado para implementar um
programa de simulação usando métodos de integração numérica (2). Após realizar um trânsito
de energia, que fornece os valores pré-defeito das tensões na rede e potências geradas,
16
efectua-se o aumento e posterior redução da matriz de admitâncias nodais, calculando-se então
os valores pré-defeito das forças electromotrizes dos geradores. Entra-se seguidamente num
ciclo durante o qual o tempo é incrementado até se atingir o tempo máximo de simulação
estipulado 𝑡𝑡𝑚𝑎𝑥. Em cada iteração deste ciclo altera-se a matriz de admitâncias caso tenham
ocorrido alterações da topologia da rede, calculando-se de seguida os novos valores das
variáveis.
Trânsito de Energia
Alteração e redução da matriz de admitâncias
Cálculo das f.e.m. dos geradores
t = 0
Alteração da rede?
Alteração e redução da matriz de admitâncias
Calcular os valores das variáveis
t = t + ∆t
t ≤ tmax
Fim
Sim
Sim
Não
Não
Início
Figura 3.1 - Fluxograma para simulação da estabilidade transitória
Um factor fundamental na análise da estabilidade transitória de um SEE é o tempo crítico
de actuação das protecções, 𝑡𝑡𝑐𝑟, ou seja, o período máximo a que o sistema pode estar sujeito
ao defeito mantendo a estabilidade. A determinação do tempo crítico pode ser efectuada
utilizando um algoritmo de bissecção que efectua simulações consecutivas, com diferentes
17
tempos de actuação das protecções, 𝑡𝑡𝑐𝑙, até se atingir a convergência. Na Figura 3.2 é
apresentado um fluxograma para a determinação do tempo crítico 𝑡𝑡𝑐𝑟, utilizando uma estimativa
inicial 𝑡𝑡𝑐𝑙0, que corresponde a uma situação instável (3).
tinst = tcl0test = 0tcl = tcl0
Integração Numérica
Estável?
tcr = test
Fim
SimNão
Início
tinst - test < ɛ
test =tcltinst =tcl
Não
Sim
tcl = (tinst + test)/2
Figura 3.2 - Fluxograma do método da bissecção para cálculo do tempo crítico
Como se pode observar pelo fluxograma, os valores do tempo de actuação das
protecções correspondentes às situações instáveis e estáveis, respectivamente 𝑡𝑡𝑖𝑛𝑠𝑡 e 𝑡𝑡𝑒𝑠𝑡, são
actualizados à medida que se realizam novas iterações do ciclo. Para cada iteração, para além
da primeira, o tempo de actuação das protecções é a média entre 𝑡𝑡𝑖𝑛𝑠𝑡 e 𝑡𝑡𝑒𝑠𝑡. À medida que o
ciclo corre observa-se que tanto 𝑡𝑡𝑖𝑛𝑠𝑡 como 𝑡𝑡𝑒𝑠𝑡 convergem para o valor do tempo crítico a
determinar, sendo a simulação interrompida quando a diferença entre estes dois tempos for
menor que uma tolerância 𝜀𝜀.
3.2. Métodos Directos
A maioria dos métodos directos para análise da estabilidade transitória dos SEE baseia-
se no segundo teorema de Lyapunov. Apesar de, aquando da sua publicação, os resultados de
Lyapunov não terem recebido grande atenção, são hoje em dia uma área de investigação e com
resultados garantidos na análise da estabilidade transitória dos SEE (5).
Estes métodos caracterizam-se por apenas ser necessário integrar as equações
diferenciais que descrevem o sistema durante o período de permanência no defeito, o que
18
permite uma redução assinalável dos tempos de computação. Apresentam contudo a
desvantagem de apresentarem problemas de modelação do sistema, apenas permitindo uma
modelação pouco pormenorizada, e de os métodos desenvolvidos não serem totalmente fiáveis,
especialmente quando o sistema opera próximo dos seus limites (6).
Com base no teorema de Lyapunov foram desenvolvidos, ao longo do tempo, outros
métodos directos, com o objectivo de contrariar limitações e dificuldades de aplicação dos
métodos existentes. Surge então, a seguir ao Método de Lyapunov, o Método da Função de
Energia Transitória, seguido dos Métodos da Superfície Limite de Energia Potencial e do Ponto
de Equilíbrio Instável de Controlo.
Seguidamente, para além destes métodos serem apresentados, é também apresentado
o Método das Áreas Iguais. Para alguns dos métodos é ainda apresentado um exemplo de
aplicação.
3.2.1. Método de Lyapunov
O Método de Lyapunov (3), apresentado por A. M. Lyapunov na sua dissertação de
doutoramento, estipula que a estabilidade de um sistema físico, de dimensão 𝑚𝑚, descrito por
�̇�𝑚 = 𝑓(𝑚𝑚), 𝑓(0) = 0 (3.3)
pode ser verificada sem integração numérica sendo apenas necessário garantir que a energia
total 𝐸(𝑚𝑚) do sistema se mantenha continuamente decrescente no tempo. É então necessário
que, excepto no(s) ponto(s) de equilíbrio 𝑚𝑚𝑒, a derivada temporal da energia 𝐸(𝑚𝑚) seja negativa.
Esta condição pode ser formalizada matematicamente. Usando a função de Lyapunov
𝑉(𝑚𝑚) para representar 𝐸(𝑚𝑚), para garantir que o sistema é estável é necessário garantir que:
a) a função 𝑉(𝑚𝑚) é definida positiva na vizinhança do ponto de equilíbrio 𝑚𝑚𝑒, ou seja,
𝑉(𝑚𝑚) > 0 excepto para 𝑚𝑚 = 𝑚𝑚𝑒 em que 𝑉(𝑚𝑚𝑒) = 0;
b) a derivada temporal de 𝑉(𝑚𝑚), �̇�(𝑚𝑚), é semi-definida negativa anulando-se para
𝑚𝑚 = 𝑚𝑚𝑒.
Caso �̇�(𝑚𝑚) seja definida negativa o sistema considera-se assimptoticamente estável.
A aplicação do método de Lyapunov, para avaliar a estabilidade de um sistema, consiste
em quatro etapas distintas:
a) cálculo da região de estabilidade em torno do ponto de equilíbrio 𝑚𝑚𝑒 para a
configuração pós-defeito do sistema;
b) formulação da função de Lyapunov 𝑉(𝑚𝑚) para a configuração pós-defeito do
sistema. Esta função é, geralmente, a soma das energias cinética e potencial do
sistema na configuração pós-defeito, 𝑉(𝑚𝑚) = 𝑉𝐾𝐸 + 𝑉𝑃𝐸, pois é a formulação que
apresenta melhores resultados;
19
c) integração das equações que descrevem o sistema durante o período de defeito,
até ao instante de eliminação deste, calculando o valor da função de Lyapunov
nesse ponto;
d) cálculo do valor crítico da função de Lyapunov, 𝑉𝑐𝑟, para a configuração de
defeito do sistema.
Se o valor da função no momento de eliminação do defeito for menor ou igual ao valor
crítico, isto é, 𝑉𝑐𝑙 ≤ 𝑉𝑐𝑟, o sistema é considerado estável, caso isso não se verifique o sistema é
considerado instável. Esta formulação é genérica aos métodos baseados no Método de
Lyapunov, os quais usam funções de energia.
Como referido anteriormente, o Método de Lyapunov apresenta, assim como todas as
formulações directas, a vantagem de não ser necessário realizar a integração numérica das
equações do sistema no período pós-defeito, ou seja, não é necessário conhecer quer a
evolução temporal dos ângulos rotóricos, quer as velocidades angulares dos geradores, para se
poder concluir se o sistema é estável ou não.
Contudo, este método apresenta algumas desvantagens. Para sistemas de grandes
dimensões, como é o caso dos SEE, é difícil a determinação de todos os pontos singulares das
equações que descrevem o modelo sendo, por vezes, necessário recorrer a métodos numéricos.
Além disso, existe ainda a dificuldade para formular a função de Lyapunov a partir das equações
que descrevem o comportamento do SEE, devido à falta de procedimentos formais. Outras das
desvantagens é o facto de que a verificação da estabilidade transitória do sistema pelo Método
de Lyapunov é uma condição suficiente mas não necessária, já que este método dá uma noção
pessimista acerca da estabilidade do sistema ao cobrir apenas uma secção da região de
estabilidade. Assim sendo, um sistema considerado instável pelo Método de Lyapunov pode não
o ser. Refere-se ainda que esta formulação não permite obter uma solução ultra-rápida.
Para contrariar as dificuldades apresentadas por este método, novos métodos foram
desenvolvidos, entre os quais o Método da Função de Energia Transitória.
3.2.2. Método da Função de Energia Transitória (TEF)
Devido às semelhanças, o Método da Função de Energia Transitória tornou-se, na
literatura sobre SEE, um sinónimo do Método de Lyapunov. Este método, que se pode
considerar um caso particular do Método de Lyapunov caso se desprezem as condutâncias,
caracteriza-se por utilizar uma função, a Função de Energia Transitória, que se obtém por
integração das equações do movimento que descrevem o comportamento do sistema.
Para representar o comportamento do sistema é necessário usar dois modelos
matemáticos, um para o período de defeito e outro para o período pós-defeito. No instante
anterior ao defeito considera-se que o sistema se encontra em regime permanente.
Determinando os desvios angulares usando como referência o centro de inércia do sistema,
tem-se, para o período de defeito:
20
𝑑𝑑𝜃𝑖𝑑𝑑𝑡𝑡
= 𝜔�𝑖
𝑀𝑖𝑑𝑑𝜔�𝑖𝑑𝑑𝑡𝑡
= 𝑓𝑖𝐷(𝜃) com 𝑡𝑡 ∈ ]0, 𝑡𝑡𝑐𝑙]
(3.4)
Enquanto que para o período pós-defeito se tem:
𝑑𝑑𝜃𝑖𝑑𝑑𝑡𝑡
= 𝜔�𝑖
𝑀𝑖𝑑𝑑𝜔�𝑖𝑑𝑑𝑡𝑡
= 𝑓𝑖𝑃(𝜃) com 𝑡𝑡 ∈ ]𝑡𝑡𝑐𝑙 , +∞[
(3.5)
Em ambos os casos 𝑑𝑑 = 1, 2, … ,𝑚𝑚.
Continuando a usar como referência o centro de inércia do sistema, o ângulo rotórico da
máquina 𝑑𝑑 no ponto de equilíbrio estável, 𝜃𝑖𝑠, para a configuração pós-defeito, é determinado
resolvendo as equações algébricas não-lineares:
𝑓𝑖(𝜃) = 𝑃𝑃𝑚𝑖 − 𝐸𝑖2𝐺𝐺𝑖𝑖 −��𝐸𝑖𝐸𝑗𝐵𝑖𝑗 sin𝜃𝑖𝑗 + 𝐸𝑖𝐸𝑗𝐺𝐺𝑖𝑗 cos𝜃𝑖𝑗�
𝑚
𝑗=1𝑗≠1
−𝑀𝑖
𝑀𝑇𝑃𝑃𝐶𝑂𝐼 = 0 (3.6)
com 𝑑𝑑 = 1, 2, … ,𝑚𝑚 e 𝜃𝑖𝑗 = 𝜃𝑖 − 𝜃𝑗.
Considera-se que o integral de movimento do sistema de equações (3.5) é uma função
de energia apropriada, cuja derivada é:
𝑑𝑑𝑡𝑡 =
𝑀𝑖𝑑𝑑𝜔�𝑖𝑓𝑖(𝜃) =
𝑑𝑑𝜃𝑖𝜔�𝑖
(3.7)
Integrando as equações de cada máquina entre (𝜃𝑖𝑠, 0), ponto de equilíbrio estável do
sistema na configuração pós defeito, e (𝜃𝑖 ,𝜔�𝑖):
𝑉𝑖(𝜃,𝜔�) =
12𝑀𝑖𝜔�𝑖2 − � 𝑓𝑖(𝜃)𝑑𝑑
𝜃𝑖
𝜃𝑖𝑠
𝜃𝑖 , 𝑑𝑑 = 1, … ,𝑚𝑚 (3.8)
Esta expressão é conhecida como a função de energia de cada máquina do SEE.
Adicionando as 𝑚𝑚 funções obtém-se o integral de movimento do sistema:
𝑉(𝜃,𝜔�) =
12�𝑀𝑖𝜔�𝑖2𝑚
𝑖=1
−�� 𝑓𝑖(𝜃)𝑑𝑑𝜃𝑖
𝜃𝑖𝑠
𝜃𝑖
𝑚
𝑖=1
(3.9)
Comparando a expressão anterior com:
𝑉(𝜃,𝜔�) = 𝑉𝐾𝐸(𝜔�) + 𝑉𝑃𝐸(𝜃) (3.10)
21
conclui-se que o primeiro termo do lado direito da equação (3.9) corresponde à energia cinética
do sistema enquanto que o segundo termo corresponde à energia potencial.
Como referido anteriormente, uma das etapas nos métodos baseados no Método de
Lyapunov, é o cálculo do valor crítico da função de energia 𝑉𝑐𝑟. No Método da Função de Energia
Transitória é necessário, para determinar este valor, calcular o ponto de equilíbrio instável de
controlo 𝜃𝑢, pois
𝑉𝑐𝑟 = 𝑉𝑃𝐸(𝜃𝑢) (3.11)
Acontece que este cálculo, inserido na determinação do domínio de estabilidade, é uma
fase que exige elevado esforço computacional, o que é um inconveniente à aplicação deste
método. Para obviar a este inconveniente, foram desenvolvidos novos métodos, como o da
Superfície Limite de Energia Potencial.
3.2.3. Superfície Limite de Energia Potencial (SLEP)
O Método da Superfície Limite de Energia Potencial (7) apresenta, em relação ao método
anterior, a vantagem de não ser necessário calcular o ponto de equilíbrio instável de controlo 𝜃𝑢,
sendo apenas necessário realizar a integração das equações do modelo que descreve o sistema
no período de defeito para determinar 𝑉𝑐𝑟, facto muito útil em sistemas de grande dimensão. Por
vezes, é ainda possível evitar o cálculo do ponto de equilíbrio estável do sistema na configuração
pós-defeito, 𝜃𝑠.
Num sistema constituído por uma máquina ligada a um barramento de potência de curto-
circuito infinita, também designado barramento infinito, este método pode ser facilmente aplicado
considerando que o valor de energia crítico, correspondente ao cruzamento da SLEP, é igual ao
valor máximo que a energia potencial 𝑉𝑃𝐸 atinge durante a trajectória do período de defeito, ou
seja, 𝑉𝑐𝑟 = 𝑉𝑃𝐸𝑚𝑎𝑥(𝜃). Basta então integrar as equações que descrevem o funcionamento do
sistema até se atingir um máximo da energia potencial 𝑉𝑃𝐸. De seguida, considerando um tempo
de actuação das protecções de 𝑡𝑡𝑐𝑙, calcula-se o valor da função de energia 𝑉(𝜃𝑐𝑙 ,𝜔�𝑐𝑙) e
compara-se esse valor com 𝑉𝑐𝑟, caso 𝑉(𝜃𝑐𝑙 ,𝜔�𝑐𝑙) < 𝑉𝑐𝑟 o sistema é estável, caso contrário é
instável. Pode-se também, caso se pretenda determinar o tempo crítico 𝑡𝑡𝑐𝑟, integrar as equações
do sistema, calculando em cada passo o valor de 𝑉(𝜃𝑐𝑙 ,𝜔�𝑐𝑙) até que 𝑉(𝜃𝑐𝑙 ,𝜔�𝑐𝑙) = 𝑉𝑐𝑟.
Apresenta-se de seguida um exemplo de aplicação deste método.
Exemplo
Considere-se um sistema constituído por uma máquina síncrona, ligada a um barramento
de potência de curto-circuito infinita, como se pode observar na Figura 3.3.
Considera-se que no instante de tempo 𝑡𝑡 = 0 ocorre um curto-circuito trifásico simétrico
na linha, junto ao barramento 2, e a linha é retirada de serviço por abertura dos disjuntores. No
instante de tempo 𝑡𝑡𝑐𝑙 religam-se os disjuntores e considera-se que o defeito foi eliminado.
22
Dados: • Gerador
• Transformador:
• Linha
𝑋𝑋𝑑′ = 0,35 pu 𝑀 = 0,02 s2
𝑋𝑋𝑐𝑐 = 0,15 pu
𝑋𝑋𝑙 = 0,09 pu
Sabendo que, no instante 𝑡𝑡 = 0, a tensão e potência recebida no barramento 3 são,
respectivamente, 𝑽𝑽𝟑 = 1,0 pu e 𝑺𝟑 = 0,8 + 𝑗𝑗0,4 pu, pretende-se saber qual o tempo crítico de
religação dos disjuntores.
GT
L1 2 3
Figura 3.3 - Rede utilizada no exemplo do Método da SLEP
Convém referir desde já que, dado o exemplo considerar um barramento de potência de
curto-circuito infinita com uma constante de inércia infinita, as variáveis correspondentes ao
centro de inércia do sistema têm os mesmos valores que as variáveis do barramento de potência
de curto-circuito infinita, ou seja, 𝛿𝛿0 = 0. Consequentemente 𝜔�1 = 𝑑𝛿1𝑑𝑡
e 𝛿𝛿1 = 𝜃1.
Partindo do estado pré-defeito calculam-se as variáveis consideradas constantes no
tempo, a potência mecânica e o módulo da força electromotriz, e ainda o valor inicial do ângulo
da força electromotriz.
𝑃𝑃𝑚 = 0,8 pu
𝑰𝑰 = �𝑺𝟑𝑽𝑽𝟑�∗
= 0,8 − 𝑗𝑗0,4 pu
𝑬𝑬𝟏𝟏 = 𝑽𝑽𝟑 + 𝑗𝑗(𝑋𝑋𝑑′ + 𝑋𝑋𝑐𝑐 + 𝑋𝑋𝑙)𝑰𝑰 = 1,323𝑒𝑒𝑗0,364 ⇒ �𝐸1 = 1,323 pu𝛿𝛿0 = 0,364 rad
�
A partir da equação (3.6) calcula-se o ponto de equilíbrio estável da situação pós-defeito.
Assumindo que para o período pós-defeito a potência fornecida pela máquina é dada por
𝑃𝑃𝑒𝑃 = 𝑃𝑃𝑒𝑃 𝑚𝑎𝑥 sin 𝛿𝛿1 =𝐸1𝑉3
𝑋𝑋𝑑′ + 𝑋𝑋𝑐𝑐 + 𝑋𝑋𝑙sin𝛿𝛿1
a equação (3.6) assume a forma
𝑃𝑃𝑚 −𝐸1𝑉3
𝑋𝑋𝑑′ + 𝑋𝑋𝑐𝑐 + 𝑋𝑋𝑙sin𝛿𝛿1𝑠 = 0
obtendo-se
𝛿𝛿1𝑠 = 0,364 rad
Agora é necessário encontrar a expressão para o cálculo da energia potencial. Usando o
segundo termo da equação (3.8) tem-se
𝑉𝑃𝐸(𝜃1) = −� 𝑓1(𝜃1)𝜃1
𝜃1𝑠
𝑑𝑑𝜃1 = −� (𝑃𝑃𝑚 − 𝑃𝑃𝑒𝑃 𝑚𝑎𝑥 sin𝜃1)𝜃1
𝜃1𝑠
𝑑𝑑𝜃1
23
e obtém-se
������� = −���� − ���� −
�����cos �� − cos ����
Calculando ������� para diferentes valores de �� em torno do ponto de equilíbrio estável
obtém-se o gráfico da Figura 3.4, cujo máximo é igual a ������ = 2,260 e, portanto, o valor de
energia crítico é ��� = ������ = 2,260.
Figura 3.4 - Gráfico de �������
O último passo para o cálculo do tempo crítico ��� é integrar as equações que descrevem
o funcionamento do sistema, durante o período de defeito, até que o valor da função de energia
����, !�� seja igual a ���. A função de energia é obtida através da equação (3.8) assumindo a
forma de
����, !�� =1
2#� !�
$ − ���� − ���� −
�����cos �� − cos ����
onde se usou o resultado obtido anteriormente neste exemplo para o cálculo da energia
potencial. As equações que descrevem o funcionamento do sistema no período de defeito,
considerando que a potência eléctrica fornecida é nula, são
!���� =�
#�
%��� =1
2
��$
#+ %�0�
sendo obtidas por integração da equação de oscilação, que no período de defeito se resume a
#'$%
'�$= �
Calculando o valor da função de energia, ao longo do processo de integração, obtém-se
o gráfico da Figura 3.5, de onde se conclui que ����, !�� = ��� para � ≅ 253ms.
24
���
Figura 3.5 - Gráfico de �����
O valor do tempo crítico é então � ≅ 253ms. ∎
Num sistema multimáquina, cada perturbação dá origem a um modo de instabilidade,
pois, consoante a localização e natureza da perturbação, a perda de sincronismo pode ser
originada por uma ou mais máquinas se tornarem instáveis. Associado a cada modo de
instabilidade está um ponto de equilíbrio instável, denominado ponto de equilíbrio instável de
controlo para essa perturbação.
Em torno do ponto de equilíbrio estável para a configuração pós-defeito do sistema há
vários pontos instáveis de equilíbrio, os quais são a solução das equações (3.6). Após a
eliminação de uma perturbação, caso o sistema seja instável, a sua trajectória, partindo do ponto
de equilíbrio estável pré-defeito, vai aproximar-se de um determinado ponto de equilíbrio instável,
sendo esse o ponto de equilíbrio instável de controlo para essa perturbação.
Na Figura 3.6 encontra-se representada, para um sistema constituído por três máquinas,
a evolução da energia potencial ������ em função dos ângulos ��e �� de duas das máquinas
referidos ao centro de inércia. Observam-se, além das superfícies equipotenciais, três pontos de
equilíbrio instável, ��, �� e ��, em torno do ponto de equilíbrio estável, ��, cada um associado a
uma determinada perturbação, e os quais podem ser pontos de sela (�� e ��) ou máximos
relativos (��). A linha a tracejado que passa pelos três pontos de equilíbrio instável e é
perpendicular às superfícies equipotenciais denomina-se Superfície Limite de Energia Potencial.
Para uma dada perturbação, se no momento em que esta for eliminada a trajectória do
sistema tiver ultrapassado a SLEP, o sistema é instável. Se a perturbação for eliminada
suficientemente cedo então a trajectória oscilará em torno do ponto de equilíbrio estável, ��, sem
ultrapassar a SLEP.
Para calcular o valor de �, a trajectória do sistema no período de defeito é monitorizada
até interceptar a SLEP no ponto �∗. Tal como no Método da Função de Energia Transitória,
25
também no Método da SLEP, o cálculo do valor crítico da função de energia, 𝑉𝑐𝑟, depende do
ponto de equilíbrio instável de controlo, 𝜃𝑢, 𝑉𝑐𝑟 = 𝑉𝑃𝐸(𝜃𝑢). Acontece que, em muitos casos, o
ponto de equilíbrio instável de controlo, 𝜃𝑢, está próximo de 𝜃∗, e então 𝑉𝑐𝑟 = 𝑉𝑃𝐸(𝜃𝑢) ≈ 𝑉𝑃𝐸(𝜃∗),
sendo esta a base do Método da SLEP e daí a determinação do cruzamento da SLEP ser uma
questão chave.
Figura 3.6 - Superfície Limite de Energia Potencial
Os passos para estabelecer os domínios de estabilidade transitória do sistema para uma
determinada perturbação, usando o Método da SLEP, são os seguintes:
1. determinar o ponto de equilíbrio estável, 𝜃𝑠, a partir da equação (3.6);
2. calcular a trajectória durante o período de defeito, através de integração
numérica das equações (3.4);
3. monitorizar quando se dá o cruzamento da SLEP, altura em que os valores de
𝜃(𝑡𝑡) definem 𝜃∗ a partir do qual se calcula 𝑉𝑃𝐸(𝜃∗) que constitui uma boa
aproximação de 𝑉𝑐𝑟; o cruzamento da SLEP é verificado pela mudança de sinal
de 𝑓𝑇(𝜃)(𝜃 − 𝜃𝑠).
Para determinar o tempo crítico, 𝑡𝑡𝑐𝑟, integram-se as equações que descrevem o sistema
no período de defeito até que 𝑉(𝜃,𝜔�) = 𝑉𝑐𝑟. Para determinar se, para uma determinada situação,
o sistema é estável, seguem-se os seguintes passos:
1. calcula-se a trajectória do sistema durante o período de defeito até ao momento
de eliminação do defeito, 𝑡𝑡𝑐𝑙, obtendo-se os valores de 𝜃(𝑡𝑡𝑐𝑙) e 𝜔�(𝑡𝑡𝑐𝑙);
2. determina-se o valor da função de energia para esse ponto, 𝑉𝑐𝑙 = 𝑉𝐾𝐸(𝜔�𝑐𝑙) +
𝑉𝑃𝐸(𝜃𝑐𝑙);
3. verifica-se a estabilidade:
• se 𝑉𝑐𝑙 ≤ 𝑉𝑐𝑟, o sistema é estável;
26
• se 𝑉𝑐𝑙 > 𝑉𝑐𝑟, o sistema é instável.
Como referido anteriormente este método evita por vezes o cálculo do ponto de equilíbrio
estável, o que acontece quando 𝜃𝑠 está próximo de 𝜃0.
3.2.4. Método do Ponto de Equilíbrio Instável de Controlo (BCU)
Este método, conhecido na literatura da especialidade como Boundary Controlling
Unstable equilibrium point (BCU), consiste numa outra abordagem de aplicação de funções de
energia à análise de estabilidade de um SEE (7). A principal característica deste método consiste
no uso da noção de gradiente durante o período pós-defeito, o que permite reduzir a ordem do
sistema, fornecendo um novo algoritmo para determinar o ponto de equilíbrio instável de
controlo.
A aplicação deste método consiste nos seguintes passos:
1. determina-se o ponto de equilíbrio estável, 𝜃𝑠, a partir da equação (3.6);
2. calcula-se o ponto de equilíbrio instável de controlo da seguinte forma:
a. integra-se o sistema de equações (3.4) que descreve o sistema no
período de defeito, calculando em cada passo de integração o valor de
𝑉(𝜃,𝜔�), determinando-se, tal como no Método da SLEP, quando se dá o
cruzamento da SLEP, obtendo-se 𝜃∗.
b. após o cruzamento da SLEP, termina-se a integração do sistema de
equações que descreve o sistema no período de defeito, integrando-se
as equações dinâmicas de gradientes do sistema correspondentes ao
período pós-defeito:
�̇� = 𝑓(𝜃), 𝜃(𝑡𝑡∗) = 𝜃∗ (3.12)
as quais definem um sistema de ordem reduzida no qual apenas se
considera a evolução dinâmica de 𝜃;
c. integra-se (3.12) enquanto se monitoriza o valor de
‖𝑓(𝜃)‖ = �|𝑓𝑖(𝜃)|
𝑚
𝑖=1
(3.13)
até se atingir o primeiro mínimo, altura em que se interrompe o processo
e se obtém 𝜃 = 𝜃𝑎𝑝𝑝𝑢 . Este ponto encontra-se na vizinhança do ponto de
equilíbrio instável de controlo pelo que 𝑉𝑐𝑟 = 𝑉𝑃𝐸�𝜃𝑎𝑝𝑝𝑢 � é uma boa
aproximação do valor da energia crítica do sistema;
d. o ponto de equilíbrio de controlo exacto obtém-se a partir das equações
(3.6) tendo como ponto de partida 𝜃𝑎𝑝𝑝𝑢 , obtendo-se 𝜃𝑢;
3. calcula-se o valor da energia crítica do sistema, através de 𝑉𝑐𝑟 = 𝑉𝑃𝐸(𝜃𝑢).
27
4. determina-se o valor de ��� utilizando os valores de ���, � previamente
calculados em 2. e verifica-se o instante de tempo em que ���, � = ���,
obtendo-se ���. Se ��� < ��� o sistema é estável.
3.2.5. Método das Áreas Iguais
Este método que, como os apresentados anteriormente, é baseado em considerações
energéticas, pode ser aplicado a um sistema constituído por uma máquina síncrona ligada a um
barramento infinito (2).
Considera-se que no período pré-defeito o sistema se encontra em regime permanente,
ou seja, a potência eléctrica é igual à potência mecânica:
�� = ��� (3.14)
Sendo a potência eléctrica, no período pré-defeito, dada por:
��� = ��
���� sin � (3.15)
Enquanto que durante o período de defeito o sistema é descrito por:
����
���= �� − ��
���� sin � (3.16)
Por fim, no período pós-defeito é descrito por:
����
���= �� − ��
���� sin � (3.17)
sendo que ������, ��
���� e ������ são constantes que dependem da configuração da rede.
Esta situação está representada na Figura 3.7.
�
�
�
��
�!
�0 �#$ �! %
�&0
�&�
�&'
Figura 3.7 - Representação gráfica do critério das áreas iguais
Na figura anterior identificam-se:
28
𝑃𝑃𝑚 - potência mecânica 𝑃𝑃𝑒0 - potência eléctrica no período pré-defeito 𝑃𝑃𝑒𝐷 - potência eléctrica no período de defeito 𝑃𝑃𝑒𝑃 - potência eléctrica no período pós-defeito 𝛿𝛿0 - ângulo da máquina na situação pré-defeito 𝛿𝛿𝑐𝑙 - ângulo da máquina no instante de actuação das protecções 𝛿𝛿𝑚𝑎𝑥 - ângulo da máquina máximo atingido durante a oscilação 𝐴𝐴𝑎 - área de aceleração 𝐴𝐴𝑑 - área de desaceleração.
Inicialmente verifica-se 𝛿𝛿 = 𝛿𝛿0. Considerando que a perturbação é eliminada no instante
𝛿𝛿𝑐𝑙, podem definir-se duas áreas:
• área de aceleração cujo valor é dado por
𝐴𝐴𝑎 = � (𝑃𝑃𝑚 − 𝑃𝑃𝑒𝐷 𝑚𝑎𝑥 sin 𝛿𝛿)𝑑𝑑𝛿𝛿
𝛿𝑐𝑙
𝛿0 (3.18)
• área de desaceleração cujo valor é dado por
𝐴𝐴𝑑 = � (𝑃𝑃𝑒𝑃 𝑚𝑎𝑥 sin 𝛿𝛿 − 𝑃𝑃𝑚)𝑑𝑑𝛿𝛿
𝛿𝑚𝑎𝑥
𝛿𝑐𝑙 (3.19)
O sistema será estável quando 𝐴𝐴𝑑 ≥ 𝐴𝐴𝑎.
Como referido anteriormente, este método é aplicado a um sistema constituído por uma
máquina síncrona ligada a um barramento de potência de curto-circuito infinita. Para o aplicar a
um sistema de maiores dimensões (mais geradores) é necessário reduzir esse sistema a uma
máquina síncrona equivalente ligada a um barramento de potência de curto-circuito infinita (8).
Este método, denominado Método das Áreas Iguais Generalizado, utiliza o critério das áreas
iguais conjugado com as seguintes hipóteses:
a) sempre que ocorre a perda de sincronismo de um sistema multimáquina há uma
separação das máquinas em dois grupos;
b) a estabilidade pode ser avaliada substituindo as máquinas de cada grupo pelo
seu centro de inércia parcial;
c) a evolução temporal das duas máquinas resultantes da hipótese anterior pode
ser descrita através de séries de Taylor devidamente truncadas;
Considerando as hipóteses anteriores o método pode dividir-se nas seguintes etapas:
1. para uma determinada perturbação, o sistema multimáquina é dividido em dois
grupos: o grupo das máquinas críticas e o grupo das restantes máquinas;
2. cada grupo é reduzido a uma máquina equivalente usando o respectivo centro
de inércia parcial;
3. reduzir as duas máquinas equivalentes ao caso de uma máquina ligada a um
barramento infinito;
29
4. aplicar o critério das áreas iguais ao sistema obtido na etapa anterior, o que
permite calcular o tempo crítico de actuação das protecções e a margem de
estabilidade transitória, definida como
𝜂𝜂 = 𝐴𝐴𝑑 − 𝐴𝐴𝑎 (3.20)
5. usar as séries de Taylor convenientemente truncadas para obter expressões que
permitam determinar as medidas de estabilidade referidas anteriormente.
Exemplo
Apresenta-se agora um exemplo deste método, aplicado ao mesmo sistema utilizado no
exemplo anterior, o qual é sujeito à mesma perturbação.
Para começar, é necessário calcular as variáveis do sistema no período pré-defeito, mais
concretamente a força electromotriz e a potência mecânica:
𝑃𝑃𝑚 = 0,8 pu
𝑰𝑰 = �𝑺𝟑𝑽𝑽𝟑�∗
= 0,8 − 𝑗𝑗0,4 pu
𝑬𝑬𝟏𝟏 = 𝑽𝑽𝟑 + 𝑗𝑗(𝑋𝑋𝑑′ + 𝑋𝑋𝑐𝑐 + 𝑋𝑋𝑙)𝑰𝑰 = 1,323𝑒𝑒𝑗0,364 ⇒ �𝐸1 = 1,323 pu𝛿𝛿0 = 0,364 rad
�
Durante o período de defeito, a linha está desligada e então:
𝑃𝑃𝑒𝐷 = 0
Após a eliminação do defeito, a linha está de novo em funcionamento. Tem-se então no
período pós-defeito:
𝑃𝑃𝑒𝑃 = 𝑃𝑃𝑒𝑃 𝑚𝑎𝑥 sin𝛿𝛿 =𝐸1𝑉3
𝑋𝑋𝑑′ + 𝑋𝑋𝑐𝑐 + 𝑋𝑋𝑙sin𝛿𝛿
As diferentes curvas de potência, bem como as áreas de aceleração e desaceleração
podem ser observadas na Figura 3.8.
Aplicando então o critério das áreas iguais tem-se, para a situação limite:
𝐴𝐴𝑎 = 𝐴𝐴𝑑
� (𝑃𝑃𝑚 − 0)𝑑𝑑𝛿𝛿𝛿𝑐𝑙
𝛿0= � (𝑃𝑃𝑒𝑃 𝑚𝑎𝑥 sin 𝛿𝛿 −𝑃𝑃𝑚)𝑑𝑑𝛿𝛿
𝛿𝑢
𝛿𝑐𝑙
𝑃𝑃𝑚(𝛿𝛿𝑐𝑙 − 𝛿𝛿0) = −𝑃𝑃𝑒𝑃 𝑚𝑎𝑥(cos 𝛿𝛿𝑢 − cos 𝛿𝛿𝑐𝑙) − 𝑃𝑃𝑚(𝛿𝛿𝑢 − 𝛿𝛿𝑐𝑙)
𝑃𝑃𝑚(𝛿𝛿𝑐𝑙 − 𝛿𝛿0 + 𝛿𝛿𝑢 − 𝛿𝛿𝑐𝑙) + 𝑃𝑃𝑒𝑃 𝑚𝑎𝑥 cos𝛿𝛿𝑢 = 𝑃𝑃𝑒𝑃 𝑚𝑎𝑥 cos 𝛿𝛿𝑐𝑙
𝛿𝛿𝑐𝑙 = 1,644 rad
30
�
��
��
��
��
���
��0
��0 �
Figura 3.8 - Representação gráfica da situação estudada no exemplo do Método das Áreas Iguais
Como referido anteriormente, no período de defeito �� = 0, ficando a equação de oscilação reduzida a:
������� = ��
Sendo a derivada igual a uma constante, a equação é facilmente integrável:
�������� = ��� ⇔����� = ���
admitindo que ���� �0� = 0. Para se obter a evolução temporal de � basta integrar novamente a
equação obtendo-se:
���� = 12����� + ��0�
Substituindo por valores obtém-se o valor de ���:
������ = ��� = 12������� + ��0� ⇔ ��� ≅ 253ms∎
3.3. Métodos Híbridos
Os Métodos Híbridos resultam de uma junção entre os Métodos de Integração Numérica
e os Métodos Directos, oferecendo a grande capacidade de modelização característica dos
Métodos de Integração Directa combinada com a velocidade característica dos Métodos Directos
(3). Além de permitirem, como referido, uma grande flexibilidade de modelização aliada à rapidez
de execução, os Métodos Híbridos permitem também, para uma determinada contingência,
31
identificar o conjunto de máquinas críticas, isto é, o conjunto de máquinas responsáveis pela
perda de sincronismo do sistema.
Os métodos híbridos podem ser divididos em dois grupos: uma formulação utiliza o
sistema multi-máquina completo; outra formulação utiliza um modelo equivalente do sistema
constituído por uma máquina ligada a um barramento de potência de curto-circuito infinita (3). De
seguida serão apresentados sumariamente dois Métodos Híbridos, um de cada grupo: um
método que se baseia na Função de Energia Transitória e que utiliza o sistema multi-máquina
sem qualquer redução; um método que se baseia no Método das Áreas Iguais e que utiliza um
modelo reduzido equivalente do sistema.
3.3.1. Método híbrido baseado na função de energia transitória
Como referido anteriormente, este método, baseado na função de energia transitória,
utiliza o modelo multi-máquina completo.
O algoritmo desde método pode ser dividido em duas etapas distintas:
1. Período de defeito, durante o qual se efectua a integração numérica das
equações que descrevem o sistema, equações (3.4), e no fim do qual se calcula
o ponto de equilíbrio estável da configuração pós-defeito;
2. Período pós-defeito, durante o qual se continua a efectuar a integração numérica
das equações que descrevem o sistema, equações (3.5), mas, determinando em
cada passo de integração, se ocorre o cruzamento da SLEP ou se o sistema
atinge um máximo local da energia potencial. No primeiro caso, a situação é
considerada instável enquanto que, no segundo, após verificadas determinadas
condições, o sistema pode ser considerado estável. Caso não seja possível
determinar a estabilidade ou instabilidade da situação em estudo a integração
numérica continua a efectuar-se até tal ser possível.
Este método permite o cálculo do tempo crítico de actuação das protecções. Este é
obtido realizando sucessivas simulações com diferentes estimativas do tempo de actuação das
protecções que tendem para o tempo crítico. Estas estimativas são efectuadas utilizando as
margens de energia transitória das duas últimas simulações, considerando que o valor desta
margem varia linearmente com o tempo de actuação das protecções. O valor da margem de
energia transitória de cada simulação é obtido através da energia cinética do sistema ou da
energia potencial, caso, respectivamente, se trate de uma situação instável, para o qual é
negativo, ou uma situação estável, para o qual é positivo. O valor onde a margem se anula
corresponde ao ponto que se pretende determinar.
3.3.2. Método híbrido baseado no método das áreas iguais
Este método, como mencionado anteriormente, utiliza um modelo equivalente reduzido
do sistema constituído por uma máquina, daí ser denominado de SIME (SIngle Machine
Equivalent).
32
Tal como no método anterior, também o algoritmo deste método pode ser divido em duas
etapas distintas:
1. Período durante o qual se efectua a integração numérica das equações que
descrevem o sistema, equações (3.4) no período de defeito e equações (3.6) no
período pós-defeito. A duração do período de integração é estipulada pelo
critério das áreas iguais.
2. O sistema é reduzido a um sistema equivalente constituído por uma máquina
ligada a um barramento de potência de curto-circuito infinita. Usando o modelo
equivalente é calculada a margem de estabilidade transitória, a partir da equação
(3.20), que permite aferir acerca da estabilidade do sistema. Caso a margem de
estabilidade transitória seja positiva o sistema é considerado estável, caso
contrário é considerado instável.
Este método permite também o cálculo do tempo crítico de actuação das protecções. Tal
como no caso anterior, efectuam-se múltiplas simulações com diferentes estimativas do tempo
de actuação das protecções mas, utilizando agora, as duas últimas margens de estabilidade
transitória. Considera-se que o valor desta margem varia linearmente com o tempo de actuação
das protecções e o ponto para o qual esta se anula corresponde ao ponto que se pretende
determinar. Quando se verificar que a diferença entre duas estimativas consecutivas é inferior a
uma determinada tolerância a simulação pode ser interrompida.
3.4. Técnicas de Inteligência Artificial
As Técnicas de Inteligência Artificial são a mais recente abordagem, das apresentadas
neste trabalho, ao estudo da estabilidade transitória do SEE. Estas técnicas são aproximações
probabilísticas e usam métodos de aprendizagem automática.
Como ponto comum entre todas as Técnicas de Inteligência Artificial destaca-se a
desvantagem da necessidade de realizar um grande número de simulações em tempo diferido,
um processo bastante moroso e pesado do ponto de vista computacional. Já em tempo real
apresentam as vantagens de apresentar uma elevada eficiência computacional, reduzindo
drasticamente o tempo de cálculo, e uma grande capacidade de interpretação dos fenómenos
em análise.
De seguida serão apresentadas mais algumas características referentes às técnicas
mais representativas desta classe: Reconhecimento de Formas, Redes Neuronais e Árvores de
Decisão (3).
3.4.1. Reconhecimento de Formas
Esta técnica parte do facto de que, com base em informação previamente adquirida
sobre o SEE, é possível tirar conclusões em tempo real acerca do funcionamento deste, mesmo
para novos pontos de funcionamento não testados.
33
Para implementar este método é necessário gerar um conjunto de treino, através de
múltiplas simulações do SEE. De seguida, cada situação simulada deste conjunto treino, é
classificada por um conjunto de variáveis, denominadas variáveis primárias. Estas variáveis são
posteriormente filtradas para eliminar informação redundante, reduzindo assim a dimensão do
vector de variáveis. Posteriormente, entra-se na última fase de implementação que consiste na
síntese de um classificador que, com base no vector de variáveis que caracteriza cada ponto de
operação do SEE, seja capaz de o classificar em estável ou instável.
3.4.2. Redes Neuronais
As Redes Neuronais, tentam representar por meios computacionais, o funcionamento do
cérebro, consistindo em várias unidades de processamento, os neurónios, com um forte grau de
ligação entre elas. Esta metodologia pode ser usada como classificador de um determinado
ponto de operação do SEE devido à sua capacidade de reproduzir complexas relações
numéricas entre as variáveis recolhidas e, devido ao facto de a partir de dados recolhidos no
período pré-defeito, ser possível inferir acerca da estabilidade ou instabilidade do sistema para
uma determinada perturbação.
Tal como no método anterior, é necessário gerar inicialmente um conjunto de treino, do
qual são seleccionadas e filtradas as variáveis mais representativas, que permitem calibrar o
classificador, gerando um conjunto de regras que vão permitir classificar em tempo real cada
ponto de operação do SEE.
3.4.3. Árvores de Decisão
As Árvores de Decisão são outras das Técnicas de Inteligência Artificial aplicadas ao
estudo da estabilidade transitória do SEE. Na estrutura criada, uma árvore, os nós podem ser
considerados como estados, neste caso estável ou instável, e os ramos são o conjunto de
factores que levam a esses estados, neste caso as variáveis que descrevem o ponto de
funcionamento do SEE.
Para criar uma árvore de decisão é necessário, inicialmente, gerar uma base de dados,
de grandes dimensões, que represente os diferentes pontos de funcionamento do SEE. De
seguida procede-se à construção da árvore propriamente dita escolhendo, para cada nó, um
conjunto de atributos e um valor limite desses atributos que conduzam aos diferentes nós
seguintes. Na Figura 3.9, pode ser observado um caso simplificado com apenas um atributo e
apenas dois nós seguintes possíveis, sendo o atributo 𝑃𝑃𝐺 e o valor limite 𝑃𝑃𝐺lim.
1
𝑃𝑃𝐺𝐺 < 𝑃𝑃𝐺𝐺lim
2 3
Sim Não
Figura 3.9 - Exemplo de uma árvore de decisão
34
Durante a construção da árvore deve-se ter em atenção o compromisso entre a
complexidade e a fiabilidade para definir um ponto de paragem. Este método pode também
analisar novos estados não observados, a partir da base de dados.
3.5. Comparação dos diferentes métodos
Durante a exposição feita de cada um dos métodos foram logo referidas, na altura,
algumas das características principais de cada um deles. Efectua-se agora uma comparação
entre algumas características comuns aos vários métodos como as possibilidades de
modelização, as tarefas que precedem a sua utilização, os requisitos computacionais em tempo
real e a capacidade de fornecer medidas de controlo, a qual pode ser observada na Tabela 3.1.
Tabela 3.1 - Comparação entre os diferentes métodos
Métodos
Características
Integração Numérica
Métodos Directos
Métodos Híbridos
Técnicas de Inteligência
Artificial
Possibilidades de modelização Pormenorizada
Modelo clássico ou
pouco pormenorizada
Pormenorizada Pormenorizada
Tarefas preparatórias
(tempo deferido) Validação Validação
exaustiva Validação Geração da
base de dados e validação
Requisitos computacionais
(tempo real)
Grande esforço
Reduzido esforço
Esforço relativo
Esforço extremamente
reduzido
Medidas de controlo Não fornece Preventivo Preventivo Preventivo e
correctivo
Relativamente às possibilidades de modelização foi já referido que, com a excepção dos
Métodos Directos que usam um modelo clássico ou pouco pormenorizado, todos os outros
métodos permitem uma modelação pormenorizada.
Quanto às tarefas que precedem a sua utilização os Métodos de Integração Numérica e
Híbridos apenas necessitam de ser validados. Já os Métodos Directos precisam de uma
validação exaustiva devido às suposições e simplificações efectuadas. As Técnicas de
Inteligência Artificial para além da validação necessitam da geração de uma base de dados, que
é uma tarefa muito pesada computacionalmente.
Em tempo real assiste-se a uma inversão da carga de trabalho necessária de cada um
dos métodos, em comparação com as tarefas preparatórias referidas anteriormente, pois as
Técnicas de Inteligência Artificial e os Métodos Directos exigem um reduzido esforço de cálculo.
35
Os Métodos Híbridos e Directos permitem implementar medidas de controlo preventivo,
enquanto que as Técnicas de Inteligência Artificial permitem estabelecer medidas de controlo
preventivo mas também correctivo.
Os métodos que apresentam melhores resultados são os Métodos Híbridos e as
Técnicas de Inteligência Artificial mas, por estas últimas necessitarem da geração de uma base
de dados, tarefa morosa e que dá origem a uma grande quantidade de informação que necessita
de ser comprimida e armazenada, os Métodos Híbridos são actualmente os mais usados.
3.6. Conclusões
Neste capítulo foram apresentados os métodos mais relevantes para analisar a
estabilidade transitória de um SEE: Métodos de Integração Numérica, Métodos Directos,
Métodos Híbridos e Técnicas de Inteligência Artificial. Para cada um deles foram apresentadas
as principais características e, para alguns dos Métodos Directos, foram apresentados exemplos
de aplicação.
Por fim realizou-se uma comparação entre os diferentes métodos apresentados
relativamente a características como possibilidades de modelização, tarefas que precedem a sua
utilização, requisitos computacionais em tempo real e capacidade de ordenar contingências e
fornecer medidas de controlo. Constatou-se que os métodos que apresentam melhores
resultados são os Métodos Híbridos e as Técnicas de Inteligência Artificial, dando-se preferência
aos primeiros por serem mais fáceis de implementar pois não necessitam da geração de uma
base de dados.
37
4. Método híbrido implementado
Neste trabalho pretendeu-se implementar um método híbrido, em linguagem matlab,
capaz de, para diferentes contingências, calcular os tempos críticos de actuação das protecções.
As contingências simuladas são curto-circuitos trifásicos simétricos francos nas linhas junto aos
barramentos.
O método implementado conjuga os Métodos de Integração Numérica com o Método das
Áreas Iguais. O processo de integração numérica é interrompido antes do período total de
simulação utilizando índices de detecção de instabilidade/estabilidade, sendo de seguida o
sistema reduzido e estudado usando o Método das Áreas Iguais. Os valores dos tempos críticos
são obtidos através de um processo simples de regressão linear utilizando as margens de
estabilidade transitória.
O programa pode ser dividido em três etapas (Figura 4.1): uma primeira etapa que inclui
a introdução dos dados e o cálculo dos valores pré-defeito, a segunda etapa que inclui o ciclo
para determinação dos tempos críticos e, por fim, a última etapa que corresponde ao
armazenamento dos resultados.
Início FimProcedimentos preliminares
Ciclo para determinação do
tempo crítico Armazenamento dos resultados
Figura 4.1 - Estrutura geral do método implementado
De seguida será feita uma apresentação detalhada do modo de funcionamento do
programa, a qual está dividida pelas diferentes etapas referidas anteriormente. No fim do capítulo
é feita a ponte entre as diferentes etapas descritas e os diferentes ficheiros que contêm o código
escrito para implementar o método híbrido.
4.1. Introdução dos dados
Os dados são fornecidos ao programa através de ficheiros Excel, sendo necessário
fornecer os dados relativos à rede e aos defeitos que se pretendem simular. O próprio programa
é responsável pela sua importação para o espaço de variáveis. Os dados fornecidos estão
distribuídos por diferentes ficheiros Excel, existindo um ficheiro para cada uma das componentes
do SEE (geradores, linhas, transformadores), um ficheiro para os dados relativos aos
barramentos da rede e ainda um ficheiro para os defeitos a simular. Cada ficheiro contém uma
lista numerada com os vários elementos de cada conjunto. Todos os ficheiros referentes a uma
rede estão colocados na mesma pasta cujo caminho é comunicado ao programa, sendo este
autónomo até ao fim da simulação. A estrutura detalhada de cada ficheiro é agora apresentada.
38
4.1.1. Barramentos
Os dados relativos aos barramentos estão contidos num ficheiro denominado
“barramentos.xlsx”. Este ficheiro contém:
• numeração dos barramentos;
• tipo do barramento (referência, PV ou PQ);
• módulo da tensão para os barramentos PV e de referência;
• argumento da tensão para o barramento de referência;
• potências geradas activa para os barramentos PV e PQ e reactiva para os PQ;
• potência activa e reactiva de carga para todos os barramentos.
O ficheiro Excel tem a estrutura da Tabela 4.1.
Tabela 4.1 - Estrutura do ficheiro com os dados dos barramentos
Nº do barramento
Tipo do barramento
Módulo da
tensão
Argumento da tensão
Potência activa gerada
Potência reactiva gerada
Potência activa de
carga
Potência reactiva de carga
4.1.2. Geradores
Os dados relativos aos geradores estão contidos num ficheiro denominado
“geradores.xlsx”. Este ficheiro contém:
• numeração dos geradores;
• barramento ao qual o gerador está ligado;
• constante de inércia 𝐻;
• reactância transitória 𝑋𝑋𝑑′ ;
• coeficiente de amortecimento 𝐷𝐷.
O ficheiro Excel tem a estrutura da Tabela 4.2.
Tabela 4.2 - Estrutura do ficheiro com os dados dos geradores
Nº do gerador
Barramento de ligação
Constante de inércia 𝐻
Reactância transitória 𝑋𝑋𝑑′
Coeficiente de amortecimento 𝐷𝐷
4.1.3. Linhas de Transmissão
Os dados relativos às linhas estão contidos num ficheiro denominado “linhas.xlsx”. Este
ficheiro contém:
• numeração das linhas;
• barramentos entre os quais a linha está ligada;
• resistência e reactância longitudinal e susceptância transversal.
O ficheiro Excel tem a estrutura da Tabela 4.3.
39
Tabela 4.3 - Estrutura do ficheiro com os dados das linhas
Nº da linha
Barramento 1
Barramento 2
Resistência longitudinal
Reactância longitudinal
Susceptância transversal
4.1.4. Transformadores
Os dados relativos aos transformadores estão contidos num ficheiro denominado
“transformadores.xlsx”. Este ficheiro contém:
• numeração dos transformadores;
• barramento onde está ligado o primário;
• barramento onde está ligado o secundário;
• resistência e reactância de curto-circuito;
• relação de transformação, para os transformadores com tomadas e
desfasadores (unitária para os transformadores simples).
O ficheiro Excel tem a estrutura da Tabela 4.4.
Tabela 4.4 - Estrutura do ficheiro com os dados dos transformadores
Nº do transformador
Barramento do primário
Barramento do secundário Resistência Reactância Relação de
transformação
4.1.5. Contingências
Os dados relativos às contingências a simular estão contidos num ficheiro denominado
“contingencias.xlsx”. Este ficheiro contém:
• numeração das contingências;
• barramento junto do qual ocorre a contingência;
• linha que sai de serviço no período pós-defeito (0 para nenhuma).
O ficheiro Excel tem a estrutura da Tabela 4.5.
Tabela 4.5 - Estrutura do ficheiro com os dados dos defeitos
Nº da contingência Barramento Linha retirada de serviço
4.2. Cálculo dos valores pré-defeito
Esta fase do programa inclui a resolução de um trânsito de energia, a transformação das
cargas em impedâncias constantes e o cálculo dos valores iniciais das forças electromotrizes e
das potências mecânicas dos geradores.
O trânsito de potências é efectuado utilizando o Método de Newton-Raphson ficando
então disponíveis os módulos e argumentos da tensão para todos os barramentos da rede.
40
Relativamente à transformação das cargas em impedâncias constantes, esta é feita
tendo em conta o exposto no Capítulo 2.
Como já referido, a modelização foi feita tendo em conta o modelo clássico, o que leva a
que se considere que o módulo da força electromotriz e a potência mecânica são constantes,
permanecendo inalteráveis durante o período de simulação. Nesta etapa inicial do método
híbrido implementado são então, desde já, calculados estes valores.
Esta fase apenas necessita de ser efectuada uma vez, ficando então os valores
disponíveis para todos os defeitos que se pretendam simular.
4.3. Ciclo para determinação do tempo crítico
Após a realização dos procedimentos preliminares dá-se início ao ciclo para determinar o
tempo crítico de actuação das protecções. Esta etapa é a mais complexa do método
implementado e a que exige maior esforço computacional. Segue-se uma apresentação mais
detalhada deste ciclo cabendo de seguida, para os pontos considerados mais importantes,
secções próprias.
Como referido anteriormente, e como se pode observar na Figura 4.2, este ciclo é
iniciado com a integração numérica das equações que descrevem o sistema. É utilizado o
modelo clássico para descrever o sistema e o método de integração adoptado é o método
trapezoidal. Em cada passo de integração são calculados três índices: um de estabilidade (IDE),
um de instabilidade (IDCS) e um último que indica o tempo de observação para identificação do
conjunto de máquinas críticas (IDTO). A simulação continua até estes índices apresentarem
determinados comportamentos indicando que a simulação pode terminar.
Início
Fim
Integraçação numérica
IDEIDCSIDTO
Identificação das máquina críticas
Redução a uma máquina
equivalente
Cálculo da margem de estabildade
transitória
Nova estimativa do tempo crítico
∆𝑡𝑡𝑐𝑐𝑐𝑐 < 𝜀𝜀
Sim
Não
Figura 4.2 - Diagrama do ciclo para o cálculo do tempo crítico
Terminado o processo de integração numérica procede-se à identificação do conjunto
das máquinas críticas. Esta identificação é feita recorrendo a dois critérios de identificação, mas
41
apenas sendo utilizado um de cada vez. O critério a utilizar é escolhido tendo em conta a
estabilidade ou instabilidade da situação analisada, a qual é determinada pelos índices
calculados durante o processo de integração numérica.
Após identificadas as máquinas críticas pode-se então reduzir o sistema a um sistema
equivalente constituído por uma máquina ligada a um barramento de potência de curto-circuito
infinita. Este processo é feito considerando dois conjuntos: o conjunto das máquinas críticas e o
conjunto constituído pelas restantes máquinas. Começa-se por reduzir cada um destes conjuntos
a um conjunto equivalente para depois se obter o modelo reduzido do sistema.
Obtido o modelo equivalente do sistema calcula-se o valor da margem de estabilidade
transitória, necessário para estimar o valor do tempo crítico. Dado que o processo de integração
numérica é interrompido, devido aos índices utilizados, antes do tempo total de simulação, não é
possível obter a curva de potência da máquina equivalente para todos os valores do ângulo de
potência equivalente 𝛿𝛿𝑒𝑞. É então necessário usar métodos para estimar os valores em falta,
partindo dos valores já conhecidos, sendo então utilizados dois métodos: aproximação por uma
função trigonométrica e aproximação por um polinómio de segundo grau. A opção pelo método
de aproximação a utilizar é tomada tendo em conta a proximidade ao valor final do tempo crítico
a determinar.
Após obtida a margem de estabilidade para determinada iteração, faz-se uma nova
estimativa do valor do tempo crítico utilizando a margem de estabilidade calculada nesta iteração
e a margem calculada na iteração anterior. Esta estimativa baseia-se no facto de se considerar
que o valor da margem de estabilidade 𝜂𝜂 varia linearmente com o tempo de actuação das
protecções 𝑡𝑡𝑐𝑙.
Segue-se agora uma exposição dos pontos considerados mais importantes.
4.3.1. Índices de estabilidade e instabilidade
Os métodos híbridos caracterizam-se por interromperem a integração das equações
diferenciais que descrevem o comportamento do sistema antes do tempo total de simulação.
Para determinar o momento em que a integração pode ser interrompida usam-se índices que,
com base nas variáveis do sistema, permitem chegar a uma conclusão acerca da estabilidade do
sistema antes do tempo total de simulação ser atingido.
No método implementado são utilizados dois índices para determinar a estabilidade do
sistema, sendo o valor destes calculado em cada iteração do período pós-defeito. Como referido
anteriormente existe ainda um terceiro índice mas, por estar relacionado com a identificação das
máquinas críticas, será apresentado posteriormente.
O primeiro índice apresentado, denominado IDCS, permite determinar a instabilidade do
sistema (3). O seu valor é calculado por
42
𝐼𝐷𝐷𝐶𝐶𝑆 = �𝑓𝑖2
𝑚
𝑖=1
(4.1)
onde 𝑓𝑖 é a potência de aceleração da máquina 𝑑𝑑, dada pela equação (2.13). Monitorizando o
valor de IDCS é possível declarar a instabilidade do sistema quando este índice atingir um
mínimo relativo.
O segundo índice, denominado IDE e que permite determinar a estabilidade do sistema
(3), é determinado a partir da expressão
𝐼𝐷𝐷𝐸 = �𝜔�𝑖�𝜃𝑖 − 𝜃𝑖𝑐𝑙�
𝑚
𝑖=1
(4.2)
onde 𝜔�𝑖 - velocidade angular da máquina 𝑑𝑑, referida ao centro de inércia do sistema 𝜃𝑖 - ângulo rotórico da máquina 𝑑𝑑, referido ao centro de inércia do sistema 𝜃𝑖𝑐𝑙- ângulo rotórico da máquina 𝑑𝑑 no instante de eliminação do defeito, referido ao centro de inércia do sistema.
Na forma matricial este índice calcula-se através da equação
𝐼𝐷𝐷𝐸 = [𝜔�]𝑇([𝜃] − [𝜃𝑐𝑙]) (4.3)
O valor deste índice é monitorizado em cada passo de integração até se observar uma
mudança de sinal, de valor positivo para negativo, o que permite concluir que o sistema é
estável. Refira-se que é necessário acontecer uma mudança de valor positivo para negativo, ou
seja, se este índice apresentar, logo na primeira iteração em que é calculado, um valor negativo,
o sistema não é considerado estável e o processo de integração prossegue.
A interacção destes índices com o processo de integração numérica é facilmente
explicável. Como referido, em cada passo de integração, no período pós-defeito, o seu valor é
calculado. O primeiro índice a apresentar um comportamento que permita tirar conclusões
(mínimo para IDCS, mudança de sinal para IDE) é o índice válido para essa situação e o
processo de integração numérica é interrompido de imediato (9).
Nas figuras seguintes pode ser observado o que acabou de ser explicado. Na Figura 4.3,
obtida utilizando um tempo de actuação das protecções de 𝑡𝑡𝑐𝑙 = 0,4 s pode-se observar que o
índice IDCS atinge primeiro um mínimo sendo então o processo de integração interrompido e a
situação considerada instável. Na Figura 4.4, onde o tempo de actuação das protecções definido
foi de 𝑡𝑡𝑐𝑙 = 0,3 s a situação já é considerada estável pois, como se observa, o índice IDE sofre
primeiro uma mudança de sinal.
43
Figura 4.3 - Comportamento dos índices IDCS e IDE, situação instável
Figura 4.4 - Comportamento dos índices IDCS e IDE, situação estável
4.3.2. Identificação das máquinas críticas
A identificação do conjunto de máquinas críticas é um passo fundamental no método
implementado. Esta identificação é feita a partir de dois critérios de ordenação das máquinas
críticas sendo, para cada situação, apenas utilizado um critério, o qual é escolhido tendo em
conta a informação da estabilidade fornecida pelos índices IDCS e IDE apresentados
anteriormente. A utilização dos critérios de ordenação implica ainda a determinação do tempo de
observação óptimo das variáveis que caracterizam o estado dos geradores, sendo os critérios
avaliados para este instante.
O tempo de observação ideal para a identificação do conjunto de máquinas críticas é
determinado a partir de um índice baseado nas variáveis que descrevem os geradores (3). Este
índice, denominado IDTO, é calculado a partir da expressão:
𝐼𝐷𝐷𝑇𝑇𝑂 = �𝑓𝑖𝜔�𝑖
𝑚
𝑖=1
(4.4)
onde 𝑓𝑖 - potência de aceleração da máquina 𝑑𝑑, dada pela equação (2.13) 𝜔�𝑖 - velocidade angular da máquina 𝑑𝑑, referida ao centro de inércia do sistema.
O índice IDTO é calculado em cada passo do processo de integração das equações que
descrevem o sistema, no período pós-defeito. O instante em que este índice apresenta uma
44
mudança de sinal corresponde ao instante ideal para a determinação do conjunto de máquinas
críticas (Figura 4.5).
Figura 4.5 - Comportamento do índice IDTO
Após estar determinado o instante ideal para a identificação do conjunto das máquinas
críticas procede-se à ordenação das máquinas segundo o critério escolhido. Quando a simulação
é estável utiliza-se um índice baseado na variação incremental dos ângulos rotóricos. Caso
contrário, ou seja, quando a simulação é instável, utiliza-se um método denominado ordenação
das máquinas críticas (CMR – Critical Machines Ranking).
O índice correspondente ao primeiro critério (3), baseado na variação incremental dos
ângulos rotóricos, é calculado em relação ao centro de inércia do sistema, para cada máquina do
sistema, pela expressão
𝐼𝐴𝐴𝐶𝑂𝐼𝑖 = � |𝜃𝑖(𝑡𝑡 + ∆𝑡𝑡) − 𝜃𝑖(𝑡𝑡)|
𝑡𝑜𝑏𝑠
𝑡=𝑡𝑐𝑙
(4.5)
onde 𝜃𝑖 - ângulo rotórico da máquina 𝑑𝑑, referido ao centro de inércia do sistema ∆𝑡𝑡 – passo do processo de integração.
De seguida os índices para todas as máquinas são ordenados por ordem decrescente do
seu valor, calculando-se, após a ordenação, a maior diferença entre dois valores consecutivos. O
conjunto de máquinas situado acima do maior intervalo constitui o conjunto de máquinas crítico.
O critério utilizado quando a simulação é instável, o CMR, baseia-se no facto de
considerar que o grau de criticalidade de uma máquina é directamente proporcional à amplitude
do seu ângulo rotórico, observado num instante apropriado (10). A aplicação deste método
inicia-se ordenando, no instante determinado pelo índice IDTO, as máquinas do sistema por
ordem decrescente em função do seu ângulo rotórico 𝛿𝛿𝑖. De seguida é necessário considerar
vários conjuntos de máquinas críticas determinando-se qual o mais desfavorável. Para isso,
partindo-se do conjunto já ordenado, procede-se da seguinte forma:
45
1. forma-se um conjunto constituído pelas primeiras 𝑑𝑑 máquinas da lista ordenada,
em que 𝑑𝑑 = 1, 2, … ,𝑚𝑚− 1, sendo 𝑚𝑚 o número total de geradores;
2. para o primeiro conjunto, considerado o das máquinas críticas, calcula-se o
ângulo rotórico equivalente procedendo-se da mesma forma para o conjunto das
restantes máquinas. Calcula-se a diferença entre estes dois valores;
3. repetir os passos 1 e 2 até que as primeiras 𝑚𝑚− 1 máquinas tenham sido
incluídas no conjunto das máquinas críticas.
Os dois conjuntos estabelecidos que correspondam à maior diferença entre ângulos
equivalentes correspondem ao conjunto das máquinas críticas e ao conjunto das restantes
máquinas.
O ângulo rotórico equivalente para cada conjunto, referido no passo 2, é calculado
recorrendo às expressões (2.7) e (2.8). Sendo C o conjunto das máquinas críticas e R o conjunto
das restantes máquinas tem-se para o conjunto C:
𝛿𝛿𝐶 =1𝑀𝐶
�𝑀𝑘𝛿𝛿𝑘𝑘∈𝐶
(4.6)
𝑀𝐶 = �𝑀𝑘𝑘∈𝐶
(4.7)
Enquanto que para o conjunto R se tem:
𝛿𝛿𝑅 =1𝑀𝑅
�𝑀𝑗𝛿𝛿𝑗𝑗∈𝑅
(4.8)
𝑀𝑅 = �𝑀𝑗𝑗∈𝑅
(4.9)
Nas expressões anteriores tem-se
𝛿𝛿𝐶 - ângulo rotórico equivalente do conjunto C 𝑀𝑘 - coeficiente de inércia da máquina 𝑘, pertencente ao conjunto C 𝛿𝛿𝑘 - ângulo rotórico da máquina 𝑘, pertencente ao conjunto C 𝑀𝐶 - coeficiente de inércia equivalente do conjunto C 𝛿𝛿𝑅 - ângulo rotórico equivalente do conjunto R 𝑀𝑗 - coeficiente de inércia da máquina 𝑗𝑗, pertencente ao conjunto R 𝛿𝛿𝑗 - ângulo rotórico da máquina 𝑗𝑗, pertencente ao conjunto R 𝑀𝑅 - coeficiente de inércia equivalente do conjunto R.
4.3.3. Redução do sistema a uma máquina equivalente
Para poder aplicar o Método das Áreas Iguais a um sistema multi-máquina é necessário
reduzir este último a um modelo equivalente constituído por uma máquina ligada a um
barramento infinito. Este redução é efectuada tendo em conta o conjunto de máquinas críticas
identificado no passo anterior. O processo inclui a representação do conjunto de máquinas
críticas por uma máquina equivalente, efectuando-se o mesmo processo para o conjunto das
46
restantes máquinas. De seguida estas duas máquinas são reduzidas a uma única máquina
equivalente, a qual descreve de forma correcta o sistema e permite uma análise facilitada
utilizando o Método das Áreas Iguais (3).
O primeiro passo, que consiste na redução do conjunto de máquinas críticas (conjunto C)
e do conjunto das restantes máquinas (conjunto R) a máquinas equivalentes, máquina C e
máquina R respectivamente, inicia-se utilizando as equações (4.7) e (4.9). De seguida, o ângulo
rotórico de cada máquina equivalente, em cada instante de tempo, é obtido utilizando equações
equivalentes às equações (4.6) e (4.8), procedendo-se de forma equivalente para as velocidades
angulares. Tem-se então para o ângulo rotórico do conjunto C:
𝛿𝛿𝐶(𝑡𝑡) =1𝑀𝐶
�𝑀𝑘𝛿𝛿𝑘(𝑡𝑡)𝑘∈𝐶
(4.10)
enquanto que para a velocidade angular se tem
𝜔𝐶(𝑡𝑡) =1𝑀𝐶
�𝑀𝑘𝜔𝑘(𝑡𝑡)𝑘∈𝐶
(4.11)
Para o conjunto R, por sua vez, tem-se para o ângulo rotórico
𝛿𝛿𝑅(𝑡𝑡) =1𝑀𝑅
�𝑀𝑗𝛿𝛿𝑗𝑗∈𝑅
(𝑡𝑡) (4.12)
enquanto que para a velocidade angular
𝜔𝑅(𝑡𝑡) =1𝑀𝑅
�𝑀𝑗𝜔𝑗(𝑡𝑡)𝑗∈𝑅
(4.13)
Nas expressões anteriores tem-se 𝜔𝐶 - velocidade angular equivalente do conjunto C 𝜔𝑘 - velocidade angular da máquina 𝑘, pertencente ao conjunto C 𝜔𝑅 - velocidade angular equivalente do conjunto R 𝜔𝑗 - velocidade angular da máquina 𝑗𝑗, pertencente ao conjunto R.
Após obtidos os modelos equivalentes para os dois conjuntos realiza-se o segundo
passo, o qual permite obter o modelo equivalente constituído por uma máquina ligada a um
barramento de potência de curto-circuito infinita. O ângulo rotórico equivalente é obtido através
de
𝛿𝛿𝑒𝑞(𝑡𝑡) = 𝛿𝛿𝐶(𝑡𝑡) − 𝛿𝛿𝑅(𝑡𝑡) (4.14)
enquanto que a velocidade angular equivalente se obtém por
𝜔𝑒𝑞(𝑡𝑡) = 𝜔𝐶(𝑡𝑡) − 𝜔𝑅(𝑡𝑡) (4.15)
47
O coeficiente de inércia total obtém-se com base nos coeficientes de inércia calculados
através das expressões (4.7) e (4.9):
𝑀𝑇 = 𝑀𝐶 + 𝑀𝑅 (4.16)
O coeficiente de inércia do sistema equivalente, obtém-se também com base nos
coeficientes de inércia calculados através das expressões (4.7) e (4.9) e é dado por
𝑀𝑒𝑞 =𝑀𝐶𝑀𝑅
𝑀𝐶 + 𝑀𝑅 (4.17)
Para estarem determinadas todas as variáveis relativas à maquina equivalente falta
apenas determinar as potências eléctrica e mecânica equivalentes. Para a potência mecânica
equivalente tem-se
𝑃𝑃𝑚 𝑒𝑞(𝑡𝑡) =
1𝑀𝑇
�𝑀𝑅�𝑃𝑃𝑚𝑘(𝑡𝑡)𝑘∈𝐶
− 𝑀𝐶�𝑃𝑃𝑚𝑗(𝑡𝑡)𝑗∈𝑅
� (4.18)
enquanto que a potência eléctrica equivalente é obtida através de
𝑃𝑃𝑒 𝑒𝑞(𝑡𝑡) =
1𝑀𝑇
�𝑀𝑅�𝑃𝑃𝑒𝑘(𝑡𝑡)𝑘∈𝐶
− 𝑀𝐶�𝑃𝑃𝑒𝑗(𝑡𝑡)𝑗∈𝑅
� (4.19)
Nas expressões anteriores tem-se 𝑃𝑃𝑚𝑘 - potência mecânica da máquina 𝑘, pertencente ao conjunto C 𝑃𝑃𝑚𝑗 - potência mecânica da máquina 𝑗𝑗, pertencente ao conjunto R 𝑃𝑃𝑒𝑘 - potência eléctrica da máquina 𝑘, pertencente ao conjunto C 𝑃𝑃𝑒𝑗 - potência eléctrica da máquina 𝑗𝑗, pertencente ao conjunto R.
Dispõe-se agora de um modelo equivalente que será usado para calcular a margem de
estabilidade transitória do sistema.
4.3.4. Cálculo da margem de estabilidade transitória
O cálculo da margem de estabilidade transitória é importante, no método implementado,
pois é através do seu valor que são estimados os valores de tempo crítico. Este cálculo, tendo
disponível as curvas de potência, é relativamente fácil, sendo efectuado através da expressão
𝜂𝜂 = � �𝑃𝑃𝑒 𝑒𝑞 − 𝑃𝑃𝑚 𝑒𝑞�
𝛿𝑒𝑞𝑢
𝛿𝑒𝑞0
𝑑𝑑𝛿𝛿𝑒𝑞 (4.20)
na qual 𝑃𝑃𝑒 𝑒𝑞 - potência eléctrica fornecida pela máquina equivalente 𝑃𝑃𝑚 𝑒𝑞 - potência mecânica da máquina equivalente 𝛿𝛿𝑒𝑞0 - ângulo rotórico da máquina equivalente no instante pré-defeito 𝛿𝛿𝑒𝑞𝑢 - ângulo rotórico da máquina equivalente, de valor superior a 𝛿𝛿𝑒𝑞0 , para o qual a potência eléctrica volta a igualar a potência mecânica.
48
Acontece que no método implementado, devido à interrupção do processo de integração
numérica antes do tempo total de simulação, não se dispõe da totalidade da curva de potência
no período pós-defeito sendo necessário estimar os seus valores para depois utilizar a expressão
(4.20). Como já foi referido, esta estimativa é feita utilizando duas aproximações diferentes,
dependendo a opção por uma delas da proximidade ao valor final do tempo crítico que se
pretende determinar. Esta aproximação é verificada analisando os intervalos entre dois valores
consecutivos da margem de estabilidade transitória. Quanto mais próximo do valor final, menores
serão os intervalos. De seguida serão apresentadas as duas aproximações utilizadas e também
a razão porque se optou por utilizar duas aproximações.
A primeira aproximação utiliza uma função trigonométrica, ou seja, pretende-se
representar a curva de potência da máquina equivalente por uma função do tipo
𝑃𝑃𝑒 𝑒𝑞�𝛿𝛿𝑒𝑞� = 𝑃𝑃𝑒 𝑒𝑞𝑚𝑎𝑥 sin�𝛿𝛿𝑒𝑞� (4.21)
na qual 𝑃𝑃𝑒 𝑒𝑞�𝛿𝛿𝑒𝑞� - potência eléctrica fornecida pela máquina equivalente 𝑃𝑃𝑒 𝑒𝑞𝑚𝑎𝑥 - potência eléctrica máxima fornecida pela máquina equivalente
𝛿𝛿𝑒𝑞 - ângulo rotórico da máquina equivalente.
Utilizando um ponto conhecido da curva de potência (𝛿𝛿𝑒𝑞, 𝑃𝑃𝑒 𝑒𝑞) calcula-se o valor de
𝑃𝑃𝑒 𝑒𝑞𝑚𝑎𝑥, podendo-se então determinar o valor da potência eléctrica fornecida pela máquina
equivalente para qualquer valor do ângulo rotórico da máquina equivalente, 𝛿𝛿𝑒𝑞.
A segunda aproximação é uma aproximação polinomial. Esta aproximação utiliza um
polinómio do 2º grau, sendo a curva de potência da máquina equivalente aproximada por uma
expressão do tipo
𝑃𝑃𝑒 𝑒𝑞�𝛿𝛿𝑒𝑞� = 𝑐𝑐1𝛿𝛿𝑒𝑞2 + 𝑐𝑐2𝛿𝛿𝑒𝑞 + 𝑐𝑐3 (4.22)
na qual 𝑃𝑃𝑒 𝑒𝑞�𝛿𝛿𝑒𝑞� - potência eléctrica fornecida pela máquina equivalente 𝑐𝑐1, 𝑐𝑐2 e 𝑐𝑐3 - constantes da aproximação polinomial 𝛿𝛿𝑒𝑞 - ângulo rotórico da máquina equivalente.
As constantes da aproximação polinomial são calculadas através da função polyfit,
disponibilizada pelo matlab sendo depois utilizada a função polyval para determinar o valor do
polinómio, que corresponde à potência eléctrica fornecida, para os diferentes valores do ângulo
rotórico 𝛿𝛿𝑒𝑞.
A aproximação trigonométrica é utilizada, nas primeiras iterações, até que a variação
verificada nas sucessivas estimativas do valor do tempo crítico seja menor que uma tolerância
definida no início da simulação, a qual é um múltiplo 𝛼𝛼 da tolerância 𝜀𝜀 necessária para parar o
ciclo completo (Figura 4.6).
49
Início
Aproximação polinomial
Fim
Sim
Não
Aproximação trigonométrica
∆𝑡𝑡𝑐𝑐𝑐𝑐 < 𝛼𝛼 .𝜀𝜀 ∆𝑡𝑡𝑐𝑐𝑐𝑐 < 𝜀𝜀
Sim
Não
Figura 4.6 - Escolha do método de aproximação da curva de potência a utilizar
A ideia de utilizar duas aproximações para a curva de potência, durante o ciclo iterativo,
nasceu durante o desenvolvimento do código. Mais concretamente, verificou-se que a
aproximação polinomial conduzia a resultados mais próximos dos reais (verificados através de
programas de integração numérica) mas que era pouco robusta quando a estimativa inicial do
tempo crítico ficava afastada do valor real. Ao invés, a aproximação por uma função
trigonométrica era mais robusta relativamente à estimativa inicial mas apresentava resultados
menos exactos. Optou-se então por iniciar o processo com a aproximação trigonométrica,
utilizando a aproximação polinomial após algumas iterações.
4.3.5. Estimativa dos valores do tempo crítico
Como já foi referido, a estimativa dos valores do tempo crítico é efectuada com base nos
valores da margem de estabilidade transitória, considerando-se que existe uma relação linear
entre o tempo de actuação das protecções 𝑡𝑡𝑐𝑙 e o valor da margem de estabilidade transitória 𝜂𝜂,
como pode ser observado na Figura 4.7 (3).
𝑡𝑡𝑐𝑐𝑐𝑐 𝑡𝑡𝑐𝑐𝑐𝑐𝑑𝑑 𝑡𝑡𝑐𝑐𝑐𝑐𝑑𝑑−1
𝜂𝜂
𝜂𝜂𝑑𝑑
𝜂𝜂𝑑𝑑−1
Figura 4.7 - Relação linear entre 𝜼 e 𝒕𝒄𝒄𝒍
Utilizando os valores da margem de estabilidade transitória da iteração actual e da
iteração anterior (𝜂𝜂𝑖 e 𝜂𝜂𝑖−1), cujos tempos de actuação das protecções correspondentes (𝑡𝑡𝑐𝑙𝑖 e
𝑡𝑡𝑐𝑙𝑖−1) são conhecidos, utiliza-se um processo de regressão linear que permite obter a equação da
recta que descreve a relação entre o tempo de actuação e a margem de estabilidade transitória.
Esta equação é do tipo
50
𝜂𝜂(𝑡𝑡𝑐𝑙) = 𝑚𝑚𝜂𝑡𝑡𝑐𝑙 + 𝑏𝜂 (4.23)
sendo 𝑚𝑚𝜂 o declive dado por
𝑚𝑚𝜂 =
𝜂𝜂𝑖 − 𝜂𝜂𝑖−1
𝑡𝑡𝑐𝑙𝑖 − 𝑡𝑡𝑐𝑙𝑖−1 (4.24)
e 𝑏𝜂 a ordenada na origem calculada através da equação (4.23) utilizando um dos pares (𝑡𝑡𝑐𝑙, 𝜂𝜂)
disponíveis.
O ponto que se pretende determinar corresponde ao valor do tempo de actuação das
protecções para o qual a margem de estabilidade transitória é nula. Utilizando novamente a
equação (4.23), fazendo 𝜂𝜂(𝑡𝑡𝑐𝑙) = 0, determina-se a nova estimativa de 𝑡𝑡𝑐𝑟 que será utilizada na
iteração seguinte.
4.4. Resultados fornecidos
Como referido, o programa permite determinar os tempos críticos de actuação das
protecções para determinadas contingências, especificadas pelo utilizador no ficheiro de entrada.
Os resultados da simulação são guardados num ficheiro Excel que, para além dos já referidos
tempos críticos, inclui também a numeração das contingências e o número de iterações
realizadas pelo algoritmo responsável por determinar os tempos críticos. Este ficheiro é
denominado “resultados.xlxs” e é criado na mesma directoria dos outros ficheiros referentes à
rede em estudo. Este ficheiro tem a estrutura da Tabela 4.6.
Tabela 4.6 - Estrutura do ficheiro com os resultados da simulação
Nº da contingência Tempo crítico Número de iterações efectuadas
4.5. Ficheiros de código
Os algoritmos que constituem o método híbrido implementado estão distribuídos por
cinco ficheiros de código. Para correr o método é chamado um dos ficheiros que depois se
encarrega de chamar os outros ficheiros quando necessário.
O ficheiro a ser chamado é denominado “TC”. Neste ficheiro, para além de estarem
definidos pelo utilizador parâmetros como o passo de integração, tempo total de simulação,
tolerância, etc., são realizadas as acções preliminares e as finais. Nas acções preliminares,
inclui-se o cálculo de um trânsito de energia, cujo código está contido no ficheiro “TC_TE”.
A base do ciclo para determinar o tempo crítico está contida no ficheiro “TC”. Mais
concretamente, este ficheiro contém os algoritmos que estimam os novos valores do tempo
crítico baseados nos valores da margem de estabilidade e é o responsável por chamar os
restantes ficheiros que contêm as rotinas necessárias ao ciclo. Estes ficheiros incluem: o ficheiro
51
“TC_IN”, onde se processa a integração numérica das equações que descrevem o sistema; o
ficheiro “TC_MC” que contém o código para identificar o conjunto das máquinas críticas e
posterior redução a uma máquina equivalente; o ficheiro “TC_ME”, no qual é calculado o valor da
margem de estabilidade transitória.
A interacção que se procurou agora explicar pode ser observada graficamente na Figura
4.8, onde se pode observar o ficheiro principal TC, com as suas diferentes fases, nas quais são
chamados os programas secundários.
• Acções preliminares
• Ciclo para determinação do tempo crítico
• Acções finais
TC
TC_TE
TC_IN TC_MC TC_ME
Figura 4.8 - Interacção dos diferentes ficheiros de código no método híbrido implementado
4.6. Conclusões
Neste capítulo foi apresentado em detalhe o método híbrido implementado, mais
concretamente a sua estrutura geral e pontos fundamentais.
Relativamente à sua estrutura, referiu-se que o método pode ser dividido em três etapas
diferentes: procedimentos preliminares, ciclo para determinação do tempo crítico e
armazenamento dos resultados. Nos procedimentos preliminares são efectuadas a leitura dos
ficheiros de dados e o cálculo dos valores pré-defeito. O ciclo para determinação do tempo crítico
é a etapa mais complexa do método e é a responsável por obter os tempos críticos para as
diferentes contingências. Na última etapa apenas é efectuado armazenamento dos resultados
obtidos durante a etapa anterior.
De seguida foram descritos em detalhe os pontos fundamentais do método desenvolvido:
índices de estabilidade e instabilidade, identificação do conjunto das máquinas críticas, redução
do sistema a uma máquina equivalente, cálculo da margem de estabilidade transitória e
estimativa dos valores do tempo crítico. Quanto aos índices de estabilidade e instabilidade foram
apresentadas as expressões para o seu cálculo e qual a sua interacção com o programa. De
seguida, apresentaram-se os critérios utilizados na determinação do conjunto das máquinas
críticas tendo-se verificado que se usa um critério para situações estáveis e outro para situações
instáveis. A redução do sistema a uma máquina equivalente é efectuada com base num conjunto
de expressões matemáticas que foram apresentadas. Para o cálculo da margem de estabilidade
52
transitória é necessário efectuar uma aproximação para determinar os valores em falta da curva
de potência utilizando-se inicialmente uma aproximação trigonométrica e posteriormente uma
aproximação polinomial de grau 2. Os valores de tempo crítico são estimados considerando uma
relação linear entre o tempo de actuação das protecções e o valor da margem de estabilidade
transitória sendo esta relação descrita por uma recta cuja expressão é determinada por um
processo simples de regressão linear.
Por fim apresentaram-se os diferentes ficheiros pelos quais está distribuído o código
necessário para implementar o método híbrido descrito. Verificou-se que há um ficheiro principal
a partir do qual são chamadas as outras rotinas, contidas noutros ficheiros, necessárias para o
cálculo do valor do tempo crítico de actuação das protecções.
53
5. Resultados computacionais
Após implementado, o método híbrido foi testado para avaliar a validade dos resultados
fornecidos. Os testes efectuados consistiram na simulação de um variado conjunto de
perturbações, na rede de teste da CIGRE, com o objectivo de comparar os tempos críticos
obtidos através deste método com os fornecidos por programas que utilizam Métodos de
Integração Numérica (3).
Neste capítulo serão apresentados os resultados dos testes realizados, uma análise dos
erros e possíveis correcções, e ainda um exemplo de aplicação para uma perturbação específica
da rede.
5.1. Testes computacionais
Como já referido, foi utilizada a rede de teste da CIGRE, com 7 geradores e restantes
características disponíveis em anexo. As perturbações simuladas consistiram em curto-circuitos
trifásicos simétricos, junto aos barramentos, cuja eliminação acontece, após o período de defeito,
espontaneamente ou pela retirada de serviço de uma linha.
5.1.1. Resultados obtidos
Para cada perturbação pretendeu-se determinar o tempo crítico de actuação das
protecções. Os resultados obtidos bem como os fornecidos por programas que utilizam Métodos
de Integração Numérica são apresentados na Tabela 5.1. Na primeira coluna encontra-se o
número da perturbação. A segunda coluna contém o barramento junto ao qual se dá o defeito
enquanto que na terceira está a linha retirada de serviço. O tipo de instabilidade, denominado em
(3) como primeira oscilação (PO) ou oscilação inversa (OI), é listado na quarta coluna. Na quinta
e sexta colunas encontram-se os tempos críticos fornecidos, respectivamente, pelo método
híbrido e por Métodos de Integração Numérica. Por fim, nas sétima e oitava colunas, é
apresentada respectivamente, em valor absoluto e percentagem, a diferença entre os valores
obtidos pelos dois métodos, correspondendo 𝑡𝑡𝑐𝑟 ℎ𝑖𝑏 ao método híbrido e 𝑡𝑡𝑐𝑟 𝑖𝑛 aos métodos de
integração numérica.
A diferença em valor absoluto foi obtida por
∆𝑡𝑡𝑐𝑟𝑚𝑒𝑡 = 𝑡𝑡𝑐𝑟 ℎ𝑖𝑏 − 𝑡𝑡𝑐𝑟 𝑖𝑛 (5.1)
A diferença em percentagem foi calculada a partir da expressão
%∆𝑡𝑡𝑐𝑟𝑚𝑒𝑡 =
∆𝑡𝑡𝑐𝑟𝑚𝑒𝑡
𝑡𝑡𝑐𝑟 𝑖𝑛 (5.2)
54
Tabela 5.1 - Tempos críticos de actuação das protecções para a rede de teste da CIGRE
N.º contingência Barramento
Linha fora de serviço
Tipo de instabilidade
𝑡𝑡𝑐𝑟 método híbrido
[ms]
𝑡𝑡𝑐𝑟 integração numérica
[ms]
∆𝑡𝑡𝑐𝑟𝑚𝑒𝑡 [ms] %∆𝑡𝑡𝑐𝑟𝑚𝑒𝑡
1 1 - PO 356 356 0 0 2 1 1 PO 345 345 0 0 3 1 2 PO 348 347 1 0,29 4 2 - PO 412 412 0 0 5 2 3 OI 401 369 32 8,67 6 2 4 OI 412 386 26 6,74 7 3 - PO 390 390 0 0 8 3 1 PO 379 379 0 0 9 3 3 PO 391 391 0 0 10 3 5 PO 387 387 0 0 11 3 6 PO 390 390 0 0 12 4 - PO 496 496 0 0 13 4 2 PO 482 482 0 0 14 4 5 PO 490 489 1 0,20 15 4 7 PO 0 0 0* 0* 16 4 8 OI 531 458 73 15,94 17 4 9 PO 494 497 -3 -0,60 18 4 10 PO 494 493 1 0,20 19 5 - PO 353 353 0 0 20 5 7 PO 0 0 0 0 21 6 - PO 486 487 -1 -0,21 22 6 8 PO ? 437 ? ? 23 6 11 PO 479 482 -3 -0,62 24 7 - PO 340 341 -1* -0,29* 25 7 12 PO 0 0 0 0 26 8 - PO 480 480 0 0 27 8 12 PO 0 0 0 0 28 8 11 PO 441 441 0 0 29 8 13 PO 478 478 0 0
Analisando os resultados conclui-se que, das 29 situações simuladas, 23 conduziram a
resultados concordantes com os obtidos por Métodos de Integração Numérica. Verifica-se que o
erro absoluto máximo observável é de 3 ms, o que corresponde a 0,62%. As restantes 6
situações, assinaladas a negrito na tabela, serão de seguida analisadas com mais pormenor.
5.1.2. Análise dos erros
As 6 situações analisadas marcadas como erradas podem ser agrupadas e estudadas
em conjunto: as perturbações 5, 6 e 16 são do tipo oscilação inversa; a perturbação 22 não
conduziu a qualquer resultado; as situações 15 e 24 apresentam valores concordantes com os
reais mas estes não foram obtidos de forma autónoma.
• Perturbações 5, 6 e 16
Como referido estas situações são do tipo oscilação inversa, ou seja, assiste-se a uma
desaceleração das máquinas pertencentes ao conjunto das máquinas criticas, o que corresponde
a uma diminuição do ângulo rotórico, ao invés de uma aceleração, correspondente a um
aumento do ângulo rotórico, como é habitual nos casos de estabilidade transitória analisados.
55
Uma das possíveis explicações para os resultados errados apresentados pelo método
implementado é a de que os critérios de identificação das máquinas críticas utilizados não têm
em conta a existência deste tipo de instabilidade, estando apenas preparados para os casos em
que se assiste a um aumento do ângulo rotórico da máquina.
Outro facto, muito relevante, verificado através dos resultados de simulação é que nas
perturbações 5 e 6, consoante o tempo de actuação das protecções definido, se assistem a
diferentes tipo de instabilidade, ou seja, tanto se observam situações em que o ângulo rotórico
aumenta, o que corresponde a uma situação de instabilidade da primeira oscilação, como
situações em que diminui, o que corresponde a uma situação de instabilidade de oscilação
inversa.
• Perturbação 22
O método híbrido não consegue convergir quando analisa a perturbação 22, tendo-se
concluído que o problema reside no cálculo da margem de estabilidade transitória. Observando
os dados obtidos na simulação desta contingência constata-se que o problema reside no facto de
a aproximação polinomial não descrever de forma correcta a curva de potência da máquina
equivalente em virtude esta apresentar uma forma muito diferente da habitual semelhante a uma
sinusoidal, como se pode observar na Figura 5.1.
Figura 5.1 - Potência eléctrica gerada (-.), aproximação polinomial da potência eléctrica gerada (linha contínua) e
potência mecânica (--) da máquina equivalente para a perturbação 22
Como se observa a aproximação polinomial é monotonamente crescente, não
intersectando a curva da potência mecânica equivalente, ponto onde se terminava o cálculo da
área de desaceleração para posterior cálculo da margem de estabilidade transitória. Esta forma
pouco habitual tem origem no facto de a evolução no tempo do ângulo rotórico equivalente e da
potência eléctrica equivalente gerada não apresentarem o comportamento sinusoidal que seria
de esperar, como se pode observar na Figura 5.2.
56
Figura 5.2 - Ângulo rotórico e potência eléctrica gerada da máquina equivalente para a contingência 22
Conclui-se então que o programa não está habilitado para lidar com situações como a
que foi agora descrita.
• Perturbações 15 e 24
Como já foi referido, o método implementado conseguiu obter valores correctos para as
contingências 15 e 24, mas estes não foram obtidos de forma autónoma. Quer isto dizer que foi
necessária uma acção externa para obter os resultados correctos, acção esta que se traduziu
pela especificação do conjunto das máquinas críticas ou por uma estimativa inicial mais próxima
do tempo crítico de actuação das protecções a determinar. Caso contrário, verificou-se que ou o
método não convergia ou convergia para valores errados.
5.1.3. Correcções a implementar
As situações erradas analisadas permitem identificar onde residem os problemas do
método implementado: identificação do conjunto das máquinas críticas, aproximação da curva de
potência da máquina equivalente no período pós-defeito e estimativa inicial do tempo crítico de
actuação das protecções. Cada um destes factos será de seguida abordado individualmente.
• Identificação do conjunto das máquinas críticas
Como verificado, os critérios utilizados para identificar o conjunto das máquinas críticas
não oferecem resultados inequívocos.
Esta limitação pode ser corrigida recorrendo a outros critérios de identificação das
máquinas críticas como o Critério das acelerações ou critérios baseados na Função de Energia
Transitória que, conjugados com os critérios já utilizados, vão permitir identificar o conjunto
pretendido.
No código desenvolvido estão já escritos os algoritmos para aplicar os Critérios das
acelerações e baseados na FET, estando em falta o necessário para a sua integração com os
critérios de momento utilizados.
57
• Aproximação da curva de potência da máquina equivalente
O problema referente à aproximação da curva de potência da máquina equivalente
poderá ser contornando utilizando condições de despiste que impeçam o programa de entrar
num ciclo infinito. Estas condições poderão verificar tanto o comportamento da aproximação
obtida como o próprio comportamento da curva de potência da máquina equivalente que, como
foi apresentado, apresenta por vezes formas diferentes das habituais.
• Estimativa inicial do tempo crítico de actuação das protecções
No método implementado, a estimativa inicial do tempo crítico de actuação das
protecções é definida pelo utilizador, cabendo depois ao programa convergir para o valor final.
Contudo, como se constatou, para algumas das situações analisadas uma estimativa inicial
afastada do valor final faz com que o método não convirja.
Para resolver este problema, pode ser determinada uma estimativa inicial do tempo
crítico, no início do programa, recorrendo ao critério das áreas iguais generalizado, apresentado
no Capítulo 3.
5.2. Exemplo de aplicação
Apresenta-se de seguida uma exemplo de aplicação do método híbrido implementado.
Como foi já explicado, o método efectua vários cálculos pré-defeito como, por exemplo, um
trânsito de energia (resultado no Anexo 1) e o cálculo dos valores iniciais das forças
electromotrizes e também algumas acções finais, como o armazenamento dos dados. No
entanto, este exemplo está centrado no ciclo responsável por determinar o tempo crítico, dado
que tanto as acções preliminares como as finais são já conhecidas, não contribuindo com
nenhum facto novo.
A situação analisada foi um curto circuito trifásico simétrico, na rede de teste da CIGRE,
junto ao barramento 1, que foi eliminado pela retirada de serviço da linha 1, o que corresponde à
perturbação 2 da Tabela 5.1. O ficheiro de entrada “contingencias.xlsx” tem o conteúdo da
Tabela 5.2.
Tabela 5.2 - Conteúdo do ficheiro de entrada utilizado no exemplo de aplicação
2 1 1
Após a leitura de todos os dados dos ficheiros de entrada e calculados todos os valores
pré-defeito o ciclo para determinar o tempo crítico é iniciado com uma estimativa inicial do tempo
crítico de 𝑡𝑡𝑐𝑟0 = 0,6 s. Para esta primeira iteração obtêm-se as evoluções temporais dos ângulos
rotóricos da Figura 5.3, onde se podem observar os ângulos rotóricos reais (em cima) e os
referidos ao centro de inércia do sistema (em baixo). Apesar de não ser visível na figura
considerou-se que o defeito acontece 0,1 s depois do início da simulação, ou seja, 𝑡𝑡𝑑𝑒𝑓 = 0,1 s.
Em qualquer dos casos é possível ver claramente que a máquina 1 (linha a tracejado) constitui o
58
conjunto das máquinas críticas, facto que é corroborado pelos algoritmos de determinação do
conjunto das máquinas críticas, como se verá posteriormente. O processo de integração
numérica é interrompido mais cedo devido à existência dos índices de estabilidade e
instabilidade mas, neste caso, optou-se por deixar correr a simulação durante o tempo total para
melhor se observarem as curvas.
Figura 5.3 - Ângulos rotóricos para a perturbação 2, com 𝒕𝒄𝒄𝒍 = 𝟎,𝟔 s
Para além disso, observa-se que o tempo de actuação das protecções especificado
conduz a uma situação instável, o que está de acordo com os dados obtidos através dos índices
de estabilidade e instabilidade. Como se observa na Figura 5.4, o índice IDCS atinge primeiro um
mínimo, o que corresponde a uma situação instável e permite parar o processo de integração
numérica, sendo apenas necessário efectuar a integração numérica das equações que
descrevem o sistema num período de aproximadamente 0,7 s.
Figura 5.4 - Índices IDCS e IDE para a perturbação 2, com 𝒕𝒄𝒄𝒍 = 𝟎,𝟔 s
Como a situação é considerada instável, o critério utilizado para determinar o conjunto
das máquinas críticas é o CMR. Da aplicação deste critério obtêm-se os resultados da Tabela
5.3, de onde se constata que apenas a máquina 1 pertence ao conjunto das máquinas críticas.
restantes máquinas
restantes máquinas
máquina 1
máquina 1
59
Tabela 5.3 - Identificação do conjunto das máquinas críticas para a perturbação 2, com 𝒕𝒄𝒄𝒍 = 𝟎,𝟔 s
Máquina Ângulos rotóricos 𝛿𝛿𝑖 [rad] Intervalo [rad] 1 7,012 6,362 3 0,700 2,855 7 0,699 2,046 5 0,666 1,581 2 0,628 1,367 6 0,618 1,121 4 0,600 0,000
Após a identificação do conjunto das máquinas críticas, o sistema é reduzido a uma
máquina equivalente, cujas curvas temporais do ângulo rotórico e potência eléctrica gerada se
encontram na Figura 5.5.
Figura 5.5 - Curvas temporais do ângulo rotórico e potência eléctrica gerada da máquina equivalente, para a
perturbação 2, com 𝒕𝒄𝒄𝒍 = 𝟎,𝟔 s
O cálculo da margem de estabilidade transitória, é efectuado utilizando a aproximação
sinusoidal, como se observa na Figura 5.6. Repare-se que devido às expressões utilizadas a
máquina equivalente apresenta, no período pós-defeito, valores de potência eléctrica gerada
negativos.
Figura 5.6 - Curva de potência da maquina equivalente (a tracejado) e aproximação sinusoidal (a cheio), para a
perturbação 2, com 𝒕𝒄𝒄𝒍 = 𝟎,𝟔 s
60
Obtém-se o valor de 𝜂𝜂0 = −6,7809, cuja polaridade está de acordo com o facto de a
situação analisada ser instável.
Dado que só se dispõe de um valor da margem de estabilidade transitória, não é possível
fazer uma nova estimativa do tempo crítico por um processo de regressão linear. O valor
necessário para a próxima iteração é calculado por 𝑡𝑡𝑐𝑟1 = 0,9 × 𝑡𝑡𝑐𝑟0 = 0,54 s, sendo a constante
0,9 definida por via experimental.
A segunda iteração é então iniciada com uma estimativa para o valor do tempo crítico de
𝑡𝑡𝑐𝑟1 = 0,54 s. Como os passos até à obtenção do valor da margem de estabilidade em cada
iteração são semelhantes, a sua repetição é desnecessária, apresentando-se já o valor da
margem de estabilidade transitória. Na segunda iteração obtém o valor de 𝜂𝜂1 = −5.2493. Usando
os valores já disponíveis da margem de estabilidade transitória e tempos de actuação das
protecções, faz-se então uma nova estimativa do tempo crítico, obtendo-se 𝑡𝑡𝑐𝑟2 = 0.334 s, o que
corresponde a uma variação em relação ao valor anterior de ∆𝑡𝑡𝑐𝑟 = 0,206 s. Este processo
repete-se até que a variação entre duas estimativas consecutivas seja menor que uma tolerância
𝛼𝛼. 𝜀𝜀 = 10 × 0,001 s = 0,01 s. Os dados relativos a estas iterações estão na Tabela 5.4.
Tabela 5.4 - Dados relativos às iterações realizadas com a aproximação sidusoidal
Iteração Tempo crítico
estimado 𝑡𝑡𝑐𝑟𝑖
Margem de estabilidade
transitória 𝜂𝜂𝑖 Variação ∆𝑡𝑡𝑐𝑟
- 0,600 -6,781 - 1 0,540 -5,249 0,060 2 0,334 9,293 0,206 3 0,466 -1,252 0,131 4 0,450 0,160 0,016 5 0,452 - 0,002
Como se observa, na 5ª iteração a variação entre dois valores consecutivos das
estimativas do tempo crítico é menor que a tolerância 𝛼𝛼. 𝜀𝜀 definida anteriormente o que significa
que, daqui em diante, no cálculo da margem de estabilidade transitória é utilizada a aproximação
polinomial de segundo grau. Para a primeira iteração com a aproximação polinomial, apenas se
dispõe de um valor calculado da margem de estabilidade transitória, não sendo por isso possível
estimar um novo valor de tempo crítico usando uma regressão linear. Devido a este facto, o
próximo valor é calculado por 𝑡𝑡𝑐𝑟6 = 0,98 × 𝑡𝑡𝑐𝑟5 = 0,443 s, sendo a constante 0,98 também definida
por via experimental. Os valores relativos a estas iterações estão disponíveis na Tabela 5.5.
Tabela 5.5 - Dados relativos às iterações realizadas com a aproximação polinomial
Iteração Tempo crítico
estimado 𝑡𝑡𝑐𝑟𝑖
Margem de estabilidade
transitória 𝜂𝜂𝑖 Variação ∆𝑡𝑡𝑐𝑟
5 0,452 -0,606 - 6 0,443 0,149 0,009 7 0,445 0,010 0,002 8 0,445 - 0,000
61
Observa-se que o valor que seria utilizado para a oitava iteração já apresenta uma
variação menor que a tolerância 𝜀𝜀, definida como 𝜀𝜀 = 0,001 s, em relação ao valor de tempo
crítico anterior. Por esse facto não é necessário realizar a 8ª iteração e o valor de tempo crítico
estimado, em 7 iterações, para a perturbação 2 é, descontando o tempo de defeito 𝑡𝑡𝑑𝑒𝑓 referido
anteriormente, igual a 𝑡𝑡𝑐𝑟 = 0,445 − 𝑡𝑡𝑑𝑒𝑓 = 0,345 s.
5.3. Conclusões
Neste capítulo apresentaram-se os resultados dos testes computacionais do método
implementado, a partir dos quais se retiraram conclusões acerca das limitações do método.
Apresentou-se também um exemplo detalhado de aplicação do método híbrido implementado.
Relativamente aos resultados dos testes computacionais observou-se que para 6 das 29
situações testadas, não se obtiveram resultados correctos. Verificou-se que estas falhas se ficam
a dever a três secções do método: identificação do conjunto das máquinas críticas, cálculo da
margem de estabilidade transitória e estimativa inicial do valor do tempo crítico. Para cada uma
delas foram sugeridas medidas que podem ajudar a resolver os problemas detectados.
Quanto ao exemplo detalhado, foi dado especial destaque ao ciclo para determinar o
tempo crítico de actuação das protecções. Foi possível observar, para uma determinada
perturbação, a evolução quer das estimativas do tempo crítico, quer dos valores calculados da
margem de estabilidade transitória. Observou-se também o instante para o qual se altera o
método utilizado para aproximar a curva de potência da máquina equivalente.
63
6. Conclusões e Propostas para trabalhos futuros
Neste capítulo é apresentada uma síntese do trabalho realizado, percorrendo todos os
capítulos até agora apresentados. Inclui-se também uma perspectiva sobre trabalhos futuros que
tenham como ponto de partida o trabalho aqui apresentado.
6.1. Conclusões
Nesta dissertação desenvolveu-se um método híbrido para o estudo da estabilidade
transitória do SEE, em linguagem matlab. Este método, que junta as vantagens dos Métodos
Directos, os quais permitem uma análise rápida, com as vantagens dos Métodos de Integração
Numérica, os quais permitem uma grande capacidade de modelização do sistema, é capaz de,
para uma determinada rede, calcular o tempo crítico de actuação das protecções para diferentes
perturbações especificadas pelo utilizador.
O método caracteriza-se por utilizar um modelo equivalente do sistema, constituído por
uma máquina ligada a um barramento de potência de curto-circuito infinita. Este modelo
equivalente é calculado após a determinação do conjunto de máquinas críticas, o qual é obtido
usando índices que, com base nos valores obtidos por integração numérica que descrevem o
funcionamento temporal das máquinas síncronas, identificam as máquinas responsáveis pela
perda de sincronismo do sistema. O sistema equivalente obtido é depois estudado utilizando o
Método das Áreas Iguais, o qual permite obter a margem de estabilidade transitória para a
perturbação em estudo, a partir da qual se vai determinar o tempo crítico de actuação das
protecções para esta perturbação.
Como se verificou este objectivo foi parcialmente atingido pois o método desenvolvido
mostrou oferecer resultados para a maior parte dos casos testados, revelando no entanto
problemas na análise de certas situações.
Faz-se agora uma breve passagem sobre os diferentes capítulos apresentados.
No Capítulo 1, a Introdução, fez-se uma contextualização sobre o tema desta tese, onde
se caracterizou o SEE actual, se apresentaram conceitos relativos ao tema e se mostrou a
necessidade de existirem métodos que permitam uma rápida análise da estabilidade transitória
dos diferentes pontos de funcionamento do sistema.
No Capítulo 2, designado Modelização do Sistema de Energia Eléctrica, foram
apresentados os modelos utilizados, nos métodos dedicados ao estudo da estabilidade
transitória, para descrever os diferentes componentes da rede e ainda as contingências
simuladas. Foi utilizado o modelo clássico para descrever os diferentes componentes do SEE e
as contingências simuladas são apenas curto-circuitos trifásicos simétricos junto aos
barramentos.
64
No Capítulo 3 apresentaram-se os diferentes métodos utilizados no estudo da
estabilidade transitória do SEE. Como se verificou, estes podem dividir-se em quatros grupos:
Métodos de Integração Numérica, Métodos Directos, Métodos Híbridos e Técnicas de
Inteligência Artificial. Este capítulo foi designado Métodos de análise da estabilidade transitória
do SEE. Os Métodos de Integração Numérica caracterizam-se por oferecer uma modelização
detalhada mas elevado esforço computacional enquanto que os Métodos Directos, por outro
lado, oferecem uma modelização pouco detalhada mas são muito eficientes do ponto de vista
computacional. Os Métodos Híbridos conjugam as capacidades de modelização oferecidas pelos
Métodos de Integração Numérica com a eficiência computacional oferecida pelos Métodos
Directos. As Técnicas de Inteligência Artificial são uma abordagem mais recente ao estudo da
estabilidade transitória do SEE e caracterizam-se por exigirem um esforço computacional muito
elevado na fase de implementação mas posteriormente possibilitarem uma análise rápida em
tempo real.
O Capítulo 4, designado Método híbrido implementado, foi dedicado à apresentação do
método desenvolvido. Procurou-se explicar a estrutura geral do método, bem como as suas
particularidades consideradas mais importantes. Mostrou-se que o método está divido em três
etapas: acções preliminares, ciclo para determinação do tempo crítico, acções finais. Como foi
apresentado, o método caracteriza-se por interromper o processo de integração numérica antes
do tempo total de simulação e por utilizar um modelo equivalente do sistema, constituído por uma
máquina equivalente ligada a um barramento de potência de curto-circuito infinita, o qual é
estudado recorrendo ao Método das Áreas Iguais. Este fornece o valor da margem de
estabilidade transitória e, com base em regressões lineares que utilizam os tempos de actuação
das protecções e respectivos valores da margem de estabilidade transitória, são estimados os
valores de tempo crítico para as perturbações definidas pelo utilizador. Terminou-se fazendo a
correspondência entre as diferentes etapas do método e os ficheiros de código.
O Capítulo 5, denominado Resultados computacionais, dividiu-se em duas partes
distintas. A primeira, os Testes computacionais, na qual se pôs à prova o método desenvolvido,
procurando identificar os erros obtidos nos resultados e quais as suas causas e possíveis
resoluções. Como se verificou, o método apresenta ainda, em situações pontuais, resultados
menos correctos, sendo as principais causas detectadas o facto de o tipo de instabilidade ser de
oscilação inversa, a identificação do conjunto das máquinas críticas não ser inequívoca e a
aproximação à curva de potência da máquina equivalente não se revelar adequada. Uma
segunda parte, Exemplo de aplicação, na qual se mostrou detalhadamente o funcionamento do
método na análise de uma determinada perturbação, dando especial ênfase ao ciclo para
determinação do tempo crítico de actuação das protecções.
Por fim, pode-se afirmar que, apesar de o método desenvolvido não se ter mostrado apto
para analisar todo o tipo de perturbações para as quais foi testado, o trabalho desenvolvido pode
ser um ponto de partida para futuros trabalhos. Esta afirmação baseia-se no facto de, apesar de
o método desenvolvido não estar completamente funcional, terem sido apresentados conceitos
65
fundamentais para o desenvolvimento de métodos deste tipo. Entre estes conceitos podem-se
referir por exemplo os índices usados para identificar o conjunto de máquinas críticas ou o tempo
de paragem da simulação.
6.2. Propostas para trabalhos futuros
Como já foi referido, o trabalho desenvolvido pode ser um bom ponto de partida para
outros trabalhos, essencialmente pelos conceitos aqui apresentados. Devido ao facto de o
método não ser completamente fiável, uma linha de desenvolvimento a seguir pode ser essa.
Numa fase posterior poderão ser incorporadas novas funcionalidades.
Relativamente ao facto de o método não ser completamente fiável, foi mostrado no
Capítulo 5 quais os erros que o método apresentava. Procurou-se também identificar quais as
suas causas sendo de seguida apresentadas várias hipóteses que poderiam solucionar os erros
detectados. Trabalhar sobre essas hipóteses é portanto uma linha a seguir por possíveis
trabalhos que tenham este como ponto de partida.
Quanto à incorporação de novas funcionalidades, pode-se considerar que esta linha de
desenvolvimento apresenta um grande conjunto de hipóteses. Primeiro, porque o programa
apenas apresenta um tempo crítico de actuação das protecções. Em trabalhos futuros o
programa pode também servir para calcular limites de potência ou para classificar e ordenar
contingências. Segundo, porque as perturbações consideradas foram curto-circuitos trifásicos
simétricos que ocorrem junto ao barramento. Trabalhos futuros podem dar a possibilidade de
simular defeitos não simétricos ou que não ocorram na extremidade das linhas. Terceiro porque
foi utilizado o modelo clássico para caracterizar os componentes do sistema. Trabalhos futuros
podem oferecer uma modelização mais detalhada dos componentes do SEE.
67
Bibliografia
1. Kundur, P. Power System Stability and Control. USA : McGraw-Hill, Inc., 1994.
2. Sucena Paiva, J. P. Redes de Energia Eléctrica: uma análise sistémica. 2ª edição.
Lisboa : IST Press, 2007.
3. Machado Ferreira, C. M. Análise da Estabilidade Transitória de Sistemas Eléctricos
de Energia utilizando Formulações Híbridas. Porto : Departamento de Engenharia Electrotécnica
e de Computadores, Faculdade de Engenharia da Universidade do Porto, Dissertação de
Doutoramento, 2005.
4. Holmes, M. H. Introduction to Numerical Methods in Differential Equations. New York :
Springer Science+Business Media, LLC, 2007.
5. Vittal, Vijay. Direct Stability Methods. [autor do livro] Leonard L. Grigsby. Electric
Power Engineering Handbook. Second Edition. Power System Stability and Control. Boca Raton,
USA : CRC Press - Taylor & Francis Group, 2007.
6. Machowski, Jan, Bialek, Janusz W. e Bumby, James R. Power System Dynamics:
Stability and Control. Second Edition. Chichester, UK : John Wiley & Sons, Ltd., 2008.
7. Sauer, Peter W. e Pai, M. A. Power System Dynamics and Stability. New Jersey,
USA : Prentice-Hall, Inc., 1998.
8. Xue, Y. e al., et. Extended equal area criterion revisited. IEEE Transactions on Power
Systems. 7, No. 3, August 1992.
9. Machado Ferreira, C. M., Dias Pinto, J. A. e Maciel Barbosa, F. P. Transient
Stability Assessment on an Electric Power System using Trajectory Sensitivity Analysis. 39th
International Universities Power Engineering Conference. September 2004.
10. Chan, K. W., Zhou, Q. e Chung, T. S. Transient Stability Margin Assessment for
Large Power System using Time Domain Simulation Based Hybrid Equal Area Criterion Method.
Proceedings of the 5th International Conference on Advances in Power System Control,
Operation and Management. October 2000.
11. MATLAB Documentation. [Online] http://www.mathworks.com/help/techdoc/.
69
Anexo 1
Neste anexo é apresentada a rede utilizada para avaliar o desempenho do método
desenvolvido. A rede é a rede de teste da CIGRE, constituída por 7 geradores, 10 barramentos e
13 linhas, como se observa na Figura A1.1. A potência base são 100 MVA e a frequência
nominal 60 Hz.
Figura A1.1 - Rede de teste da CIGRE
De seguida serão apresentadas as características dos diferentes componentes da rede
(geradores e linhas), seguidos das características dos barramentos e o resultado do trânsito de
potências.
Na Tabela A1.1 estão disponíveis as características das linhas que constituem a rede.
Tabela A1.1 - Características das linhas da rede de teste da CIGRE
Nº da linha
Barramento 1
Barramento 2
Resistência longitudinal [p.u.]
Reactância longitudinal [p.u.]
Susceptância transversal [p.u.]
1 1 3 0,00988 0,0484 0,2025 2 1 4 0,00988 0,0484 0,10126 3 2 3 0,04504 0,12365 0,2025 4 2 10 0,0164 0,0638 0,30376 5 3 4 0,01185 0,07803 0,30376 6 3 9 0,01146 0,05531 0,30376 7 4 5 0,00395 0,01975 0,10126 8 4 6 0,00751 0,01975 0,6075 9 4 9 0,0488 0,19161 0,30376 10 4 10 0,0164 0,06519 0,405 11 6 8 0,01877 0,06282 0,2025 12 7 8 0,01185 0,07803 0,30376 13 8 9 0,0488 0,19161 0,2025
70
Na Tabela A1.2 podem observar-se as características referentes aos 7 geradores.
Tabela A1.2 - Características dos geradores da rede de teste da CIGRE
Nº do gerador
Barramento de ligação
Constante de inércia 𝐻 [s]
Reactância transitória 𝑋𝑋𝑑′ [p.u.]
1 1 11,35 0,074 2 2 7,75 0,118 3 3 14,31 0,062 4 4 17,98 0,049 5 5 11,35 0,074 6 6 12,76 0,071 7 7 10,71 0,087
Na Tabela A1.3 mostram-se os dados referentes aos barramentos e os resultados da
resolução do trânsito de energia.
Tabela A1.3 - Resultados do trânsito de energia da rede de teste da CIGRE
Nº do barramento
Tipo do barramento
Módulo da
tensão [p.u.]
Argumento da tensão
[graus]
Potência activa gerada [MVA]
Potência reactiva gerada [MVA]
Potência activa de
carga [MVA]
Potência reactiva de carga
[MVA] 1 Ref. 1,052 0,000 0,00 68,08 0,00 0,00 2 PV 1,006 -7,166 120,00 117,21 200,00 120,00 3 PV 1,037 -1,523 256,00 58,01 0,00 0,00 4 PV 1,01 -3,642 300,00 151,72 650,00 405,00 5 PV 1,037 -1,372 230,00 94,48 0,00 0,00 6 PV 1,015 -3,235 160,00 19,45 80,00 30,00 7 PV 1,011 -0,746 174,00 35,02 90,00 40,00 8 PQ 1 -4,411 0,00 0,00 100,00 50,00 9 PQ 1 -6,641 0,00 0,00 230,00 140,00 10 PQ 1 -7,029 0,00 0,00 90,00 45,00
