Técnicas de solução de sistemas de equações diferenciais e algébricas: aplicação em sistemas de energia elétrica

SISTEMAS DE POTÊNCIA

Técnicas de solução de sistemas de equações diferenciais e algébricas: aplicação em sistemas de energia elétrica

José E. O. Pessanha^I; Carlos Portugal^I; Alex A. Paz^II

^IUFMA-CPGEE-DEE-UFMA, Avenida dos Portugueses s/n, Campus do Bacanga, São Luís, Ma, - Brasil -65080-040, [e-mail], [e-mail]

^IIPUC/RJ-DEE, Rua Marquês de São Vicente, 225, Gávea, Rio de Janeiro, RJ - Brasil - 22453-900, [e-mail]

RESUMO

O presente trabalho investiga e compara o desempenho de técnicas numéricas aplicadas na solução de Equações Diferenciais e Algébricas (EDAs) representando sistemas elétricos e seus componentes. Entre os métodos considerados, destaca-se o método conhecido como Método de Diferenciação Regressiva Modificado Estendido (MEBDF), por este apresentar propriedades de estabilidade numérica que os Métodos de Diferenciação Regressiva (BDF) convencionais não apresentam, sendo que estas propriedades podem evitar problemas numéricos durante a solução das EDAs. As técnicas são comparadas através de testes envolvendo simulações computacionais no domínio do tempo de fenômenos de estabilidade de curta-, e de longa-duração (angular e de tensão, respectivamente) em sistemas-teste do IEEE. O objetivo principal foi verificar a eficiência dessas técnicas numéricas sob o ponto de vista computacional, sem comprometer a precisão. O aspecto computacional está relacionado com o tempo de cpu gasto nas simulações.

Palavras-chave: Equações Diferenciais e Algébricas, Análise Numérica, Estabilidade de Sistemas Elétricos, Simulação de Sistemas Elétricos, Análise no Domínio do Tempo.

ABSTRACT

The present work investigates and compares the computational performance of numerical techniques applied to the solution of algebraic and differential equations (DAEs) representing power systems and their components. Among the considered techniques, emphasis was given to the method known as Modified Extended Backward Differentiation Formulae (MEBDF), due to its numerical stability properties that conventional Backward Differentiation Formulae (BDF) methods do not have, and these properties may avoid numerical troubles during DAEs solution. The comparison is made through computational tests related to simulations of power system short- and long-term stability phenomena (transient angular stability and long-term voltage stability), using IEEE system-test. The main objective was to check the computational efficiency of these numerical techniques. The computational aspect is related to the simulation CPU time.

Keywords: Differential algebraic equations, Numerical analysis, Power system stability, Power system simulation, Time domain analysis.

1 INTRODUÇÃO

Simulação no domínio do tempo é uma forma de análise muito útil para as concessionárias de energia elétrica desenvolverem estudos computacionais envolvendo fenômenos de estabilidade em sistemas de energia elétrica. Um programa computacional para análise de estabilidade capaz de simular fenômenos de curta- e de longa-duração requer a solução de sistemas de Equações Diferenciais e Algébricas (EDAs) não lineares, rígidas e de grande porte (Astic et al., 1994), processo que pode exigir elevado tempo computacional. Este esforço depende principalmente das características dos métodos numéricos usados, da complexidade dos modelos matemáticos implementados, da dimensão do sistema elétrico simulado, das constantes de tempo envolvidas, da velocidade do fenômeno simulado (curta- ou longa-duração), da capacidade do computador e do tempo total da simulação (Paz, 2004). As análises e simulações podem se tornar ainda mais complexas se diferentes formas de estabilidade se manifestarem simultaneamente.

O uso de modelos simplificados pode resultar numa redução no tempo de simulação computacional, entretanto, dependendo do cenário simulado, pode reduzir também a precisão dos resultados. Por outro lado, modelos mais complexos irão exigir um considerável esforço computacional na solução dos sistemas de equações diferenciais e algébricas que descrevem o comportamento dinâmico do sistema e de seus componentes, sendo que esta complexidade pode ser desnecessária para certos cenários (Paz, 2004).

A implementação num mesmo programa computacional de modelos matemáticos representando dispositivos de controle de resposta rápida, e também de resposta lenta, permite quando necessário, a "captura" de efeitos inerentes a fenômenos de curta- e de longa-duração. Para simular fenômenos rápidos e lentos num mesmo cenário, além da necessidade de disponibilizar modelos com essas características, é relevante que o método numérico usado na resolução dos sistemas de equações diferenciais e algébricas seja eficiente. A eficiência de um método numérico neste caso pode ser avaliada em termos de precisão/confiabilidade de resultados e de tempo de processamento gasto (CPU). Uma boa precisão é garantida, ou não, pela capacidade em manter sob controle os erros associados às técnicas de aproximação usadas pelo método. Já a sua eficiência computacional é medida em termos de tempo de CPU, quanto menor, melhor sua eficiência computacional.

Nota-se um crescente interesse pelo desenvolvimento/adaptação de metodologias numéricas apilcadas na solução de sistemas de EDAs (Ascher e Petzold, 1998; Brenan et al., 1996; Cash, 2000; Mazzia e Iavernaro, 2003; Petzold e LI, 2000). Este interesse tem produzido e disponibilizado novas alternativas de solução e o presente trabalho procura integrar e aproveitar adequadamente estes avanços com interesse específico em fenômenos de estabilidade em sistemas elétricos.

As técnicas numéricas abordadas neste trabalho são: Método de Diferenciação Regressiva Modificado Estendido Implícito (Cash, 1983), Método de Diferenciação Regressiva Modificado Estendido com Esparsidade (Abdulla et al., 2000), Método de Diferenciação Regressiva (Petzold, 1983) e Runge-Kutta 4ª ordem. Na literatura esses métodos são identificados como MEBDF, MEBDFSD, BDF e RK4. Há um interesse particular nas diferentes estratégias usadas na avaliação do passo, da ordem, no controle de erro e propriedades de estabilidade numérica destacadando-se as vantagens, desvantagens e diferenças em cada um dos estágios do processo de integração.

2 EQUAÇÕES DIFERENCIAIS E ALGÉBRICAS

Um sistema de EDAs corresponde a um tipo de Equações Diferenciais Ordinárias (EDOs) cuja solução deve satisfazer restrições impostas por equações algébricas não-lineares. Pode-se então considerar as EDAs como uma extensão do sistema de EDOs com restrições, representadas na forma semi-explícita pelas Equações (1) e (2), onde y, y', e f são vetores de dimensão n, com y representando as n variáveis diferenciais, z as m variáveis algébricas e p os k parâmetros do sistema (no caso dos sistemas de energia, esses parâmetros definem uma configuração específica e condição de operação).

O interesse por modelos matemáticos complexos baseados em sistemas de EDAs de grande porte com índices superiores (maiores que um) tem crescido bastante. Deve-se ter cuidado e atenção na formulação das EDAs representando modelos muito complexos, uma vez que sistemas de índices superiores apresentam uma grande tendência à instabilidade (Ascher e Petzold, 1998).

Um sistema de EDAs pode ser convertido à um sistema de EDOs, sendo que este procedimento pode ser muito complexo. É possível reduzir o índice de um sistema de EDAs através de sucessivas diferenciações. Um sistema de índice diferencial igual a zero é considerado como um sistema de equações diferenciais ordinárias e a partir deste ponto aplica-se métodos de solução bem conhecidos (Hairer e Wanner, 1996). Embora existam metodologias e estratégias eficientes de solução das EDOs, existem situações onde é preferível resolver diretamente as EDAs sem convertê-las inicialmente a um sistema de EDOs. Uma delas está relacionada com o uso de uma metodologia que permite a criação e incorporação de modelos pré-definidos para gerar automaticamente o sistema de EDAs e iniciar o processo de solução. Mantendo o sistema de EDAs inalterado, são mantidas as características esparsas do sistema e o significado físico das variáveis, evitando a perda de informações importantes (Secanell e Córcoles, 2002).

2.1 Formas de Representação

Embora a forma geral implícita do sistema de EDAs possa representar uma variedade de condições, existem casos especiais onde esse tipo de representação pode resultar em dificuldades no processo de solução se não forem identificados seus índices (medida da singularidade do sistema). Um sistema com índices maiores que 1 geralmente apresenta dificuldades durante o processo de solução. Felizmente, a maioria dos problemas encontrados na prática podem ser expressos como uma combinação de estruturas de EDOs acopladas a equações de restrição, que podem ser representados pela Equação (3).

Uma característica interessante das EDAs com índices superiores é que apresentam restrições ocultas, restrições estas que podem ser encontradas por diferenciação e manipulação algébrica no processo de estimação do índice (Ascher e Petzold, 1998; Fábián et al., 2000; Tarraf e Asada, 2002).

2.2 Rigidez das EDAs

A rigidez está relacionada com grandes diferenças nas constantes de tempo de um sistema de EDAs, ou seja, com a taxa de decaimento da solução. A definição matemática correta sobre rigidez de um sistema de EDAs ainda gera polêmica entre os especialistas. A opinião mais pragmática e também a mais antiga (Curtiss & Hirschfelder 1952) é: equações rígidas são equações onde certos métodos implícitos, em particular o BDF, apresentam melhor desempenho, geralmente muito melhor em comparação aos explícitos. (Hairer e Wanner, 1996).

Um sistema de EDAs possui diferentes taxas de decaimento por cada variável de estado, e estão relacionadas com os autovalores da matriz Jacobiana (¶f/¶y). Esses autovalores desempenham uma função importante nesta decisão, mas quantidades como a dimensão do sistema, a suavidade da solução ou o intervalo de integração também são importantes. As variáveis com taxa de decaimento lenta são responsáveis pela determinação do erro de truncamento, mas existem variáveis com taxa de decaimento muito rápida, que diminuem até atingir valores muito pequenos (Figura 1). Estas são responsáveis pela determinação de comprimentos de passo muito pequenos que tentam assegurar a estabilidade numérica da metodologia de integração, tornando o processo muito demorado em termos de tempo computacional. Este problema é conhecido como problema rígido (Gear, 1971). Portanto, pode-se resumir que um sistema de EDAs é considerado rígido se um método numérico é obrigado a usar passo de integração muito pequeno em relação à suavidade da solução exata do problema no intervalo em questão (Jardim, 1997).

É recomendável o uso de metodologias implícitas que possuam melhores propriedades de estabilidade sem muitas restrições na escolha do comprimento do passo. Atualmente, as metodologias de passo variável BDF e MEBDF permitem alcançar ordens de convergência menores que 2 e 4, respectivamente, mantendo sua característica A-estável. Portanto, estas são consideradas na literatura como as técnicas mais adequadas para a solução numérica de problemas de valor inicial de sistemas rígidos de EDAs (Ascher e Petzold, 1998; Hairer e Wanner, 1996).

3 MÉTODOS DE SOLUÇÃO

3.1 Método de Diferenciação Regressiva

Métodos BDF multipasso são muito utilizados para resolver sistemas de EDAs rígidas (Hairer e Wanner, 1996). As fórmulas BDF são construídas de acordo com um processo de interpolação dos pontos solução (y_{j – k + 1},..., y_j e y_{j + 1}) através de em interpolação polinomial e baseado em combinação linear dos polinômios de Lagrange. A Fórmula geral de Diferenciação Regressiva está representada pela Equação (4) onde o valor do coeficiente a_r é calculado pela Equação (5).

A metodologia baseada nas equações anteriores foi desenvolvida para resolver problemas de valor inicial na seguinte forma implícita (Petzold, 1983):

F, y e y' são vetores de dimensão n. A derivada presente na Equação (6) é substituída por diferenças regressivas, processo conhecido como discretização matemática, e o sistema de equações não-lineares resultante é solucionado por uma variante do Método de Newton. Por exemplo, substituindo a derivada da Equação (6) pela primeira diferença, a seguinte fórmula implícita de Euler é obtida:

A expressão Dt_{j + 1} = t_{j + 1} – t_j representa o comprimento de passo no intervalo t_{j + 1}. Ao invés de usar uma fórmula de primeira ordem, a técnica aproxima a derivada usando uma fórmula da família BDF de ordem de convergência q = k. Em cada intervalo de integração é escolhida a ordem k e o comprimento de passo Dt_{j + 1}, baseado no comportamento da solução.

Resultado do processo de discretização, a Equação (8) representa o sistema de equações não-lineares que deve ser resolvido a cada intervalo de tempo, onde a é uma constante que varia quando o comprimento do passo, ou a ordem, variam, e b é um vetor que depende da solução nos intervalos de tempo anteriores. Portanto t, y, a e b são avaliados no instante t_{j + 1}.

O método BDF utiliza uma fórmula explícita para o estágio previsor e uma fórmula implícita BDF no corretor. O estágio previsor utiliza uma fórmula de extrapolação explícita (ou polinômio de diferenças divididas) que interpola os pontos solução y_{j + 1 – k} nos últimos k intervalos de tempo. O previsor é responsável pela primeira aproximação y(0)_{j + 1} através de uma simples avaliação no intervalo de tempo t_{j +1} dada pela Equação (9).

Onde as diferenças divididas são definidas por:

Esta primeira aproximação é usada no estágio corretor para a determinação da solução final no intervalo t_{j + 1}. Da mesma forma, o vetor y'(0)_{j + 1} é obtido diferenciando-se o polinômio previsor (Equação 9) no instante t_{j + 1}.

Na etapa de correção, algumas hipóteses apresentadas por (Byrne e Hindmarch, 1975; Jackson e Sacks-Davis, 1980) permitem o conhecimento implícito do vetor solução y_{j + 1} no instante t_{j + 1} através da relação com os valores aproximados na etapa previsão. O polinômio corretor interpola o polinômio previsor em k pontos igualmente espaçados anteriores a t_{j + 1} (onde k é a ordem das fórmulas BDF) obtendo-se o sistema de equações algébricas não-lineares representadas por:

Aqui, o estágio corretor é o mais relevante dentro do esquema previsor-corretor, com maior interesse na escolha e implementação do método de solução do sistema de equações não-lineares, onde o parâmetro a é função do comprimento do passo de integração determinado em t_{j + 1} e permanece fixo enquanto não houver variação no comprimento do passo e/ou na ordem do método. Da mesma forma, b permanece fixo durante todo o processo iterativo, uma vez que as funções y(0)_{j + 1} e y'(0)_{j + 1} são calculadas pelo polinômio previsor de acordo com a Equação (9).

Seleção da Ordem e do Passo de Integração

Neste trabalho, a metodologia BDF considerada foi desenvolvida por (Brenan et al., 1996) e utiliza uma estratégia conhecida como Coeficiente Fixo Direcionado (CFD) para decidir qual o comprimento de passo e a ordem para avançar a solução de um intervalo para o próximo. O CFD é uma estratégia híbrida, incluindo as melhores características das estratégias de coeficiente fixo e de coeficiente variável, com a finalidade de melhorar a eficiência do processo de variação da ordem do método e do comprimento de passo de integração, permitindo adaptar métodos BDF de passo fixo para métodos de passo variável. O processo de seleção de ordem e passo de integração é feito através de estratégias baseadas no cálculo dos erros, fixando-se o comprimento dos passos de integração dos últimos intervalos...,t_{j – 1}, t_j e ordens ..., k_{j – 1} , k_j.

Após satisfazer a condição de convergência, inicia-se outro processo de avaliação correspondente a aceitação do comprimento de passo Dt_{j + 1} = t_{j + 1} – t_j através de uma estratégia de cálculo de erro. Esta etapa corresponde a um dos aspectos mais relevantes da metodologia.

Os dois tipos de erros envolvidos no processo de aceitação ou rejeição do passo estão relacionados ao erro de truncamento local do método e ao processo de interpolação polinomial, sendo este último composto de duas partes (Brenan et al., 1996). Se considerados estes erros, é elaborada a seguinte condição de decisão:

O comprimento do passo é rejeitado se a condição acima não for satisfeita e o cálculo da ordem do método independe da aceitação do passo de integração.

Após o teste de verificação do comprimento do passo e do cálculo da ordem do método, se aceitos, segue-se calculando novos passos, sendo realizada uma nova tentativa se fracassar na determinação do comprimento do passo que satisfaça as condições estabelecidas (Shampine e Gordon, 1975). O erro na nova ordem k é estimado como se os últimos j + 1 intervalos fossem calculados com o novo passo. É escolhido o novo passo de integração de forma conservativa dado pela Equação (13) para que o erro seja a metade do valor da tolerância relativa ao erro de integração desejado.

Est representa o erro estimado em função da ordem do método.

Quando o passo é aceito, o seu comprimento aumenta apenas se for possível dobrá-lo. Se for necessária uma redução, este é reduzido de um fator mínimo r = 0,9, e máximo r = 0,5. Por outro lado, quando o comprimento do passo é rejeitado, sendo esta a primeira tentativa desde a última aceitação, o fator r é calculado pela Equação (13) e multiplicado por 0,9. Após uma segunda tentativa sem êxito, o passo de integração é reduzido por um fator 0,25. Depois de três tentativas sem sucesso, a ordem é reduzida a 1 e o comprimento do passo é reduzido por um fator de 0,25 sempre que o teste falhar. Se o comprimento do passo for tal que t + Dt » t, o processo é interrompido. Este comprimento do passo Dt_min, representa o comprimento da condição de t + Dt_{j + 1}» t, e é calculado através da seguinte equação:

Onde representa o erro de arredondamento e T_Final o tempo final desejado.

3.2 Método de Diferenciação Regressiva Modificado Estendido

Os métodos BDF têm sido considerados os mais adequados para a solução numérica de problemas de valor inicial de sistemas de EDAs rígidas, mas existem na literatura novas metodologias que podem oferecer vantagens sobre os métodos BDF convencionais. Uma dessas metodologias é o método BDF modificado e estendido (MEBDF), apresentando características particularmente eficientes durante o processo de solução de sistemas de EDAs com índice menor e igual a três (Cash, 2000).

Uma das características mais relevantes da metodologia MEBDF está relacionada com a sua formulação modificada e estendida, relativa à metodologia BDF convencional. A formulação da metodologia MEBDF é proposta para alcançar ordens de convergência menores ou iguais a quatro, mantendo sua característica A-estável. Portanto, a metodologia é considerada A(a)-estável para ordens menores ou iguais a nove e Rígida-estável para ordens maiores que nove. Estas propriedades de estabilidade são consideravelmente melhores que as atingidas pelos métodos baseados nas fórmulas BDF convencionais (Cash, 1983; Cash e Considine, 1992; Cash, 2000).

O método MEBDF também apresenta a estratégia previsor-corretor, com duas etapas no estágio previsor e uma no estágio corretor. Na primeira etapa do previsor é calculado o ponto

_{j + k} utilizando-se uma fórmula BDF dada pela Equação (15) de ordem k baseada na interpolação dos pontos solução y_j, y_{j + 1},..., y_{j + k –}, sendo esta uma solução prévia usada nas etapas seguintes.. Esta etapa é muito importante porque constitui uma excelente aproximação para se conseguir a solução na etapa corretor. A matriz de iteração do esquema Newton Modificado tem a forma apresentada na Equação (16), onde I é a matriz identidade e J representa a matriz Jacobiana.

Na segunda etapa previsor é usada uma fórmula BDF convencional dada pela Equação (17) para a interpolação dos pontos solução y_{j + 1}, y_{j + 2}, ..., y_{j + k – 1}, _{j + k} (k pontos), e calcula um ponto super-futuro _{j + k + 1} e sua derivada '_{j + k + 1} . Esta solução servirá como uma aproximação ótima para o intervalo seguinte (t_{j + k + 1}) com relação a solução da etapa um do esquema previsor. A matriz iteração G^Pdo Newton Modificado será a mesma da primeira etapa. A matriz Jacobiana J é calculada e formada quantas vezes for necessário para se obter uma convergência ótima, embora o ideal seja mantê-la constante sempre que possível, para ser usada nos intervalos de tempo subseqüentes e desta forma amenizar o esforço computacional.

Com os pontos solução y_j, y_{j + 1}, ..., y_{j + k – 1}, _{j+ k}, _{j + k + 1}, '_{j + k + 1} obtidos no estágio previsor, calcula-se no estágio corretor a solução aproximada y_{j + k} para o ponto t_{j + k,} utilizando uma nova fórmula BDF modificada e estendida de ordem de convergência k + 1 dada pela Equação (18). Uma representação gráfica da estratégia previsor-corretor é mostrada na Figura 2, e uma comparação das regiões de estabilidade no plano complexo entre os métodos BDF e MEBDF de ordem 2 e 3 é apresentada na Figura 3.

A estabilidade e a precisão desta metodologia dependem muito das fórmulas usadas no estágio previsor, especificamente na primeira etapa. As fórmulas previsor devem ser pelo menos de ordem k se todo o processo previsor-corretor for de ordem de convergência k + 1.

Seleção do Passo de Integração

A metodologia mantém um registro dos comprimentos dos passos prévios que foram aceitos e que foram bem sucedidos no processo de integração. Esta informação é usada para decidir qual procedimento a ser empregado quando um passo de integração for recusado. Pode acontecer que, ao tentar incrementar o comprimento do passo, este processo tenha falhado duas vezes consecutivas e não consiga passar pelo teste de erro. Então, a estratégia assumida é calcular um novo comprimento de passo ótimo e continuar tentando passar no teste (Cash e Considine, 1992).

Para evitar a repetição de uma nova rejeição do comprimento de passo, a metodologia adota uma estratégia diferente que consiste em escolher um novo passo calculado de acordo com a Equação (19), na qual Dt é novo passo, Dt_ótimo é o passo estimado e Dt_êxito é o último passo bem-sucedido.

4 MÉTODO BASEADO NA FÓRMULA DE QUADRATURA RADAU

Rodolphe Radau (Radau, 1880) desenvolveu as conhecidas fórmulas de quadratura, enfatizando os métodos Gauss, Lobatto e Chebyshev, conhecidos atualmente como "fórmulas RADAU". As fórmulas de quadratura foram estendidas para serem usadas na solução de sistemas de EDAs.

As referências (Crank e Nicolson, 1947; Curtiss e Hirschfelder, 1952; Fox e Goodwin, 1949; Loud, 1949; Dahlquist, 1963) apresentam estudos das principais equações diferenciais e algébricas rígidas, e as referências (Butcher, 1964; Ehle, 1969; Axelsson, 1969) desenvolveram novas metodologias de integração numérica baseadas na teoria dos métodos de passo simples RUNGE-KUTTA totalmente implícitos. Especificamente, (Butcher, 1964) aborda métodos implícitos RUNGE-KUTTA baseado nas fórmulas de quadratura de Radau.

Estas fórmulas de quadratura tornaram possível a evolução dos métodos RUNGE-KUTTA que não eram totalmente implícitos (só permitiam que apenas o primeiro ou o último estágio fosse implícito). Originou-se uma nova família de métodos chamada de "Processos II", e posteriormente (Ehle, 1969) e (Axelsson, 1969) foram responsáveis pelas mofificações que tornaram os métodos RUNGE-KUTTA implícitos baseados nas fórmulas de quadratura RADAU, em métodos com propriedades A-estável, conhecidos na literatura como métodos RADAU IIA. (Hairer e Wanner, 1999).

Estratégias de Seleção de Ordem e Controle de Passo de Integração

A metodologia RADAU IIA apresenta duas estratégias de mudança de passo de integração baseadas na teoria de extrapolação de Richardson (Hairer e Wanner, 1996). Primeiro é calculado uma aproximação

_n de ordem S através da Equação (20), a partir dos últimos resultados (Y₁, Y₂, ..., Y_S) obtidos após a convergência do Newton Simplificado.

O erro é calculado pela Equação (21), onde o valor de corresponde ao autovalor real da matriz A^–1 para que a fatoração da matriz M – h.g₀ .J esteja disponível do esquema Newton Simplificado previamente convergido e possa ser utilizada no processo de seleção do comprimento de passo. A matriz A representa a matriz de coeficientes a_ij e M é uma matriz constante (possivelmente singular) e torna o sistema explícito quando diferente da matriz identidade.

A primeira estratégia para o cálculo do comprimento do passo necessita do erro e do comprimento do passo do intervalo já convergido (n-ésimo passo) sendo calculado pela Equação (22). O valor do fator fac depende do número máximo de iterações Newton Simplificado (k_max) e do número de iterações do último Newton Simplificado convergido (New).

A segunda estratégia calcula o passo utilizando informação dos últimos dois passos convergidos. Esta estratégia é chamada de "controle de passo com memória" e o passo é calculado de acordo com a Equação (23).

A metodologia RADAU IIA avalia as duas estratégias, calcula os dois possíveis comprimentos do passo e escolhe o menor, mas quando os passos prévios apresentam uma tendência a aumentar, o novo passo é calculado usando a Equação (22), e quando a tendência é a diminuir, o passo é calculado usando a Equação (23).

Quando uma metodologia oferece a possibilidade de mudança de ordem é necessário implementar uma estratégia para realizar eficientemente esta função. Basicamente, consiste em selecionar a ordem de tal forma que o erro por unidade de comprimento de passo seja mínimo. O cálculo do erro de truncamento é uma tarefa muito complicada e de difícil implementação. Um comportamento ineficiente da estratégia poderia afetar o cálculo do novo comprimento de passo uma vez que o número de iterações do esquema Newton Simplificado é dependente da ordem escolhida.

A solução de sistemas de EDAs de índice 1 para os métodos RUNGE-KUTTA implícitos não apresenta nenhum tipo de dificuldade, mas para sistemas de EDAs de índices superiores são necessários métodos com propriedades de mudança de ordem em cada estágio para facilitar a convergência, como é o caso dos métodos RADAU IIA (Hairer et al., 1989).

5 APLICAÇÃO EM SISTEMAS ELÉTRICOS

É apresentada uma estrutura geral do modelo representando o sistema de energia elétrica para a análise de fenômenos de estabilidade no domínio do tempo. O conjunto equações que descreve o modelo (as variáveis estão descritas no ^Apêndice) para estudos de curta e de longa-duração pode ser expresso na seguinte forma (Pessanha, 1997; Cutsen, 1993; Hiskens, 1995; Löf, 1995; Cutsen et al., 1994; Quintana e Vargas, 1994; Carvalho Mendes, 1997):

Reduções e simplificações podem ser realizadas nestas equações em função do objetivo do estudo bem como dos mecanismos envolvidos. A solução requer um tratamento desacoplado do sistema de EDAs, reduzindo-o a um sistema explícito de EDOs podendo comprometer a característica esparsa do sistema e a simetria da matriz Jacobiana. A formulação implícita das EDAs tem sua justificativa no bom desempenho do tratamento de matrizes não-simétricas, característica própria da estrutura de EDAs após a discretização do sistema de EDOs.

Para efeitos da aplicação e compatibilidade das formulações matemáticas dadas pelo conjunto de Equações (24), estas podem ser reescritas na seguinte forma implícita:

A estrutura do sistema de EDAs resultante da análise de estabilidade apresenta características não-simétricas por natureza, cuja configuração tem uma disposição tipo blocos diagonais. Os blocos contendo as equações das máquinas síncronas são acoplados com a rede de transmissão do sistema. Entretanto, a matriz de admitância da rede é de grande-porte, complexa e altamente esparsa.

5.1 Testes Computacionais

Quatro técnicas numéricas são testadas através de simulação computacional no domínio do tempo de fenômenos de estabilidade em sistemas de energia elétrica, sendo estas: MEBDFSD (com esparsidade), MEBDF (implícito), BDF (DASSL) e RADAU IIA. Estas técnicas estão disponíveis sob a forma de solvers numéricos (códigos computacionais). (Mazzia e Iavernaro, 2003).

A fim de capturar nas simulações os efeitos relevantes aos fenômenos de interesse, foram considerados os seguintes modelos (alguns foram retirados da literatura e outros foram desenvolvidos (Paz, 2004)):

• Máquina Síncrona - Modelos Clássico e IV.

• Regulador Automático de Tensão.

• Sinais Estabilizadores.

• Limitador de Sobreexcitação.

• Transformador de Tape Variável (tipo contínuo).

• Carga Estática tipo polinomial (parcela de potência, corrente e impedância constante).

• Carga dinâmica tipo potência constante.

• Acréscimo contínuo de carga ativa, reativa ou aparente.

• Acréscimo discreto de carga ativa, reativa ou aparente.

Caso I: Estabilidade Transitória Angular

Para avaliar a eficiência das técnicas de solução das EDAs para um fenômeno dinâmico de curta-duração utilizou-se o sistema-teste IEEE 118 barras com 54 geradores. O modelo IV de máquina síncrona, que incluí efeitos transitórios e subtransitórios, amortecimento e saturação, é usado para representar cada um dos 54 geradores síncronos, considerando-se também quatro modelos diferentes de curvas de saturação. São incluídos também dispositivos de controle como reguladores automáticos de tensão e sinais estabilizadores.

O evento consiste na aplicação de um curto-circuito trifásico sólido numa determinada barra em t=1 seg sendo removido 7 ms após sua aplicação. O tempo total de simulação é de 10 segs. A Tabela 1 ilustra as informações numéricas relacionadas a esta simulação.

Os resultados obtidos estão ilustrados sob a forma de gráficos nas Figuras 4 (a)-(c) oferecendo uma comparação entre os resultados obtidos por cada técnica. As grandezas selecionadas para comparação são aquelas de maior interesse e que ajudam na validação da eficiência numérica das metodologias. Todas as curvas apresentam uma excelente concordância sendo que a metodologia MEBDFSD é mais eficiente no aspecto esforço computacional de acordo com a Figura 5 que mostra um diagrama de barras para o tempo de CPU.

Imediatamente após o distúrbio, os comprimentos do passo de integração de todas as metodologias apresentam um comportamento moderado. A ordem varia com maior freqüência para alguns métodos a fim de se obter precisão e, assim, representar o comportamento oscilatório do sistema de forma adequada.

Quando o sistema tende a atingir uma condição de equilíbrio, as metodologias iniciam o processo de aumento do comprimento do passo de integração, de acordo com cada estratégia de controle de mudança de passo. O DASSL (BDF) apresenta um comportamento mais conservador, mantendo quase constante o comprimento do passo durante o restante da simulação. Já o RADAU, baseado no RK4, aproveita ao máximo sua estratégia de variação de passo. O MEBDF e o MEBDFSD devido a estratégia baseada no registro dos comprimentos dos passos de integração prévios, é ainda mais conservadora, fixando o tamanho de passo calculado no último passo bem sucedido sem arriscar a convergência do processo de solução, garantindo dessa forma a fidelidade e precisão da solução.

Caso II: Estabilidade de Tensão de Longa-Duração

Para avaliar a eficiência das técnicas de solução das EDAs para fenômenos dinâmicos de longa-duração, é utilizado o sistema-teste IEEE 150 barras com 50 geradores. Foram incluídos nesta configuração cinco transformadores de tape variável e os sistemas de controle de excitação dos geradores (regulador automático de tensão) passaram a ser monitorados por limitadores de sobrexcitação. Com isso, pretende-se capturar os efeitos relevantes ao fenômeno da estabilidade de tensão de longa-duração e aumentar a complexidade do processo de simulação.

O sistema é submetido a uma primeira contingência, caracterizada pela inserção de um reator de 20 Mvar em uma determinada barra no instante t = 5 segundos, resultando numa redução no perfil de tensão nesta barra, provocando a ação do transformador de tape variável 35 segundos após seu dispositivo de controle detectar níveis de tensão fora da faixa desejada. As Figuras 6(a) e ^(b) ilustram, respectivamente, a dinâmica do tape restaurando a tensão da carga para níveis pré-distúrbio e as dinâmicas lentas associadas a ação dos limitadores de corrente de campo. Já as Figuras 6(c) e 6(d) ilustram o comportamento do passo de intergração e da ordem de cada método, respectivamente.

No instante t = 200 segundos, o sistema é submetido a um aumento progressivo de potência reativa em cinco barras do sistema com diferentes taxas de incremento. A medida que a demanda de potência reativa aumenta, os geradores vão atingindo seus limites de capacidade de geração de potência reativa, limite este imposto pelo limitador de sobreexcitalção. Outros geradores vão atingindo seus limites e o sistema passa o enfrentar problemas de instabilidade de tensão, onde, por conseqüência, a ocorrência do colapso de tensão é iminente.

A Tabela 2 apresenta informações para uma análise do desempenho computacional das metodologias, destacando-se as informações inerentes ao número de avaliações da matriz Jacobiana e o tempo de CPU. A Figura 7 possibilita uma comparação entre os tempos de CPU gastos por cada metodologia. Observa-se a ineficiência computacional do RADAU para este caso em particular. Já as técnicas BDF e MEBDF apresentam um bom desempenho, destacando-se o MEBDFSD entre todas as testadas.

5.2 Análise do Desempenho Computacional das Técnicas

A análise comparativa de desempenho está baseada em três indicadores relativos as seguintes habilidades: avaliação das funções F, avaliação da matriz Jacobiana, fatoração da matriz Jacobiana e um último processo conhecido na literatura como Backward-Fordward Substitution. Esta última corresponde a solução do sistema de equações lineares da forma A.x = b. Associados a estas habilidades, têm-se: tempo gasto por execução do processo, o número de chamadas do processo, e a porcentagem do tempo gasto do processo na simulação global.

O Tempo Gasto por Processo (TGP) indica o tempo computacional médio gasto por uma metodologia para avaliar um determinado processo. Este indicador é calculado através da Equação (26).

T_i representa o tempo parcial gasto por um determinado processo a cada avaliação. Assume-se por cada j-ésimo intervalo os tempos parciais das "i" avaliações aceitas ou rejeitadas, NA_j representa o número de avaliações do processo por intervalo de tempo e NI é o número de intervalos de tempo avaliados. O número de chamadas indicará quantas vezes uma metodologia precisou realizar determinado processo durante toda a simulação. A porcentagem sobre o tempo total gasto é um índice de eficiência de cada técnica.

As Tabelas 3 e 4 apresentam os indicadores quantitativos resultantes da simulação computacional do Caso I e do Caso II, respectivamente, referentes a cada técnica numérica.

Como se pode observar, em todos os casos, o processo de fatoração exigiu maior esforço computacional em comparação aos demais. Para as metodologias BDF e MEBDF, o processo backward-forward substitution representa no resultado global o processo que consome maior porcentual do tempo total da simulação devido ao grande número de avaliações. Na metodologia RADAU o processo mais exigente é a fatoração da matriz Jacobiana, sendo ineficiente para simular sistemas de energia elétrica considerados de médio porte ou grande porte, devido à hiper-dimensão da matriz Jacobiana (5, 9 ou 13 vezes maior em comparação aos demais).

Observa-se que o número de avaliações nos processos é maior para o MEBDF em relação ao BDF, sendo esta última, portanto mais eficiente. Entretanto, o MEBDFSD continua sendo a melhor alternativa em termos de eficiência, apresentando tempos gastos por processo muito pequenos.

6 CONCLUSÕES

Este trabalho investigou o desempenho de técnicas numéricas aplicadas na solução de sistemas de equações diferenciais e algébricas, enfocando o desempenho computacional de cada uma delas em simulações computacionais de fenômenos de estabilidade em sistemas de energia elétrica (transitória angular e de tensão de longo-termo).

A metodologia RADAU baseada nas fórmulas de quadratura apresentou o maior tempo gasto durante o processo de fatoração da matriz Jacobiana aumentando o esforço computacional e tornando a metodologia ineficiente em termos de tempo de processamento.

As estratégias implementadas na metodologia BDF resultaram num melhor desempenho em comparação à metodologia MEBDF, realizando um número menor de processos backward-forward substitution, de avaliação e fatoração da matriz Jacobiana, sendo adequada para estudos de fenômenos de estabilidade em sistemas de energia elétrica onde o tempo de processamento é uma das principais preocupações.

As estratégias do MEBDF são mais conservadoras em relação aos métodos BDF uma vez que priorizam a convergência e precisão da solução. Estas estratégias combinadas com técnicas de esparsidade (MEBDFSD) as tornam muito robustas, eficientes e atraentes para serem usadas em estudos de estabilidade no domínio do tempo.

Das técnicas numéricas investigadas aqui, conclui-se que o MEBDFSD, MEBDF e o BDF são as mais indicadas para solucionar sistemas rígidos de EDAs, destacando-se o MEBDFSD, principalmente se o tempo de processamento for o principal interesse. Esta metodologia foi muito eficiente em termos de rapidez computacional em comparação as demais.

Artigo submetido em 01/07/2005

1a. Revisão em 11/09/2005;

Aceito sob rec. do Ed. Assoc. Prof. Carlos Alberto de Castro Jr.

Y: Vetor m-dimensional contendo as variáveis de estado dinâmicas, por exemplo; deslocamentos angulares dos rotores dos geradores, velocidades angulares, tensões de saída dos reguladores automáticos de tensão, escorregamentos dos rotores dos motores de indução, etc.

X: Vetor n-dimensional contendo as variáveis de estado algébricas, por exemplo; componentes de eixo direto e em quadratura do estator da máquina síncrona, potência ativa e reativa injetadas, amplitude e ângulo das tensões nas barras do sistema, etc.

Z: Vetor q-dimensional contendo as variáveis que sofrem alterações através de passos discretos, por exemplo; relação de espiras dos transformadores com troca automática ou manual de tapes, corrente de campo das máquinas síncronas monitorada pelo limitador de sobre-excitação (oxl), etc.

W: Vetor r-dimensional contendo as variáveis que evoluem ao longo do tempo, como por exemplo; potência ativa e reativa consumida pelas cargas, potência ativa programada para as unidades geradoras, potências ativas de intercâmbio entre áreas, etc.

U: Vetor s-dimensional contendo as variáveis de controle independentes, como por exemplo: potência ativa dos geradores, tensão de referência dos reguladores de tensão, tape dos transformadores controladores, etc.

P: Vetor s-dimensional correspondente aos parâmetros do sistema, como por exemplo: constantes de tempo e reatâncias das unidades geradoras, impedâncias e susceptâncias das linhas de transmissão, parâmetros dos modelos de cargas, etc.

F: Vetor de funções não-lineares que descreve as equações dinâmicas das máquinas síncronas, dos reguladores automáticos de tensão, dos sistemas de excitação, dos reguladores de velocidade, das cargas, dos compensadores estáticos, etc.

G: Vetor de funções não-lineares que descreve as equações algébricas de restrições impostas pela rede e também às equações algébricas das máquinas síncronas, dos reguladores de tensão, dos reguladores de velocidade, dos sistemas de excitação, etc.

H: Vetor das funções que expressam a natureza discreta no tempo das variáveis z.

f: Vetor das funções não-lineares temporais associadas aos ciclos de carga, das programações de geração e de intercâmbio, etc.

K: Tempo ou situação correspondente a um evento discreto.

T: Tempo em segundos, minutos ou horas.

• Abdulla, T. J; Cash, J. R. and Diamantakis, M. T. (2000). An MEBDF Package for the Numerical Solution of Large Sparse Systems of Stiff Initial Value Problems, Computers and Mathematics with Applications 42 (2001) 121-129.
• Ascher, U. M. and Petzold, L. R. (1998). Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Philadelphia, PA: SIAM Press.
• Astic, J. Y; Bihain, A. and Jerosolimski. (1994). The Mixed Adams-BDF Variable Step Size Algorithm to Simulate Transient and Long-Term Phenomena in Power Systems, IEEE Transactions on Power Systems, Vol.9, No. 2.
• Axelsson, O. (1969). A class of A-stable methods, BIT 9, 185-199.
• Brenan, K. E; Campbell, S. L. and Petzold, L. R. (1996) The Numerical Solution of Initial Value Problems in Differential-Algebraic Equations, SIAM Classics Series.
• Butcher, J. C. (1964). Integration Processes Based on Radau Quadrature Formulas, Math. Comput. 18, 233-244.
• Byrne, G. D and Hindmarch, A. C. (1975). A Polyalgorithm for the Numerical Solution of Ordinary Differential Equations, ACM Transactions on Mathematical Software, vol. 1, nş 1,pp. 71-96.
• Carvalho Mendes, P. P. (1997). Aplicaçăo de Redes Neurais Artificiais na Análise em Tempo Real da Estabilidade de Tensăo de Regime-Permanente de Sistemas Elétricos de Potęncia, Proposta de Tese de Doutorado, COPPE-UFRJ.
• Cash, J. R. (1983). The Integration of Stiff Initial Value Problems in EDOs Using Modified Extended Backward Differentiation Formulae, Comp. Math. Appl. 9, 645 - 660.
• Cash, J. R. and Considine, S. (1992). An MEBDF Code for Stiff Initial Value Problems, ACM Trans. Math. Software 18, 142 - 160.
• Cash, J. R. (2000). Modified Extended Backward Differentiation Formulae for the Numerical Solution of Stiff Initial Value Problems in EDOs and DAEs, Comput. Math. 125. 117 - 130.
• Crank, J. and Nicolson, P. (1947). A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat-Conduction Type, Proc. Cambridge Philos. Soc. 43, 50-67.
• Curtiss, C. F. and Hirschfelder, J. O. (1952). Integration of Stiff Equations. Proc. Nat. Acad. Sci., vol. 38, pp.235-243. [IV.1]
• Cutsen, T. V. (1993) Analysis of Emergency Voltage Situations, Proc. 11th Power Systems Computation Conference, Avignon, France, Vol.1, pp. 323-330.
• Cutsen, T. V; Jacquemart, Y; Marquet, J. N and Pruvot, P. (1994). A Comprehensive Analysis of Mid-term Voltage Stability, IEEE Trans. on Power Systems, SM 511-6, summer meeting.
• Dahlquist, G. (1963) A Special Stability Problem for Linear Multistep Methods. BIT, vol. 3, pp. 27-43.
• Ehle, B. L. (1969). On Padé Approximations To The Exponential Function and A-stable Methods For The Numerical Solution Of Initial Value Problems, Report CSRR 2010, Dept. AACS, Univ. of Waterloo, Ontario, Canada.
• Fábián, G; Van Beek, D.A. and J. E. Rooda. (2000). Substitute Equations for Index Reduction and Discontinuity Handling, In Proc. of the Third International Symposium on Mathematical Modelling, Vienna.
• Fox, L. and Goodwin, E. T. (1949). Some New Methods For The Numerical Integration Of Ordinary Differential Equations, Proc. Cambridge Philos. Soc. 45, 373-388.
• Gear, W. (1971). Numerical Initial Value Problems in Ordinary Diferential Equations.
• Hairer, E; Lubich, C and Roche, M. (1989). The Numerical Solution of Diferential-Algebraic Systems by Runge-Kutta Methods, Lecture Notes in Maths, Vol. 1409, Springer, Berlin.
• Hairer, E. and Wanner, G. (1996). Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, 2nd revised Edition, Springer Series in Comput. Math., Vol. 14, Springer, Berlin, 614pp.
• Hairer, E. and Wanner, G. (1999). Stiff differential equations solved by Radau methods, Journal of Computational and Applied Mathematics 111, 93-111 ELSEVIER.
• Hiskens, I. (1995). Analysis Tools for Power systems Contending with Nonlinearities, IEEE Proceedings, vol. 83, no.11, pp. 1573-1587.
• Jackson, K. R. and Sacks-Davis, R. (1980). An Alternative Implementation of Variable Step Size Multistep Formulas for Stiff ODEs, ACM Trans, Math. Software, 6, 295-318.
• Jardim, J. L. (1997). Utilizaçăo de Ferramentas de Simulaçăo Dinâmica de Longa Duraçăo na Análise de Fenômenos de Colapso de Tensăo e no Treinamento de Operadores, XIV SNPTEE, Belém, Pará, outubro.
• Löf, P.A. (1995). On Static Analysis of Long-Term Voltage Stability in Electric Power Systems, Ph.D. Thesis, Kungl Tekniska Högskolan - Royal Institute of Technology, Stockholm, Sweden.
• Loud, W. S. (1949). On The Long-Run Error In The Numerical Solution Of Certain Differential Equations, J. Math. Phys. 28 (1), 45-49.
• Mazzia, F. and Iavernaro, F. (2003). Test Set for Initial Value Problem Solvers, Department of Mathematics University of Bari ITALY, Report 40.
• Paz, A. R. A (2004). Implementaçăo de um Simulador Numérico Num Programa Computacional de Estabilidade, Dissertaçăo de Mestrado, CPGEE, UFMA, Fevereiro.
• Pessanha, J. E. O. (1997). Análise do Fenômeno da Estabilidade de Tensăo no Domínio do Tempo: Simulaçăo dos Períodos Transitórios e de Longo-Termo, Tese de Doutorado, Departamento de Engenharia Elétrica, PUC-Rio.
• Petzold, L. R. (1983). A Descriptions of DASSL: A differential/ algebraic system solver, in Scientific Computing, R. S. Stepleman et al., eds., North-Holland, Amsterdam, 65-68.
• Petzold, L. R. and LI, S. (2000). Software and Algorithms for Sensitivity Analysis of Large-Scale Differential Algebraic Systems, Journal of Computational and Applied Mathematics 125, 131-145, Elsevier.
• Quintana, V. H. and Vargas, L. (1994). Voltage Stability as Affected by Discrete Changes in the Topology of Power Networks, IEE Proceedings Generation, Transmission and Distribution, Vol. 141, no. 4, pp. 346-352.
• Radau, R. (1880). Étude Sur les Formulas D'approximation qui Servent ŕ Calculer la Valeur Numérique D'une Intégrale Définie, J. Math, Pures Appl. Sér. 3, 6, 283-336.
• Secanell, M. and Córcoles, F. (2002). DAEs implementation of dynamic power systems, IEEE, 663-669.
• Shampine, L. F. and Gordon, M. K. (1975). Computer Solution of Ordinary Diferential Equations, W.H. Freeman and Co., San Francisco, CA.
• Tarraf D. C. and Asada, H. H. (2002). On the Nature and Stability of Differential-Algebraic Systems, Proceedings of the American Control Conference Anchorage, AK May 8-10.

Apêndice

Datas de Publicação

• Publicação nesta coleção
  17 Abr 2006
• Data do Fascículo
  Set 2005

Histórico

• Revisado
  11 Set 2005
• Recebido
  01 Jul 2005