SciELO - Brazil - Determinação do ponto de máximo carregamento de sistemas de potência utilizando o fluxo de carga desacoplado rápido parametrizado Determinação do ponto de máximo carregamento de sistemas de potência utilizando o fluxo de carga desacoplado rápido parametrizado
Sistemas de Potência • Sba Controle & Automação 14 (2) • Jun 2003 • https://doi.org/10.1590/S0103-17592003000200007 copiar
Determinação do ponto de máximo carregamento de sistemas de potência utilizando o fluxo de carga desacoplado rápido parametrizado
L.C.P. da Silva
V.F. da Costa
Métodos de Continuação; Fluxo de Carga Desacoplado Rápido; Ponto de Máximo Carregamento
Continuation Methods; Fast Decoupled Load Flow; Maximum Loading Point
D.A. AlvesI; L.C.P. da SilvaII; C. A. CastroII; V.F. da CostaII
IUNESP - Departamento de Engenharia Elétrica, Av. Brasil, 56, Cx. P. 31, 15385-000 Ilha Solteira, Brasil
IIUNICAMP - Faculdade de Engenharia Elétrica e de Computação, Cx. P. 6101, 13083-970, Campinas, Brasil
Os métodos de fluxo de carga por Newton-Raphson e fluxo de carga desacoplado rápido convencionais são considerados inadequados para a obtenção do ponto de máximo carregamento de sistemas de potência, devido à problemas de mal-condicionamento neste ponto crítico e na sua vizinhança. Neste ponto a matriz Jacobiana do método de Newton-Raphson torna-se singular e considera-se que não são mais válidas as hipóteses de desacoplamento P-V e Q-q utilizadas para a formulação do método fluxo de carga desacoplado rápido. No entanto, mostra-se neste trabalho, que com pequenas modificações, as versões XB e BX do fluxo de carga desacoplado rápido tornam-se adequadas para a obtenção do ponto de máximo carregamento. Estas novas versões modificadas são comparadas entre si com o intuito de explicitar suas características, assim como da influência da atuação dos limites de geração de potência reativa e de tap's de transformadores. Os resultados obtidos para os sistemas testes do IEEE (14, 30, 57 e 118 barras) mostram que as características de convergência das versões originais são preservadas. Além disso, durante o traçado das curvas PV, os diversos métodos podem ser comutados entre si possibilitando o cálculo de todos os pontos da curva com um número reduzido de iterações.
Palavras-chave: Métodos de Continuação, Fluxo de Carga Desacoplado Rápido, Ponto de Máximo Carregamento.
The conventional Newton and fast decoupled load flow methods are considered to be inadequate to obtain the maximum loading point due to ill-conditioning problems at and near this critical point. At this point the Jacobian matrix of the Newton-Raphson method becomes singular, and it is believed that the P-V and Q-q decoupling assumptions made for the fast decoupled load flow formulation no longer apply. However, as shown in this paper, with small modifications the fast decoupled load flow versions (XB and BX) become adequate for the computation of the maximum loading point. These new modified versions are compared to each other with the purpose of pointing out their features, as well as the influence of reactive power and transformer tap limits. The results obtained for the IEEE systems (14, 30, 57 and 118 buses) show that the convergence characteristics of the conventional versions are preserved. In addition, it is shown that these methods can be switched during the tracing process in order to efficiently determine all the PV curve points with few iterations.
Keywords: Continuation Methods, Fast Decoupled Load Flow, Maximum Loading Point.
Os métodos de Fluxo de Carga (FC) têm sido amplamente utilizados nos estudos do planejamento e da operação de sistemas elétricos. Através destes estudos são definidas, entre outras, as características nominais e o carregamento dos equipamentos, os requisitos de suporte de reativos para a manutenção do perfil de tensão em condições normais e de emergência. Mais recentemente, o FC tem sido usado para a análise estática da estabilidade de tensão. Essa análise tem-se tornado uma necessidade crítica para a operação dos sistemas de potência na medida que estes têm sido levados a operarem próximos de seus limites em decorrência do crescimento contínuo da demanda, associado a restrições econômicas e ambientais impostas à construção de novas linhas de transmissão e usinas de geração. A implantação da nova política de desregulamentação do setor elétrico tem levado a uma operação em um ambiente competitivo, exigindo a adoção de novas estratégias operacionais. Muitas vezes isso implica em um aumento da transferência de potência entre determinadas regiões da rede, podendo levar o sistema a operar próximo de seu limite de máxima transferência de potência. Consequentemente, o conhecimento preciso de quão distante o atual ponto de operação se encontra de seu limite de estabilidade tornou-se crucial para o operador. O operador precisa saber se para um dado distúrbio, como por exemplo a saída de uma linha de transmissão, a variação súbita do carregamento do sistema, ou o aumento da transferência de potência entre áreas, ainda existirá um ponto de operação factível para as novas condições.
As equações do FC representam um limite para a região de operação sendo, portanto, essenciais para o estudo da estabilidade estática. Se para uma condição especificada elas não apresentarem solução conclui-se que a geração e a rede não são fisicamente capazes de suprir a carga especificada, o que exigirá mudanças ou na geração ou no sistema de transmissão, ou em ambas (Lesieutre et alii, 1999).
As análises estáticas de estabilidade de tensão de sistemas de potência podem ser realizadas através da obtenção do perfil de tensão das barras em função de seu carregamento (curvas PV e VQ). Estas curvas possibilitam a compreensão das condições de operação do sistemas para diferentes carregamentos, e têm sido recomendadas pelas empresas do setor elétrico nacional (Força Tarefa, 1999) e internacional (Work Group, 1998), para avaliação da estabilidade de tensão. Entre outras aplicações, estes perfis são usados para determinar os limites de transferência de potência entre as áreas de um sistema, ajustar margens, observar o comportamento das tensões das barras do sistema em análise, e comparar estratégias de planejamento. Procedimentos automáticos, utilizando método de Newton-Raphson (MNR), vêm sendo desenvolvidos para se efetuar estas análises a partir das curvas PV.
Um dos principais objetivos do estudo da estabilidade de tensão em regime permanente é a obtenção do ponto de máximo carregamento (PMC) do sistema. A obtenção do PMC é importante tanto para o cálculo de margens de estabilidade quanto para a realização da análise modal, uma vez que é neste ponto que ela fornece informações para a determinação de medidas efetivas para o reforço do sistema, já que o PMC define a fronteira entre as regiões de operação estável e instável. Tradicionalmente, a obtenção de sucessivas soluções do FC tem sido feita através da variação manual do carregamento do sistema. Este procedimento é realizado até que o processo iterativo deixe de convergir. Para fins práticos, este ponto é considerado como sendo o PMC, embora se saiba que ele está próximo do PMC.
A solução do FC dependerá das características comuns aos processos de solução de equações algébricas não-lineares, tais como da existência da solução, das múltiplas soluções existentes, do método utilizado na resolução, da existência de singularidades, da estimativa inicial, etc. Em geral, a inicialização a partir de um caso base se deve à existência de múltiplas soluções. Por isso, o processo iterativo a ser utilizado tem que ser cuidadosamente inicializado para assegurar que uma solução apropriada seja obtida, uma vez que uma inicialização inadequada poderá resultar em não convergência, ou convergência para soluções não desejadas. Portanto, uma das características desejadas para o método a ser utilizado, seria a de possuir uma região de convergência tão ampla quanto possível, evitando assim a possibilidade de erros de convergência. Problemas de convergência encontrados pelo FC convencional para a obtenção do PMC, é um exemplo claro das dificuldades numéricas relacionadas à singularidade da matriz Jacobiana. Esta instabilidade numérica é observada mesmo com o uso de cálculo com dupla precisão ou de algoritmos anti-divergentes. Assim sendo, o uso dos métodos convencionais de FC para a obtenção das curvas PV fica restrito à sua parte superior (correspondendo a operação estável). Muito embora o uso de métodos de FC convencionais possibilitem o cálculo de pontos de operação muito próximos ao PMC, sempre será necessário ponderar se os problemas de não convergência são devidos a problemas numéricos ou a limitações físicas do sistema. Em geral, as diferenças não são óbvias.
Métodos de continuação vêm sendo utilizados por longa data na análise de sistemas de equações algébricas não-lineares parametrizadas (Seydel, 1994). Mais recentemente foi proposta a sua utilização para o traçado das curvas de carregamento, obtenção de múltiplas soluções e do ponto de colapso de tensão de sistemas de potência (Kundur, 1993; Van Cutsen e Vournas, 1998). O método de continuação supera as dificuldades de obtenção de sucessivas soluções do FC possibilitando o traçado completo do perfil de tensão variando automaticamente o valor do carregamento do sistema, sem preocupação com as singularidades das equações do sistema de potência. O método de Newton (Ajjarapu e Christy, 1992; Cañizares et alii, 1992) e o método desacoplado (Bijwe e Tare, 1997) são utilizados para a solução do problema parametrizado. A parametrização fornece uma forma de identificar cada solução ao longo da trajetória a ser obtida. Entre os diversos métodos de fluxo de carga continuados descritos na literatura, o mais amplamente utilizado consiste de quatro elementos básicos: um passo preditor, um procedimento de parametrização, um controle de passo, e um passo corretor. A reformulação das equações do FC visa a eliminação da singularidade da matriz Jacobiana no PMC e, consequentemente, dos problemas numéricos que ocorrem em torno deste.
A ampla utilização do método desacoplado rápido tradicional, no que diz respeito ao cálculo de FC, se deve às suas características de velocidade de convergência e pequeno espaço de memória necessário. Em função destas vantagens, há um interesse por parte das empresas na sua utilização nos algoritmos de continuação visando, com isso, a redução do tempo de CPU necessário para o traçado das curvas PV. A não utilização dos métodos desacoplados se deve ao fato de se acreditar que o desacoplamento não possa ser usado próximo ao PMC, e que somente as soluções de alta tensão podem ser obtidas.
Bijwe e Tare (Bijwe e Tare, 1997), apresentam um método de continuação baseado no desacoplamento da matriz Jacobiana, onde os efeitos das submatrizes não diagonais são levados em conta por meio de subesquemas iterativos. Uma vez que para condições de alto carregamento estas submatrizes não podem ser completamente ignoradas, seus efeitos foram transferidos para os seus respectivos mismatches.
Este trabalho apresenta as características de novos métodos de continuação baseado em pequenas modificações do método do fluxo de carga desacoplado rápido convencional (FCDR), para o traçado das curvas PV. Os métodos de FCDR parametrizados (FCDRP) utilizam qk ou Vk, ângulo de fase e magnitude de tensão de uma barra k, respectivamente, ou o fator de carregamento . Também se apresenta a possibilidade do uso da potência reativa de uma barra PV como parâmetro. No passo preditor, os métodos usam a técnica de previsão trivial conhecida como preditor polinomial modificada de ordem zero (Chiang et alii, 1995).
Observa-se aqui, que as características inerentes dos FCDR convencionais são preservadas. A aplicação deste novos métodos aos sistemas teste do IEEE (14, 30, 57, e 118 barras) mostram que o PMC pode ser obtido com a precisão desejada sem qualquer dificuldade numérica. Além disso, mostra-se que estes métodos podem ser comutados entre si durante o traçado da curva PV, permitindo com isso um baixo número de iterações durante o traçado da curva completa. A troca de parâmetro é feita quando o número de iterações excede um determinado valor preestabelecido. O processo de troca é realizado de forma que o método que apresenta o melhor desempenho seja escolhido em cada parte da curva.
Outro aspecto importante a ser ressalvado é que o PMC obtido, para seus respectivos sistemas, foi sempre o mesmo, independentemente da metodologia utilizada. Assim, mostra-se que a diferença do ponto de divergência dos métodos, quando da obtenção do PMC, se deve mais a sua inadequação para o cálculo do respectivo ponto, do que do método em si.
2 O TRAÇADO DA CURVA PV
As curvas PV podem ser obtidas por meio de sucessivas soluções de FC, a partir de um caso base até o PMC, para incrementos graduais da carga numa direção predefinida. Obtém-se, com isso, não só a margem de carregamento, mas também, informações adicionais a respeito do comportamento das tensões das barras do sistema em análise. Dependendo da direção do aumento de carga, diferentes PMC podem ser obtidos. Em geral as equações de FC podem ser escritas como
em que q é o vetor dos ângulos de fase nodais, V é o vetor das magnitudes de tensões nodais, é o fator de carregamento, Pesp = Pgen - Pcarga é a diferença entre as potências ativas gerada e consumida para as barras de carga (PQ) e de geração (PV), e Qesp = Qgen - Qcarga para as barras PQ. O sistema de equações (1) pressupõe que o carregamento da rede é proporcional ao do caso base e considera o fator de potência constante. Pesp e Qesp podem ser, respectivamente, definidos como sendo iguais a ( KpS cosj)espe ( KqS sinj)esp. Assim, é possível realizar uma variação de carregamento individual, isto é, para cada barra do sistema, considerando para cada uma, um crescimento de carga com fatores de potência diferentes aos do caso base. Tradicionalmente, entretanto, assume-se que o aumento de carga de uma determinada área é feito com fator de potência constante e proporcional ao carregamento do caso base com modelo de carga de potência constante, visto que este fornece a condição operacional mais segura para o sistema (WSCC, 1998).
O traçado da curva PV, utilizando o método de Newton convencional para a resolução do sistema de equações (1), é realizado incrementando-se gradualmente o valor de , a partir do caso base ( = 1), até um valor para o qual não mais se obtenha solução (o processo iterativo do FC não converge). Neste caso é tratado como parâmetro no processo iterativo de Newton. Entretanto, no procedimento utilizado pelos métodos de continuação em geral, é considerado como variável dependente e, portanto, varia automaticamente. Dessa forma, o número de variáveis do sistema de equações (1), cuja dimensão é n = 2nPQ+nPV , onde nPQ e nPV correspondem ao número de barras PQ e PV, respectivamente, passa a ser n + 1. A diferença entre os métodos de continuação está na forma de tratar esta nova variável e em como contornar a singularidade da matriz Jacobiana. No contexto do traçado das curvas PV, a adição de equações parametrizadas ou normalizadas têm sido um procedimento padrão nas análises numéricas (Ajjarapu e Christy, 1992).
Como anteriormente comentado, ao incluir-se como variável o sistema de equações (1) passa a ter (n + 1) variáveis para n equações. Assim, pode-se definir qualquer uma das n + 1 variáveis como parâmetro. O valor deste parâmetro (qk, Vkou até mesmo pode ser especificado e, portanto, esta variável pode ser removida do sistema. Com isso, o sistema passa a ser novamente de n equações a n incógnitas. Portanto, não é necessariamente preciso, acrescentar uma equação. A escolha do parâmetro em função daquele que apresentar a maior variação no vetor tangente elimina a singularidade da matriz Jacobiana. Esta característica de troca de variáveis corresponde à rotação de 90o do diagrama qk versus , ou Vk versus (Seydel, 1994). Este será o procedimento utilizado nas versões parametrizadas por qk e Vk.
2.1 Métodos de Newton Parametrizados
No caso em que é tratado como variável independente no processo iterativo de Newton, isto é, quando ele for usado como parâmetro e seu valor for prefixado, a linearização da equação (1) de acordo com o método de Newton fornece:
onde as submatrizes que compõem a matriz Jacobiana J são representadas por H=¶P/¶q, N=¶P/¶V, M=¶Q/¶q e L=¶Q/¶V. DP e DQ correspondem aos mismatches de potências ativa e reativa, respectivamente, enquanto DV e Dq correspondem às correções das magnitudes e ângulos das tensões.
A figura 1(a) mostra o diagrama unifilar de um sistema de 3 barras e 2 ramos, o qual será usado como exemplo para o traçado da curva PV. Conforme se pode verificar nas figuras 1(b) e 1(c), através de incrementos sucessivos de o ponto de operação do sistema (q,V) variará. O ponto A corresponde ao PMC, para o qual a matriz J torna-se singular, conforme se pode constatar na figura 1(d), na qual é apresentado o determinante de J em função do parâmetro . Observa-se que, com o cálculo das soluções da equação (2) para diferentes valores de , só é possível a obtenção de pontos da parte superior da curva PV até pontos muito próximos do PMC (ponto A).
2.1.1 Método de Newton Parametrizado por Vk
Para se superar as dificuldades numéricas devidas à singularidade de J no PMC, é necessária uma mudança de parâmetro de forma que a nova matriz Jacobiana não apresente singularidade no PMC, e de preferência em toda sua vizinhança.
De acordo com Arya e Verma (1996), a curva PV completa pode ser obtida considerando a magnitude da tensão como parâmetro, sendo considerado como variável dependente.
Ajjarapu e Christy (1992) mostraram que os ângulos da tensão também podem ser escolhidos como parâmetro. Conforme ilustra a figura 1, na curva PV a tensão está continuamente decrescendo à medida que o carregamento se aproxima do PMC. Assim, a tensão de uma determinada barra poderia ser prefixada e a respectiva solução determinada. Para alcançar este objetivo, uma possibilidade seria a utilização da magnitude de tensão de uma barra k, Vk, como parâmetro, e nesta condição a aplicação do método de Newton à equação (1) resultaria em:
onde o vetor DV' não contém o elemento DVk, o qual é substituído por D . As submatrizes N' e L' diferem das N e L da equação (2) apenas na coluna k na qual as derivadas de P e Q com relação a Vk foram substituídas pelas derivadas com relação à nova variável . Conforme pode-se ver na figura 1(d), procedendo desta forma a nova matriz Jacobiana JV não apresentará singularidade no ponto A, resultando num método mais robusto não só no ponto em questão, mas ao longo de toda a curva.
2.1.2 Método de Newton Parametrizado por qk
Uma outra possibilidade seria o uso do ângulo da tensão de uma barra k, qk, como parâmetro. Neste caso, a aplicação do método de Newton o sistema de equações (1) resultaria em:
em que o vetor Dq' não contém o elemento Dqk, o qual é substituído por D . As submatrizes H' e M' diferem das H e M da equação (2) apenas na coluna k onde as derivadas de P e Q com relação a qk foram substituídas pelas derivadas com relação à nova variável .
Conforme pode-se ver na figura 1(d), também neste caso a nova matriz Jacobiana Jq não seria singular no ponto A, mas somente no ponto B. Sendo este ponto localizado na parte inferior da curva, a singularidade de J no ponto A é evitada uma vez mais. Entretanto, caso se queira obter mais pontos da curva, seria preciso efetuar uma troca de parâmetros na vizinhança do ponto B, pois o método apresentaria dificuldades numéricas em torno desse ponto. Observa-se que se q3 fosse adotado como parâmetro, a matriz Jq não apresentaria singularidade ao longo de toda a curva.
2.1.3 Método de Newton Parametrizado pela potência reativa de uma barra PV (QPV), m
Alves et alii (Alves et alii, 2000) mostraram que a parametrização pela potência reativa de uma barra PV pode ser utilizada para contornar a singularidade da matriz J no PMC. A injeção de potência reativa gerada de uma barra PV qualquer pode fazer parte do vetor x =[q V]T por meio da inclusão da respectiva equação de balanço de potência reativa no sistema de equações (1), tomando o cuidado para que quando o limite de reativos for atingido, a correspondente magnitude de tensão seja considerada no vetor x, e seja efetuada a escolha de uma nova barra PV.
Visto que os programas de FC podem normalmente incluir as equações das barras PV (IEEE/PES, 1999), o uso da potência reativa de uma barra PV qualquer como parâmetro apresenta-se como uma opção muito atraente, já que isto praticamente implicaria em poucas alterações nos programas de FC convencionais. Não haveria nenhum problema relacionado com perda de esparsidade, apenas a substituição de uma coluna correspondendo à nova variável . Portanto, bastaria especificar o valor de Q, através do parâmetro m, de uma barra PV que esteja com folga de reativos. Este procedimento seria válido enquanto os reativos gerados permanecessem dentro de seus limites. Quando o limite fosse atingido, haveria a necessidade da escolha de uma outra barra PV ou a troca do tipo do parâmetro, por exemplo, por Vk. Neste procedimento a seguinte equação é acrescentada ao sistema de equações (1):
Nesta equação Q corresponde à equação da potência reativa da barra de geração k escolhida. Q0 é o respectivo valor da função no caso base. Um novo parâmetro m é adicionado ao problema. Como uma equação é adicionada, pode ser tratado como uma variável dependente e m é considerado como um parâmetro. Para m = 1, a solução convergida deverá resultar em = 1. Em (5) pode-se especificar o valor desejado para a função, através do parâmetro m, e a sua solução, juntamente com o sistema de equações (1), fornecerá o ponto de operação, incluindo o fator de carregamento para o qual os valores especificados ocorrem.
A expansão do sistema de equações (1), acrescido da equação (5), em série de Taylor, incluindo somente os termos de primeira ordem, resulta em:
onde =[PespQesp]T e DG=[ DP DQ]T. Qx e correspondem às derivadas de Q em relação a x e respectivamente. DW corresponde ao mismatche de potência reativa da barra de geração k escolhida. Observe que esta equação pode ser colocada na mesma forma da equação (3). Entretanto, deve-se observar que o sistema de equações em questão, ao contrário do empregado quando do uso de qk ou Vk como parâmetro, possui (n + 1) variáveis para (n + 1) equações. A nova variável m é tratada como uma variável independente, ou seja, seu valor é prefixado.
2.2 Métodos Desacoplados Parametrizados
As simplificações introduzidas por Stott e Alsaç (Stott e Alsaç, 1974) na matriz Jacobiana do FC pelo método de Newton deram origem ao método de FCDR, descrito por:
Esta versão é atualmente conhecida como versão XB. Outras versões foram implementadas e testadas por Van Amerongen (Van Amerongen , 1989), sendo uma delas a versão BX. A diferença básica entre as versões está nos elementos das matrizes B' e B". As matrizes são constantes e de dimensões diferentes em virtude da exclusão das barras do tipo PV da matriz B".
Monticelli et alii (Monticelli et alii, 1990) apresentaram uma justificativa analítica e unificada, que possibilitou uma melhor compreensão do bom desempenho dos métodos desacoplados rápidos (BX e XB), a qual era até então desconhecida. Mostrou-se que os algoritmos das versões BX e XB podem ser deduzidos a partir da equação do método de Newton (equação (2)), assumindo-se que o sistema seja radial ou então constituído por ramos com a mesma razão rkl/xkl, e considerando que V = 1 e q = 0, condição conhecida como flat-start (estimativa inicial onde todos os ângulos são assumidos iguais a zero e as tensões iguais a 1,0 p.u. (El-Abiad, 1983)).
Mesmo quando o sistema não seja constituído por ramos com a mesma razão rkl/xkl ou não seja radial, os algoritmos BX e XB continuam a ser uma boa aproximação para fins práticos. No geral, a versão BX apresenta melhor desempenho que a XB, em virtude das aproximações extras (desconsideração dos efeitos das barras PVs e dos elementos shunts na matriz de sensibilidade P-q) necessárias para a obtenção da última (Monticelli et alii, 1990).
Os métodos propostos do FCDR parametrizados (FCDRP) foram obtidos a partir dos respectivos métodos de Newton parametrizados, descritos no item 2.1 acima, usando o mesmo procedimento apresentado em (Monticelli et alii, 1990). O sistema considerado para a obtenção dos FCDRP pode ser visto na figura 1(a). Serão consideradas três versões: BX parametrizada por Vk, XB parametrizada por qk e BX parametrizada pela potência reativa de uma barra PV, QPV.
2.2.1 Versão BX Parametrizada por Vk
Considere inicialmente o sistema de 3 barras da figura 1(a). Utilizando V3 como parâmetro e considerando V2 = 1 p.u. e q2 = q3 = 0, a equação (3) torna-se:
em que (gkl + jbkl) = (rkl + jxkl)-1. Usando o mesmo procedimento da referência (Monticelli et alii, 1990), obtém-se a matriz Leq = L' - MH-1N':
Em (Monticelli et alii, 1990) mostra-se que mesmo no caso do sistema não ser radial ou não satisfazer a condição da relação rkl/xkl constante, os elementos da Leq ainda terão valores muito próximos dos da matriz B" da versão BX. Entretanto, se a relação rkl/xkl for constante para todo o sistema:
em que (G2 + jB2) e (G3 + jB3) são os elementos da diagonal da matriz admitância. Em (Alves et alii, 2000) mostra-se que esta mesma condição é válida para um sistema não radial que apresente uma relação rkl/xkl constante. Finalmente, a versão BX do FCDRP utilizando Vk como parâmetro é dada por:
em que o elemento DVk do vetor DV' é eliminado, sendo substituído por D . A matriz B' é obtida pelo cálculo de H para o flat-start. A matriz B" é dada por Leq= L'- MH-1N' calculada para o flat-start. Além disso, os elementos fill-in são desprezados e bkl é substituído por 1/xkl. A coluna k de B" é substituída pelo seguinte vetor B":k
2.2.2 Versão XB Parametrizada por qk
Considere novamente o sistema de 3 barras da figura 1(a). Utilizando q3 como parâmetro e considerando V<2 = V3 = 1 p.u. e q2 = 0, a equação (4) torna-se:
Utilizando um procedimento similar ao da versão BX, para sistemas radiais e malhados com relação rkl/xkl constante, a versão XB parametrizada por qk do FCDRP será representada por:
onde a matriz B" é a submatriz L calculada para V=1 e q=0. Já a matriz B' é dada pela matriz equivalente Heq= H'- NL-1M' calculada para o flat-start, com as hipóteses adicionais propostas em (Monticelli et alii, 1990), isto é, desconsiderando os efeitos dos shunts, barras PV e tap de transformadores, que são representados na matriz L. Assim sendo, a coluna correspondente ao novo parâmetro escolhido, qk, é substituída pelo vetor B'k, cujos elementos são dados por:
2.2.3 Versão BX Parametrizada pela potência reativa de uma barra PV (QPV), m
Considere o sistema de 3 barras da figura 1(a), porém com a barra 3 send o agora uma barra PV. Usando a potência reativa desta barra como parâmetro e considerando V2 = 1 p.u. e q2 = q3 = 0, a equação (6) torna-se:
Utilizando um procedimento similar ao usado para obter as versões anteriores, para sistemas radiais e malhados com relação rkl/xkl constante a matriz B" desta versão BX, será dada por:
onde é a potência reativa consumida na barra de geração (PV). Assim, a versão BX do FCDRP utilizando m como parâmetro é dada por:
onde o elemento DVk do vetor DV' é eliminado, sendo substituído por D . As matrizes B' e B" são obtidas com procedimentos similares aos da versão BX do FCDRP utilizando Vk como parâmetro. A coluna de B", correspondente a barra de geração Qk cuja equação da potência reativa foi escolhida, é substituída pelo seguinte vetor B"Qk:
2.3 Passo Preditor e Controle de Passo
Uma vez que a solução do sistema de equações (1) para o caso base tenha sido encontrada (V0, q0, 0 = 1) pelo método convencional de FC, um dos métodos propostos pode ser usado para calcular soluções adicionais até que o PMC seja alcançado. O passo preditor é executado para encontrar uma estimativa para a próxima solução. Entre as diversas técnicas de previsão diferentes encontradas na literatura a do preditor secante (Chiang et alii, 1995) e a do preditor tangente (Ajjarapu e Christy, 1992) são os mais populares.
No preditor tangente, a estimativa da próxima solução pode ser encontrada dando um passo, de tamanho apropriadamente escolhido, na direção do vetor tangente à curva PV. A técnica do preditor pelo vetor tangente é usualmente mais precisa do que a da secante. Quando se está utilizando os métodos de Newton, o cálculo do vetor tangente não implicará num aumento significativo do custo computacional, já que se pode usar a última matriz Jacobiana fatorada. Já o preditor secante usa as soluções atual e anterior, para estimar a próxima. Portanto, este é uma aproximação do preditor tangente. Quanto mais próximos os dois pontos estiverem um do outro, tanto melhor será a aproximação do vetor secante com relação ao vetor tangente. Entretanto, isto implicaria no cálculo de um maior número de pontos da curva para se obter o PMC. Por outro lado, pontos muito distantes um do outro produziriam estimativas muito imprecisas, que poderiam dificultar a convergência do passo corretor.
O preditor secante tem como principais vantagens a implementação computacional relativamente mais simples e a inexistência de problemas relacionados à singularidade das matrizes Jacobianas, porém, tem como desvantagem o fato de resultar em algoritmos menos robustos, que requerem passos menores.
Por fim, uma outra técnica de previsão, a qual será adotada para o passo preditor em todos os métodos propostos nesse trabalho, é a denominada de previsão trivial ou polinomial modificada de ordem zero (Chiang et alii, 1995; Seydel, 1994). Esta técnica usa a solução atual e um incremento fixo no parâmetro ( , V, q, ou m) como uma estimativa para a próxima solução. Observa-se que a solução sucessiva do sistema de equações (1) pelo método de Newton com a prefixação de vários valores para o parâmetro em uso ( , qk, ou Vk) são casos particulares dessa técnica de previsão.
A eficiência de um método de continuação durante o traçado da curva PV está intimamente relacionada com a estratégia de controle de passo. No geral, sua escolha é dependente do sistema em estudo. Em situação de carga leve uma variação de carga resultará em uma pequena mudança no ponto de operação. Portanto, o tamanho do passo pode ser maior. Por outro lado, uma pequena mudança na carga resultaria em grandes variações do ponto de operação no caso do sistema estar altamente carregado, e nesse caso o passo deveria ser menor. Para se obter um bom desempenho global, o tamanho do passo deveria adaptar-se às condições reais de convergência.
Um método simples de controle de passo é baseado no número de iterações do passo corretor (Seydel, 1994). Se o número de iterações é pequeno, o sistema ainda está sob condições de carga leve ou normal e o tamanho do passo para o próximo ponto de operação pode ser maior. Se o número de iterações é grande, o sistema está estressado e o tamanho do passo para o próximo ponto de operação deve ser menor. Entretanto, este procedimento pode falhar na identificação da proximidade do maior carregamento já que o aumento do número de iterações do passo corretor também pode ser conseqüência da atuação de diversos limites operacionais de equipamentos.
Como se verá mais a frente, a parametrização por qk, ou Vk implica num controle de passo implícito em . Também será mostrado, que é possível deslocar-se de um ponto para outro mais eficientemente (com um número menor de iterações) através da mudança de parâmetro ao invés da redução do passo para o traçado da curva utilizando o FCDRP.
2.4 Passo Corretor
Uma vez especificado o incremento fixo no parâmetro ( , V, q, ou m) como uma estimativa para a próxima solução, torna-se necessário realizar-se a correção, a partir da solução atual, para obter-se a solução final. Em geral, o incremento adotado pelo passo preditor exige poucas iterações para que a próxima solução seja obtida dentro da precisão desejada. A solução correta é obtida pelas equações (2) e (7) no caso de ser adotado como parâmetro, pelas equações (3) e (11) no caso de Vk , pelas equações (4) e (14) no caso de qk , e pelas equações (6) e (18) no caso de potência reativa de uma barra PV (QPV) qualquer.
Para todos os testes realizados, a tolerância adotada para os mismatches é igual a 10-4 p.u.. A tensão com a menor magnitude, ou aquela com a maior taxa de decréscimo foi escolhida como parâmetro. Os passos de 0,01 e 10 foram adotados para a magnitude da tensão Vk e para o ângulo da tensão qk, respectivamente. O esquema de iteração padrão são usados em todos os exemplos que se seguem, para os métodos desacoplados. Entende-se por esquema de iteração padrão aquele que, após a convergência de um dos subproblemas, iterações sucessivas podem ser realizadas no outro subproblema, até que ambos atinjam a convergência ou, o que pode ocorrer, o subproblema convergido divirja (Stott e Alsaç, 1974).
O primeiro ponto de cada curva é obtido com o método de FC convencional. Os limites superior e inferior adotados para os taps de transformadores foram 1,05 e 0,95. O ajuste de tap nos transformadores OLTC (transformadores com comutação de tap sob carga) consiste da inclusão da posição do tap como variável dependente, ao passo que a magnitude da tensão das barras controladas são consideradas variáveis independentes (Peterson e Meyer, 1971).
O controle de limites de Q nas barras PV é o mesmo utilizado no método convencional de FC. Em cada iteração a geração de reativos de todas as barras PV são comparadas com seus respectivos limites. No caso de violação, a barra PV é alterada para tipo PQ. Estas barras podem voltar a ser PV nas iterações futuras. As violações de limite de tap também são verificadas. Observa-se que nas análises de colapso de tensão, um modelo que considere precisamente os limites de campo e de corrente de armadura da máquina síncrona deveria ser considerado. Contudo, neste trabalho o modelo proposto em (Cañizares, 1995) e adotado em (Zambroni et alii, 1997) é usado por simplicidade, isto é, as reatâncias internas dos geradores são desprezadas e a tensão terminal é assumida como sendo constante.
As cargas são modeladas como de potência constante e o parâmetro é usado para simular incrementos de carga ativa e reativa, considerando fator de potência constante. Cada aumento de carga é seguido por um aumento de geração equivalente usando . O objetivo dos testes é comparar os métodos propostos com o propósito de explicitar suas características, considerando a influência dos limites de geração de potência reativa (Q) e de tap dos transformadores.
3.1 Efeito dos Limites de Geração de Potência Reativa das Barras PV e dos tap de Transformadores no Desempenho dos FCDRP
A figura 2 ilustra os resultados dos testes para o sistema IEEE-30. Todas as curvas, exceto a da figura 2(d) foram obtidas usando a versão BX do FCDRP por Vk, usando V30 como parâmetro.
A figura 2(a) e a tabela 1 mostram o efeito dos limites de geração de potência reativa e do tap sobre o valor do fator de carregamento no PMC. Destes resultados pode-se concluir que a inclusão dos limites no controle de potência reativa, ao contrário dos limites de controle do tap, ocasionam uma grande redução no máximo carregamento. As figuras 2(b) e 2(c) mostram as curvas PV para a barra crítica e para duas barras PV, respectivamente. O efeito do limite de controle de potência reativa é evidente.
É importante notar nas figuras 2(b) e 2(d), conforme comentado no item 2.3, que a escolha de um passo constante para o parâmetro (V30 ou q30) resulta num passo variável para . Em particular, a variação diminui à medida que o carregamento se aproxima do PMC. Portanto, as parametrizações usando Vkou qk resultam num controle automático de passo para . Isto ocorre porque, em geral, um passo fixo na tensão corresponde a passos grandes na variável durante a carga leve e normal, onde a tensão varia pouco, e em passos reduzidos para altos carregamentos. Observa-se que esta redução intrínseca apenas facilita o cálculo do PMC com maior precisão, não sendo necessária ao longo de toda a curva PV para o caso do parâmetro Vk. Isto já não é o caso do parâmetro qk, que apresenta um ''nariz'' no ponto A, ver figura 2(d). Nesta figura mostra-se a curva Pq para a barra crítica obtida pela versão XB (q30 como parâmetro) e a versão BX (V30 como parâmetro). Quando q30 é utilizado como parâmetro, somente é possível obter a curva até a vizinhança do ponto A, onde a matriz Jacobiana modificada Jq torna-se singular. Neste ponto será necessária uma mudança de parâmetro se o objetivo for o de traçar a curva completa. Deve-se enfatizar que este problema também pode ocorrer quando se está usando V30 como parâmetro, isto é, a matriz JV, da mesma forma que a Jq, pode ser singular ou na parte inferior ou na superior da curva PV, dependendo da barra escolhida para o parâmetro.
3.2 Desempenho das versões XB e BX do FCDRP
A figura 3 mostra o desempenho da versão XB para o sistema IEEE-57, considerando o vetor calculado de acordo com a equação (15), e desconsiderando o termo (Gi/Bi) . O número de iterações necessárias para a convergência do FC é plotado para cada ponto da curva. A linha vertical pontilhada indica o PMC. Pode-se observar que os dois procedimentos apresentam o mesmo desempenho para praticamente todos os pontos da curva. Portanto, a forma mais simples foi adotada. Com relação à versão BX, os testes mostraram que o seu desempenho é afetado fortemente pela definição do vetor , sendo que a melhor alternativa é usá-lo como foi definido na equação (12). Observa-se que estes comportamentos se repetiram em todos os sistemas estudados.
O objetivo da figura 4 é o de comparar, para o sistema IEEE-57, o desempenho (número de iterações) das versões XB e BX do FCDRP entre si, bem como com as respectivas versões XB e BX do FC convencional, isto é, parametrizado por . Para isso, calcula-se os mesmos pontos da curva PV utilizando cada uma das versões, de forma a garantir que os mesmos possam ser corretamente comparados. Observa-se que cada ponto da curva PV corresponde a uma solução do problema de FC para um valor especificado de . Os pontos foram obtidos através da parametrização pela tensão da barra 31 (V31) com passo de 0,01 p.u usando a técnica de previsão trivial. Uma vez obtidos todos os estados, os respectivos valores dos demais parâmetros, e q31, foram usados como pontos estimados por suas respectivas versões, ou seja, usando a técnica de previsão trivial. Dessa forma garante-se que estes valores levarão o sistema de um mesmo estado inicial para o mesmo estado final. Conforme se pode ver nas figuras 4(a) e 4(b), ambas as versões do FCDRP não só obtiveram êxito na determinação do PMC, mas também permitem a obtenção de pontos além do PMC (isto é, pontos da parte inferior da curva PV) com um número pequeno de iterações, enquanto que as do FC convencional falham nesta tarefa.
Outro objetivo foi o de avaliar os efeitos da normalização dos mismatches (DP/V e DQ/V) da mesma forma que o usado no FCDR convencional (Monticelli et alii, 1990). Quando do uso deste procedimento, os vetores e também devem ser normalizados. Os resultados apresentados nas figuras 3 e 4 mostram que esta prática reduz o número de iterações próximo ao PMC e na parte inferior da curva PV, para ambas as versões. Também, pode-se observar na figura 4 que, no geral, a versão BX tem um melhor desempenho quando comparado à XB, conforme se verifica pelo número de iterações necessárias para a obtenção de cada ponto ao longo de toda curva.
Na figura 5 pode-se ver o desempenho da versão BX para o sistema IEEE-118. Verifica-se que as versões do FCDR e do Newton parametrizadas por V52, possibilitam mais uma vez o traçado da curva PV, e em particular, a determinação do PMC. Deve-se observar que a versão BX proposta, requer apenas a fatoração da matriz B" em cada ponto da curva em virtude da atualização do vetor , enquanto na versão XB a matriz B' é mantida constante. Isto representa uma vantagem evidente das versões propostas para o FCDR sobre os métodos de Newton parametrizados, que requerem a fatoração da matriz Jacobiana a cada iteração.
A tabela 2 ilustra alguns resultados para o sistema IEEE-14 (considerando os tap fixos e sem limites de geração de potência reativa), onde um ponto ( max) da curva muito próximo ao do PMC foi obtido usando como estado inicial o estado do ponto anterior. O desempenho (número de iterações) de diferentes versões podem ser comparadas, confirmando o melhor desempenho dos métodos parametrizados em relação ao dos convencionais. Da tabela também pode-se observar que os valores de max e Vcrítico encontrado pelas diferentes versões são praticamente os mesmos.
É fato conhecido que a solução de um conjunto de equações algébricas não-lineares com múltiplas soluções requer uma estimativa inicial suficientemente próxima de uma solução desejada. O objetivo da tabela 3 é o de explicitar a robustez da versão XB proposta, com relação à estimativa inicial.
Todos os pontos da tabela 3 foram obtidos com inicialização flat-start. O primeiro ponto da tabela corresponde à um ponto muito próximo do PMC ( = 1,8664 da tabela 1), enquanto que os demais correspondem à pontos da parte inferior da curva PV. Já a versão BX não apresentou o mesmo desempenho com relação a inicialização a partir do flat-start.
3.3 Efeito do chaveamento de parâmetro no método do FCDRP
A figura 6 e a tabela 4 mostram a relação entre o tamanho do passo e o número de iterações para as diferentes versões do FCDRP. As magnitude e o ângulo da tensão da barra 52 do sistema IEEE-118, para 4 pontos de operação (A a D) são mostrado na figura 6. Partindo do ponto A, é possível atingir o ponto B usando ou (método convencional), ou qk, ou Vk como parâmetro. O mesmo se aplica para a obtenção dos demais pontos (B, C, e D).
O número de iterações necessários para a obtenção de um ponto de operação a partir do ponto anterior (p.ex., o ponto C a partir do ponto B) estão apresentados na tabela 4. Para obter-se o ponto B a partir do ponto A, 14 iterações Pq e 13 QV são necessárias quando é o parâmetro, enquanto 6 iterações Pq e 6 QV e 10 iterações Pq e 10 QV são necessárias quando q52 e V52 são usados como parâmetros, respectivamente. Também pode-se observar que há uma diferença de 0,01 p.u. entre as tensões dos pontos de operação adjacentes. As diferenças entre os ângulos das tensões também são praticamente constantes. Contudo, as diferenças entre os fatores de carregamento são variáveis, e diminuem conforme o carregamento do sistema aumenta. Entretanto, a despeito da redução do passo, o número de iterações aumenta ao invés de diminuir, quando é usado como parâmetro. Assim, é possível deslocar-se de um ponto para outro mais eficientemente (com um número menor de iterações) através da mudança de parâmetro ao invés da redução do passo.
A parametrização implica num controle de passo implícito em . A redução de passo é necessária no caso em que o ''nariz'' associado com todos os parâmetros coincidem. De fato esta é uma situação atípica. Portanto, o uso de um limite predefinido do número de iterações para estabelecer a mudança de parâmetro durante o traçado da curva PV, asseguraria que o método com o melhor desempenho seja escolhido em cada parte da curva. Assim procedendo, como será apresentado a seguir, evitar-se-ia uma redução excessiva do tamanho do passo e um número excessivo de iterações durante o traçado da curva.
As figuras 7 e 8 mostram, respectivamente, os resultados de simulação obtidos com os métodos de Newton e do FCDRP modificados, para o sistema IEEE-118. Com o objetivo de se obter todos os pontos da curva PV com poucas iterações, os métodos são chaveados durante o traçado da curva. Isto é feito sempre que o número de iterações excede um limite predefinido. Se o número de iterações excede o limite, a magnitude ou o ângulo da tensão (Vk ou qk) com a maior taxa de decréscimo é usada como parâmetro. O processo de troca é realizado de forma que o método que apresenta o melhor desempenho seja escolhido em cada parte da curva. Na figura 8 apresenta-se uma tabela onde constam os limites de iterações predefinido, bem como do passo adotado para cada um dos parâmetros. Como se pode ver, com este procedimento o desempenho será bom para todos os pontos da curva. Um ponto importante a ressaltar, conforme mostram as figuras 8(d) e 8(e), é que o PMC não é tão susceptível ao tamanho do passo.
4 DESEMPENHO DA VERSÃO BX DO FCDRP PARAMETRIZADO PELA POTÊNCIA REATIVA DE UMA BARRA PV, QPV
O objetivo deste teste é o de mostrar a possibilidade do uso da potência reativa de uma barra PV como parâmetro de continuação para o FCDRP. Neste procedimento, somente é necessário especificar o passo no valor da potência reativa QPV (através do parâmetro m) de uma barra PV particular, que está dentro de seus limites de geração de potência reativa. O comportamento deste parâmetro é similar ao do parâmetro Vk, está continuamente crescendo à medida que o carregamento se aproxima do PMC. Assim, o valor de m de uma determinada barra PV poderia ser prefixado e a respectiva solução determinada. Observa-se que somente no sistema IEEE 118 foi possível o uso deste parâmetro para obter o PMC. Em geral, para os demais sistemas IEEE as barras PV atingem seus limites antes ou coincidentemente com o PMC. Isto mostra que as equações de potência reativas numa barra PV, nem sempre são os mais apropriados parâmetros para obter o PMC.
No início, qualquer potência reativa de uma barra PV que está dentro de seus limites poderia ser usada como um parâmetro. De qualquer forma, à medida que o sistema se aproxima do PMC, várias barras vão alcançando seus respectivos limites de geração de potência reativa. Portanto, é necessário monitorar as barras PV (isto já é feito usualmente durante o processo iterativo), e escolher aquela que sempre está dentro de seus limites.
A figura 9 mostra o desempenho do parâmetro m de três barras QPV do sistema IEEE-118 e o tradicional FCDR. Os pontos foram obtidos através da parametrização pela tensão da barra 52 (V52) com passo de 0,01 p.u usando a técnica de previsão trivial. Uma vez obtidos todos os estados, os respectivos valores dos parâmetros m (barras PV: 40, 42, 116) e o , foram usados como pontos estimados pela suas respectiva versões, ou seja, utilizou-se a técnica de previsão trivial com os pontos fornecidos pela versão BX parametrizada por tensão (V52). Como esperado, nas vizinhanças do PMC o método FCDRP encontra uma solução, enquanto a versão convencional, usando como parâmetro, falha. Observa-se que poder-se-ia traçar a curva partindo-se do valor de Q0 obtido no caso base, e, utilizando-se a técnica de previsão trivial, predefinindo um passo fixo para o parâmetro m. A figura 10 mostra o desempenho do sistema IEEE-118 adotando-se a potência reativa da barra 40, iniciando de Q0=0.2 p.u., considerando incrementos sucessivos iguais a 1 para o parâmetro m e os efeitos da normalização dos mismatches (DP/V e DQ/V).
Os métodos apresentados permitem o traçado completo das curvas PV através de simples modificações dos métodos convencionais de Newton e FCDR. Além disso, mantêm suas vantagens características. Dos resultados dos testes pode-se concluir que a versão BX do FCDR parametrizado é mais adequada para o traçado completo da curva PV. Por outro lado, a versão XB é mais robusta quando do cálculo de pontos individuais a partir do flat-start. Além disso, as versões parametrizadas do FCDR também permitem obter soluções na parte inferior da curva PV, contradizendo a noção de que as versões desacopladas sempre fornecem as soluções da parte superior da curva PV. O chaveamento entre as versões apresentadas, durante o traçado das curvas PV, possibilitam a determinação eficiente de todos os pontos da curva com um reduzido número de iterações.
A sugestão com relação às versões XB e BX apresentadas é de que elas devem ser usadas quando as convencionais falham; enquanto que as versões de Newton parametrizadas podem ser usadas como último recurso. O custo de chaveamento para versões alternativas quando as outras falham é mínimo no FCDRP quando comparado ao chaveamento para o método de Newton parametrizado.
Também se mostra a possibilidade de uso da potência reativa de uma barra PV como parâmetro da continuação do fluxo de carga desacoplado rápido.
Artigo submetido em 21/02/02
1a. Revisão em 29/05/02; 2a. Revisão em 23/10/02
Lesieutre B. C.; P. W. Sauer e M. A. Pai (1999). Existence of Solutions for the Network/Load Equations in Power Systems, IEEE Trans. on Circuits and Systems-I: Fundamental Theory and Applications, CAS-I-46 (8), pp. 1003-1011.
Força Tarefa Colapso de Tensão (1999). Critérios e Metodologias Estabelecidos no âmbito da Força - Tarefa Colapso de Tensão do GTAD / SCEL / GCOI para Estudos de Estabilidade de Tensão nos Sistemas Interligados Norte/Nordeste, Sul/Sudeste e Norte/Sul Brasileiros, XV SNPTEE, Foz do Iguaçu, PR, Brasil.
WSCC-Reactive Power Reserve Work Group (RRWG) (1998). Final Report, Voltage Stability Criteria, Undervoltage Load Shedding Strategy, and Reactive Power Reserve Monitoring Methodology, Maio.
Seydel, R. (1994). From Equilibrium to Chaos: Practical Bifurcation and Stability Analysis, 2a ed., Springer - Verlag, New York.
Kundur, P. (1993). Power System Stability and Control, McGraw - Hill, New York.
Van Cutsem, T. e C. Vournas (1998). Voltage Stability of Electric Power System, Kluwer Academic Publishers, Boston.
Ajjarapu, V. e C. Christy (1992). The Continuation Power Flow: a Tool for Steady State Voltage Stability Analysis, IEEE Trans. on Power Systems, Vol. 7, No 1, Fevereiro, pp. 416-423.
Cańizares, C. A.; F. L. Alvarado; C. L. DeMarco; I. Dobson e W. F. Long (1992). Point of Collapse Methods Applied to AC/DC Power Systems, IEEE Trans. on Power Systems, Vol. 7, No 2, Maio, pp. 673-683.
Bijwe, P. R. e R. S. Tare (1997). An Efficient Continuation Power Flow Method for Voltage Stability Analysis, Electric Machines and Power Systems, Vol. 25, pp.927-938.
Chiang, H. D.; A. Flueck; K. S. Shah e N. Balu (1995). CPFLOW: A Practical Tool for Tracing Power System Steady State Stationary Behavior Due to Load and Generation Variations, IEEE Trans. on Power Systems, Vol. 10, No 2, Maio, pp. 623-634.
Arya, L. D. e H. K. Verma (1996). A Method for Tracing PV Curve for Voltage Stability Analysis with Voltage Dependent Loads, Electric Machines and Power Systems, Vol. 24, pp.583-596.
Alves, D. A.; L. C. P. da Silva; C. A. Castro e V. F. da Costa (2000). Continuation Load Flow Method Parameterized by Power Losses, Proceedings of the IEEE Power Engineering Society Winter Meeting 2000, Singapore, Janeiro.
Draft of IEEE/PES (1999) Voltage Stability Report, Chapter 4, disponível em http://www.uwaterloo.ca/
Stott B. e O. Alsac (1974). Fast Decoupled Load Flow, IEEE Trans. on PAS, vol. 93, pp. 859-869.
Van Amerongen R. A. M. (1989). A General-Purpose Version of the Fast Decoupled Load Flow, IEEE Trans. on PAS, Vol. 4, No 2, pp. 760-770.
Monticelli A., A. Garcia, e O. R. Saavedra (1990). Fast Decoupled Load Flow : Hypothesis, Derivations and Testing, IEEE Trans. on Power Systems, Vol. 5, No 4, Novembro, pp. 1425-1431.
El-Abiad, A. H. (1983). Power Systems Analysis and Planning, McGraw - Hill, New York.
Peterson N.M. e W.S. Meyer (1971). Automatic Adjustment of Transformer and Phase-Shifter Taps in the Newton Power Flow, IEEE Trans. on PAS, Vol. PAS-90, No 1, Fevereiro, pp. 103-108.
Cańizares, C. A. (1995). Conditions for a Saddle-Node Bifurcations in AC/DC Power Systems, Electric Power & Energy Systems, Vol.17, No 1, pp. 61-68.
Zambroni de Souza, A. C.; C. A. Cańizares e V. H. Quintana (1997). New Techniques to Speed Up Voltage Collapse Computations Using Tangent Vectors. IEEE Trans. on Power Systems, Vol. 12, No 3, August, pp. 1380-1387.